A comprehensive analytical theory is presented for non-hysteretic RF SQUIDs operating in the adiabatic mode in the presence of large thermal fluctuations. When beta << 1 (beta = 2 pi LIc/Phi(0) is the hysteresis parameter, L is the SQUID inductance, I-c is the critical current of the Josephson junction, and Phi(0) is the flux quantum) the theory is applicable also for RF SQUIDs operating in the non-adiabatic mode. In contrast to previous theories in which the noise is treated perturbatively and which therefore are applicable only if the product beta Gamma << 1 (Gamma = 2 pi k(B)T/Phi(0)I(c) is the noise parameter, k(B) is the Boltzmann constant, and T is the absolute temperature)-the case of small thermal fluctuations-the present theory is valid for beta Gamma around unity or higher. In the limit beta Gamma --> 0 the theory reproduces the results of small thermal fluctuations theories. It has been found that in the presence of large thermal fluctuations the screening current in the SQUID inductance is suppressed by a factor that increases with increasing beta Gamma. Taking into account this new basic fact, all SQUID characteristics (output signal, transfer function, noise spectral density and energy sensitivity) have been recalculated and a good agreement with experimental data has been obtained. It has been also found that RF SQUIDs can be operated with substantially higher values of the inductance and of the noise parameter than DC SQUIDs. These two aspects, which are of particular importance at liquid nitrogen temperature, make high T-c RF SQUIDs very attractive.