Program gives different results only on one computer - problem with FPU?

[ZlamalJakub]

I have problems with my notebook. My program gives on it different results. Sometimes correct, sometimes there appears NaN in different places of code and different variables (when I detect NaN and run this part of code again it returns correct result). On another computers are results consistent.

I think there is problem with FPU. Is there any thorough test of FPU? I tested RAM and it is OK.


Thanks for suggestions

Jakub

[jimdempseyatthecove]

Jakub,

The first thing to do is to enable runtime diagnostics for the Debug configuration to make access of uninitialized variables checks. This is the most likely cause for errors like this.

You cannot assume that variables are pre-initialized to 0 prior to use.

Also, use IMPLICIT NONE
This will catch a lot of typographical errors that result in new (unintended) variables to be created (and thus unintendedly not initialized).

Jim Dempsey

[ZlamalJakub]

My code is compiled with CVF 6.6C (it is old version of my program) and I think there I cannot use runtime diagnostics.

Problem occurs in simple routine calculating field derivatives like

subroutine calculate(ix,iy,iz,type,Value)
use FieldModule ! defined NX,NY,NZ,Field
integer*4 ix,iy,iz,type
real*8 value

ipass=0
10 continue
ipass=ipass+1
select case (type)
...

case (3)
if (iz < NZ) then
Value=Field(ix,iy,iz+1)-Field(ix,iy,iz)
endif
...



NZ (field bounds) and Field are declared in module.

I check every array bound, ix,iy,iz varies in bounds of field and are allways defined.

I hope all in this routine is initialised to correct value


when I put at the end of code

if (isnan(Value) .AND. ipass==1) then
goto 10
endif

in pass 2 is not Value=NaN.

[DavidWhite]

If the routine is not called with Type = 3, is Value set anywhere else. Likely that it has not been initialized - cannot asuume variables are initialized, so need to initialize every variable.

David

[mecej4]

Given how little of the code has been shown, one can only guess as to the reasons, and here is one such guess.

If you added the "if (isnan..." test assuming that Value would have the value NaN if it had not been assigned a value before, you are mistaken. An uninitialized variable may have any value.

That is why one of the usual symptoms of the presence of uninitialized variable is the obtaining of different results in different runs with the same input data.

[jimdempseyatthecove]

Not seeing all of the code it is hard to diagnose the problem

...
case (3)
if (iz < NZ) then
Value=Field(ix,iy,iz+1)-Field(ix,iy,iz)
endif
...

Assuming ix, iy, iz (args passed into routine)are in bounds of the array Field
Then when iz == NZ the variable Value is not set (remains prior state of Value) in the above case(3)

Value may or may not have been written to elsewhere in your routing (code not shown).
If Value was not written to in your routine. Then it may contain

a) Valid prior number from caller
b) Uninitialized data
b.1) Uninitialized to arbitrary butvalid number ( not NaN)
b.2) Uninitialized to arbitrary but invalid number (NaN)

Also your select case may be presented with a type that is not in one of your cases

Jim

[ZlamalJakub]

Thank for all advices I will try and will report results it may be useful for others.

Jakub

[ZlamalJakub]

[fortran]recursive subroutine CalcDerivative(ix,iy,iz,iDiff,Value)
! calculates derivative of field
use FieldData
use FieldInterpol
	real*8 Value
	real*8 Val(3)

	Value=0.D0
	ipass=0
