Precision Guidelines

Following guidelines ensure that applications created using 2DECOMP&FFT can support single/double precision switch. Most of them are applicable to general Fortran programs.

• Most importantly, when defining new floating point variables, use real(mytype) instead of real. Here mytype is the proper KIND variable defined in the base module of 2DECOMP&FFT. Normally, do not use DOUBLE PRECISION (or REAL*8, which is bad style) directly to define variables unless required by external libraries.
• If complex variables are needed, defined them as complex(mytype).
• However, using mytype to compute memory displacement or file size is not portable and not reliable - compilers do not necessarily represent double precision as KIND 8. Use global variable mytype_bytes instead for such computations.
• Pay particular attention to the MPI data-type constants when calling MPI routines from applications. Use real_type instead of MPI_REAL or MPI_DOUBLE_PRECISION; use complex_type instead of MPI_COMPLEX or MPI_DOUBLE_COMPLEX. Both new types are defined in the base 2DECOMP&FFT module.
• Be very careful with explicit data type conversions in applications and make sure there is no loss of precision.
• Some Fortran intrinsic functions have built-in type conversion capability, which should be used wherever possible. For example:

        OMEGA = CMPLX(1.0, 0.0, kind=mytype)

• Some algorithms are sensitive to precisions. For example, in the following code if one does not use the "_mytype" part, a compiler might create a single-precision constant and then promote it to double before assigning the value to the variable, which may lose accuracy.

        REAL(mytype) :: two_pi
        two_pi = 6.28318530717958_mytype

• It is safest to define all floating point constants like above (that is what the Fortran standard says). However this makes code untidy in some situations. Use the technique at key places at least.
• If applications rely on external libraries, consider the single/double precision support of the library in question. If needed use the DOUBLE_PREC directives for code branching. For example in ACML FFT, work array sizes and subroutine names are dependent on precisions.

        #ifdef DOUBLE_PREC
           allocate(comm(max(3*nx+100, max(3*ny+100, 3*nz+100))))
           call dzfft(......
        #else
           allocate(comm(max(5*nx+100, max(5*ny+100, 5*nz+100))))
           call scfft(......
        #endif

• If algorithms allow, verify that machine accuracy can be obtained. For example, a forward FFT followed by a backward FFT plus proper scaling should recover the original input to machine accuracy.
• Some compilers provide flags to set the precision of floating point numbers at compile-time, such as -fdefault-real-8 in gfortran and -r8 in several other compilers. Do not rely on these flags as your code may become not portable.
• It may be necessary to define variables to satisfy the requirements of external libraries. For example, MPI_WTIME returns a double-precision result even if applications using it may be running in single precision mode.

If all the guidelines are properly followed, applications created will hopefully work for both single precision and double precision. Single precision is the default mode and switching to double can be done easily at compile time by passing the -DDOUBLE_PREC flag to the compiler.