Subject: SQRT in MD

From: [email-masked] (David van der Spoel)
Date: Sat, 7 Sep 1996 13:19:49 +0200 (MDT)

John Marelius <[email-masked]> wrote:

:The core of the problem is that square root calculations and also
:divisions take up to two orders of magnitude longer than multiplications 
:and additions to execute on most computers.

:Approaches by me and others involving integer-math, look-up tables,
:Taylor polynominals or Newton-Raphson iterations ended up being
:either slower or less accurate than the standard function.

The GROMACS MD code (http://rugmd0.chem.rug.nl/~gmx) includes
a newton rhapson iteration method for the 1.0/sqrt(x) function
that takes 5 floating point instructions, rather than 100
(two orders of magnitude ??) This enables our code
to vectorize or be pipelined on machines (compilers) that can not
inline library functions (e.g. the Portland Compiler on our 
parallel Intel-i860 box).

A description of the algorithm we use can be
found in the GROMACS User Manual, which is freely available from
our WWW page. In most cases it is considerably faster to use the
special instructions for 1.0/sqrt(x), e.g. the SGI R8000 chip has
such an instruction and the compiler detects 1.0/sqrt(x) constructs
in the code. Consequently the inner loop of our coulomb routine
(that calculates 1 coulomb interaction and forces, total 27 floating
point operations) executes in 12 clock cycles on a R8000, leading to
an inner loop performance of 168 MFlops, slightly more than
half the theoretical maximum of 300 (at 75 MHz). 

Regards,

David van der Spoel
---------------------------------------------------------
EMAIL:	[email-masked]
WWW:	http://rugmd0.chem.rug.nl/~spoel
PHONE:	31-50-3634327	FAX: 31-50-3634800
MAIL:	Nijenborgh 4, 9747 AG Groningen, The Netherlands.
---------------------------------------------------------