Compiler error?

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

This is my simple example:

  PROGRAM COMPLEX

! IMPLICIT NONE ! REAL8 XIK COMPLEX16 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

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.

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.

VZ(3)=-CMPLX(.0d0,1.d0/XIK, KIND(.0d0))

But in this sense also

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

should return the single precision result.

Quoted from dpannhorst But in this sense also

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

should return the single precision result.

I agree, it should! And its in VZ(1) where the error occurs. VZ(3) is ok.

The default behaviour of complex is to return a single precision value (default real kind value, to be precise) so that it is compatible with the behaviour in Fortran 77 programs.

The default behaviour is equivalent to the following:

Z = CMPLX(SNGL(D1), SNGL(D2))

You only get double precision by adding the optional kind parameter:

Z = CMPLX(D1, D2, KIND=KIND(1.0D0)) Z = CMPLX(D1, D2, KIND(1.0D0)) Z = CMPLX(D1, D2, 2)

Some fortran compilers (e.g. gfortran) have an additional, non-standard intrinsic DCMPLX which does this for you:

Z = DMPLX(D1, D2)

Paul, it would be good to eliminate this inconsistency. I have noticed it as well, but usually I almost always the KIND argument, so I don't see it in my codes.

David.

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 (REAL4) KIND = 2 is double precision (REAL8) KIND = 3 is extended precision

Intel's FORTRAN help says:

KIND = 4 is single precision (REAL4) KIND = 8 is double precision (REAL8) KIND = 16 is extended precision (REAL*16)

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

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.

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.