forums.silverfrost.com

dpannhorst · Joined: 29 Aug 2005 Posts: 165 Location: Berlin, Germany

I have detected a problem in calculation of complex values, which give wrong results.

This is my simple example:

PROGRAM COMPLEX
!
IMPLICIT NONE
!
REAL*8 XIK
COMPLEX*16 VZ(4)
!-----------------------------------------------------------------------
XIK=3.997186510534627D-003
VZ(1)=+CMPLX(.0d0,1.d0/XIK)
VZ(2)=+DCMPLX(.0d0,1.d0/XIK)
VZ(3)=-CMPLX(.0d0,1.d0/XIK)
VZ(4)=-DCMPLX(.0d0,1.d0/XIK)
!
OPEN(20,FILE='RESULT')
WRITE(20,*)'XIK=',XIK
WRITE(20,*)'VZ(1)=',VZ(1)
WRITE(20,*)'VZ(2)=',VZ(2)
WRITE(20,*)'VZ(3)=',VZ(3),' <---'
WRITE(20,*)'VZ(4)=',VZ(4)
CLOSE(20)
!
END

and these are the results of the four calculations:

XIK= 3.997186510535E-03
VZ(1)= (0.0000000000,250.175966862)
VZ(2)= (0.0000000000,250.175966862)
VZ(3)= (0.0000000000,-250.175964355) <---
VZ(4)= (0.0000000000,-250.175966862)

As you can see the 3rd calculation gives a wrong result.

Any ideas?

Detlef Pannhorst

PaulLaidler · Posted: Fri Apr 27, 2012 9:16 am Post subject:

There is an inconsistency here in the compiler generated code.

vz(1) is calculated directly whilst vz(3) uses a temporary intermediate store which is single precision when presumably it should be double precision.

Either way the results ought to be consistent.

I will log this for investigation.

PaulLaidler · Posted: Sun Apr 29, 2012 8:44 am Post subject:

I have looked at this again and I now think that it is OK.

The third argument of Fortran 90 intrinsic CMPLX is optional and, if you leave it out, CMPLX returns a value of the default kind which is single precision.

If you get more accuracy, as in VZ(1), then that is a bonus. In other words, VZ(1) and VZ(3) are effectively equal for the precision requested.

The following line gives you the result you require.

dpannhorst · Joined: 29 Aug 2005 Posts: 165 Location: Berlin, Germany

But in this sense also

VZ(1)=+cmplx(.0do,1.d0/XIK)

should return the single precision result.

davidb · Joined: 17 Jul 2009 Posts: 560 Location: UK

dpannhorst · Joined: 29 Aug 2005 Posts: 165 Location: Berlin, Germany

I hav large problems understanding the KIND parameter.

For example for the use of REAL values:

Silverfrost's FORTRAN help says:

KIND = 1 is single precision (REAL*4)
KIND = 2 is double precision (REAL* Cool

KIND = 3 is extended precision

Intel's FORTRAN help says:

KIND = 4 is single precision (REAL*4)
KIND = 8 is double precision (REAL* Cool

KIND = 16 is extended precision (REAL*16)

So it seems KIND is not portable??? FORTRAN is not FORTRAN???

davidb · Joined: 17 Jul 2009 Posts: 560 Location: UK

The actual KIND values do vary from compiler to compiler.

If you want the double precison kind number you can use the kind intrinsic with a double precision constant argument:

INTEGER, PARAMETER :: DP = KIND(1.0D0)
REAL(KIND=DP) :: A

Of course you can still say:

DOUBLE PRECISION A

Or you can specify what range and precision you require as your working precision:

INTEGER, PARAMETER :: WP = SELECTED_REAL_KIND(13, 307)
REAL(KIND=WP) :: A

Both of these will give you 8 byte reals with Silverfrost, Intel and many others, equivalent to REAL*8.

This is the only way to ensure code is portable.
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

PaulLaidler · Posted: Sun Apr 29, 2012 4:43 pm Post subject:

As I understand it, VZ(1) has been assigned the single precision result. In other words the result is correct to the expected number of decimal places (and more as it happens). In any calculation the truncation error will vary with the context, with the compiler and perhaps even with the processor. There is no point in the compiler writer introducing a deliberately different value for the truncation error.