10 continue
	ipass=ipass+1
	select case (iDiff)
	case (D_POT)
		Value=Field(ix,iy,iz)
	case (D_POTX)
		! derivative in x direction
		if (ix == 1) then 
			! forward
			if (mirror(1)) then	! is mirrored about axis
				Value=0.D0
			else
				Value=(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix+1,iy,iz)-Field(ix+2,iy,iz))/(2.D0)
			endif
		elseif(ix ==  Header.nx) then
			! backward
			Value=-(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix-1,iy,iz)-Field(ix-2,iy,iz))/(2.D0)
		else
			! normal
			Value=(Field(ix+1,iy,iz)-Field(ix-1,iy,iz))/(2.D0)
		endif
	case (D_POTY)
		! derivative in y direction
		if (iy == 1) then 
			! forward
			if (mirror(2)) then	! is mirrored about axis
				Value=0.D0
			else
				Value=(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix,iy+1,iz)-Field(ix,iy+2,iz))/(2.D0)
			endif
		elseif(iy ==  Header.ny) then
			! backward
			Value=-(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix,iy-1,iz)-Field(ix,iy-2,iz))/(2.D0)
		else
			! normal
			Value=(Field(ix,iy+1,iz)-Field(ix,iy-1,iz))/(2.D0)
		endif
	case (D_POTZ)
		! derivative in z direction
		if (iz == 1) then 
			! forward
			if (mirror(3)) then	! is mirrored about axis
				Value=0.D0
			else
				Value=(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix,iy,iz+1)-Field(ix,iy,iz+2))/(2.D0)
			endif
		elseif(iz ==  Header.nz) then
			! backward
			Value=-(-3.D0*Field(ix,iy,iz)+4.D0*Field(ix,iy,iz-1)-Field(ix,iy,iz-2))/(2.D0)
		else
			! normal
			Value=(Field(ix,iy,iz+1)-Field(ix,iy,iz-1))/(2.D0)
		endif
	case (D_POTXY)
		! derivative in x direction of derivative in y direction
		if (ix == 1) then 
			! forward
			if (mirror(1)) then	! is mirrored about axis
				Value=0.D0
			else
				call CalcDerivative(ix,iy,iz,D_POTY,Val(1))
				call CalcDerivative(ix+1,iy,iz,D_POTY,Val(2))
				call CalcDerivative(ix+2,iy,iz,D_POTY,Val(3))
				Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
			endif
		elseif(ix ==  Header.nx) then
			! backward
			call CalcDerivative(ix,iy,iz,D_POTY,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTY,Val(2))
			call CalcDerivative(ix-2,iy,iz,D_POTY,Val(3))
			Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
		else
			! normal
			call CalcDerivative(ix+1,iy,iz,D_POTY,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTY,Val(2))
			Value=(Val(1)-Val(2))/(2.D0)
		endif
	case (D_POTXZ)
		! derivative in x direction of derivative in z direction
		if (ix == 1) then 
			! forward
			if (mirror(1)) then	! is mirrored about axis
				Value=0.D0
			else
				call CalcDerivative(ix,iy,iz,D_POTZ,Val(1))
				call CalcDerivative(ix+1,iy,iz,D_POTZ,Val(2))
				call CalcDerivative(ix+2,iy,iz,D_POTZ,Val(3))
				Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
			endif
		elseif(ix ==  Header.nx) then
			! backward
			call CalcDerivative(ix,iy,iz,D_POTZ,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTZ,Val(2))
			call CalcDerivative(ix-2,iy,iz,D_POTZ,Val(3))
			Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
		else
			! normal
			call CalcDerivative(ix+1,iy,iz,D_POTZ,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTZ,Val(2))
			Value=(Val(1)-Val(2))/(2.D0)
		endif
	case (D_POTYZ)
		! derivative in y direction of derivative in z direction
		if (iy == 1) then 
			! forward
			if (mirror(2)) then	! is mirrored about axis
				Value=0.D0
			else
				call CalcDerivative(ix,iy,iz,D_POTZ,Val(1))
				call CalcDerivative(ix,iy+1,iz,D_POTZ,Val(2))
				call CalcDerivative(ix,iy+2,iz,D_POTZ,Val(3))
				Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
			endif
		elseif(iy ==  Header.ny) then
			! backward
			call CalcDerivative(ix,iy,iz,D_POTZ,Val(1))
			call CalcDerivative(ix,iy-1,iz,D_POTZ,Val(2))
			call CalcDerivative(ix,iy-2,iz,D_POTZ,Val(3))
			Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
		else
			! normal
			call CalcDerivative(ix,iy+1,iz,D_POTZ,Val(1))
			call CalcDerivative(ix,iy-1,iz,D_POTZ,Val(2))
			Value=(Val(1)-Val(2))/(2.D0)
		endif
	case (D_POTXYZ)
		! derivative in x direction of derivative in yz direction
		if (ix == 1) then 
			! forward
			if (mirror(1)) then	! is mirrored about axis
				Value=0.D0
			else
				call CalcDerivative(ix,iy,iz,D_POTYZ,Val(1))
				call CalcDerivative(ix+1,iy,iz,D_POTYZ,Val(2))
				call CalcDerivative(ix+2,iy,iz,D_POTYZ,Val(3))
				Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
			endif
		elseif(ix ==  Header.nx) then
			! backward
			call CalcDerivative(ix,iy,iz,D_POTYZ,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTYZ,Val(2))
			call CalcDerivative(ix-2,iy,iz,D_POTYZ,Val(3))
			Value=(-3.D0*Val(1)+4.D0*Val(2)-Val(3))/(2.D0)
		else
			! normal
			call CalcDerivative(ix+1,iy,iz,D_POTYZ,Val(1))
			call CalcDerivative(ix-1,iy,iz,D_POTYZ,Val(2))
			Value=(Val(1)-Val(2))/(2.D0)
		endif
	end select
	if (isnan(Value) .AND. ipass==1)then
		call ErrMessage("Isnan"C,"isnan"C,0)
		goto 10
	endif
	return
end		[/fortran]

This is code of my routine

Header,n,Header.ny, Header.nz is size of array Field.
Values ix,iy,iz are allways in bounds of field.
Vaue idiff is allways one of selected case.

I tried and checked everything and problem still persists. sometimes (every run at different ix,iy,iz coordinates) I obtain isnan(Value). When I trace second pass in debugger, value is correct.

I think problem is in my FPU.

Jakub

[peterklaver]

Have you tried compiling with fpe:0 and then stepping through the routine? I'm thinking in this case the debugger would break with an unhandled exception instead of just assigning NaN to Value, and then you could at least see where it is happening.

[Steven_L_Intel1]

I suggest that the FPU is the least likely culprit here.

[jimdempseyatthecove]

Jakub,

>>Values ix,iy,iz are allways in bounds of field

Have you inserted an assert to verify the bounds of ix,iy,iz?

Enable FPP (Fortran preprocessor)

Add in front of your subroutine a function that explicitly tests the array subscripts for array Field

function debugField(ix, iy, iz)
use FieldData
use FieldInterpol
real*8 :: debugField
integer :: ix, iy, iz
if((ix < 1) .or. (iy < 1) .or. (iz < 1) &
& .or. (ix > Header.nx) .or. (iy > Header.ny) .or. (iz > Header.nz)) then
Call ErrMessage("Index out of bounds")
endif
debugField = Field(ix,iy,iz)
end function debugField

Following that add an FPP directive that will change all array references to Fieldto function references to debugField.

#define Field debugField

Set your break points on all Call ErrMessage

Run your test something may show up

Note,

On your recursive calls you are calling with various indexes with +1, +2, -1 and -2 displacements from the original ix, iy or iz indexes. One of these may extend out of range of the indexes.

Although the runtime check for index out of bounds should have caught this, the above changes will catch the error when compiled in release mode with full optimizations.

If this doesn't expose the problem

Then after every "Value=(expression)" add

if(isnan(Value)) call ErrMessage("Value is NaN")

Then place break point inside ErrMessage, on break, walk up the stack one level and look at the location in your case selection plus any residual data. In release mode with full optimizationthis may be difficult. You can add additional code to help trap the bug.

Jim Dempsey

[ZlamalJakub]

I will try your suggestions.

Problem occurs in the middle of the array so index +-2 cannot cause problems (ix is allways >=2 <=Header.nx-2)

Problem with NaN values occurs only on my notebook other computers I am using have not problems.


I will report my investigation.

Jakub