non ieee floating point exception 2

[mtamimi]

Hello
I am writing a monte carlo simulation fortran 90 code. The code reads data from a file then uses these data to simulate particle transport in a meduim using random number generator. I am using fortran power station 4 (old compiler) given to me by my professor so help is limited.
The code runs fine with low number of histories. But for large number of histories it terminates without any explainations why.
The code was crashing for several reasons:
1) Flaws in the logic where wrong data is looked up --> I fixed.
2) Asking for array element that does not exist. --> fixed.
3) Rounding values. Was easily detected and fixed.

I looked for ways to trap the floating point exceptions and debug the code using FPS and sure enought it has some help. For floating point error traping and debuging It instructed me the code below:
! Exception handler routine hand_fpe
FUNCTION hand_fpe (signum, excnum)
!MS$ATTRIBUTES C :: hand_fpe
USE MSFLIB
INTEGER(2) signum, excnum
WRITE(*,*) 'In signal handler for SIG$FPE'
WRITE(*,*) 'signum = ', signum
WRITE(*,*) 'exception = ', excnum
SELECT CASE(excnum)
CASE(FPE$INVALID )
STOP ' Floating point exception: Invalid number'
CASE( FPE$DENORMAL )
STOP ' Floating point exception: Denormalized number'
CASE( FPE$ZERODIVIDE )
STOP ' Floating point exception: Zero divide'
CASE( FPE$OVERFLOW )
STOP ' Floating point exception: Overflow'
CASE( FPE$UNDERFLOW )
STOP ' Floating point exception: Underflow'
CASE( FPE$INEXACT )
STOP ' Floating point exception: Inexact precision'
CASE DEFAULT
STOP ' Floating point exception: Non-IEEE type'
END SELECT
hand_fpe = 1
END
to use this code I needed to set the error handeling handle early on the program
iret = SIGNALQQ(SIG$FPE, hand_fpe)
using this code intreduced the default error ' Floating point exception: Non-IEEE type' I can seem to know why. tracing the error lead me to the following code:
Do jj = 1, UsableNumGroups, 1

Read(FreeFile1,231) GXS(jj)
231 Format(66X, E13.5)

End Do
The only way to clear this error if I declare the GXS(jj) as double precision-Real(8). I am trying to read a set of numbers that look like " 1.96607E-03" arranged in columns. I can read them fine if I exchange the E13.5 with a13 and GXS with a character variable.


Does any one know what I am doing wrong and why???
Any help well be greatly appreciated.

Mat M Al-Tamimi
mtamimi@netzero.net

 

[xwb]

Assuming you are doing something like

program main
integer, parameter:: gxsmax = 20
real, dimension(gxsmax):: gxs
integer, parameter:: datafile=30

open (datafile,...)
read (datafile, '(5E13.5)') (gxs(ii), ii = 1, gxsmax)
close (datafile)

You can declare GXS as either real(4) or real(8). If it is real(4), the use E format. If it is real(8) use D format. Alternatively use G format which is for all forms of numbers.

 

[gummibaer]

Well, it would be better if you show us the code that produces the error together with a section of your inputfile rather than the one you are trying to handle the error with. In fact, I would forget about an error trap like you coded. When you post your code and the section of your data, include them in code-tags, that is, place the tags first by using the above icon and copy and paste your code and data between the tags. Otherwise html would not show the format properly (dropping the leading blanks).

I would take a guess, that your input format is incompatible to the dataformat in your file. Note: reading numerical data as character variables is no solution because you will encounter all sorts of errors trying to use these in any type of math.

Try to have a printout of GXS - depending on the number of elements maybe only some of these - to make sure you have got proper values.

And try to use list directed input, so you do not have to be too exacting with your format.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[mtamimi]

Hello
Thanks for the replies. The code is over 5000 lines. I an pasting the parts that I think is causing the error message.
In put file:

################################################################################
Material = Dry_Air
Material-Dry_Air|Probability [fraction (0-1)],Heat XS units (cGy - cm^2), GXS data units (1/cm)
GRP#|EBinMidPoint|AbsorptionProb|AbsorptionHeat|ScatteringProb|ScatteringHeat|TotalXs|TotalHeat
001| 1.500E+04| 0.830118| 3.09787E-10| 0.169882| 1.27139E-12| 1.96607E-03| 3.11058E-10
002| 2.500E+04| 0.542463| 1.03536E-10| 0.457537| 3.21679E-12| 5.86184E-04| 1.06753E-10
003| 3.750E+04| 0.267581| 4.18371E-11| 0.732419| 6.85515E-12| 3.20345E-04| 4.86921E-11
004| 5.250E+04| 0.117951| 1.96932E-11| 0.882049| 1.19526E-11| 2.42317E-04| 3.16454E-11
005| 6.500E+04| 0.065320| 1.21952E-11| 0.934680| 1.67657E-11| 2.17489E-04| 2.89608E-11
006| 8.500E+04| 0.030157| 6.57637E-12| 0.969843| 2.65001E-11| 1.96406E-04| 3.30751E-11
007| 1.250E+05| 0.010927| 3.00675E-12| 0.989073| 4.47206E-11| 1.76119E-04| 4.77274E-11
008| 1.750E+05| 0.003866| 1.36527E-12| 0.996134| 7.09435E-11| 1.56964E-04| 7.23040E-11
009| 2.500E+05| 0.001505| 6.56073E-13| 0.998495| 1.09652E-10| 1.39247E-04| 1.10308E-10
010| 3.500E+05| 0.000578| 3.18925E-13| 0.999422| 1.62472E-10| 1.22379E-04| 1.62791E-10
011| 4.250E+05| 0.000339| 2.15206E-13| 0.999661| 2.02150E-10| 1.13383E-04| 2.02364E-10
012| 4.800E+05| 0.000256| 1.73915E-13| 0.999744| 2.28136E-10| 1.07438E-04| 2.28310E-10
:
:
:
This data is being read by the a subroutine called from main program:
Main program:
Program Pixel8SnglBeamSimulation
Use PORTLIB
Use MSFLIB
Implicit None

!*) Establish the name of the exception handler as the function to be
! invoked if an exception happens {hand_fpe} is attached below.
!######################################################################!
INTERFACE
FUNCTION hand_fpe (sigid, except)
!MS$ATTRIBUTES C :: hand_fpe
INTEGER(4) hand_fpe
INTEGER(2) sigid, except
END FUNCTION
END INTERFACE

$Debug
Character(Len=200) TransverseXY, SagitalYZ, CoronalXZ
Character(Len=200) CAXPDDDir, PlanarResultDir, ProblemDir
Character(Len=200) BeamDir,VoxelDoseDir, Voxel2DoseDir
Character(Len=120) InfileDirLoc, InFileName, ProblemTitle
Character(Len=120) GroupName, GroupDir, DelFirstCharIfBlank
Character(Len=120) LUString, GetStringValue, BannerMessage
Character(Len=120) Message
Character(Len=24) Stamp
Character(Len=40) ProblemType, Materials)
Character(Len=15) MaterialMapType
Character(Len=10) JawsType
Character(Len=8) FieldType
Real AirTotalXs), AirSigmaTotal
Real ScattHeat,, TotalXs,, AbsorptionHeat,
Real ScattProb,, DownScatterGTGProb,, MUData,:,
Real VoxelDose,:,, Voxel2Dose,:,
Real VoxelDoseStdDev,:,
Real MaxSigma, SigmaTotal
:
:
Integers.........................
:
:
Allocatable..............

!Set the floating point error handler:
!----------------------------------------
iret = SIGNALQQ(SIG$FPE, hand_fpe)


!15) Read Air total cross section only !
!######################################################################!
LUString = "Material-" // Trim(Materials(0)) // Char(124) // &
& "Probability [fraction (0-1)],Heat XS units (cGy - cm^2), GXS data units (1/cm)"

Call ReadGXSDataOnly(InfileDirLoc,InfileName, LUString, &
& SpectrumNumGrps, AirTotalXs, nFile, WriteFlag)

:
:
:
End Program Pixel8SnglBeamSimulation

Subroutine ReadGXSDataOnly(FileDirectory,FileName,SearchTag,&
& UsableNumGroups, GXS, nFile, WriteFlag)
!
!Purpose: read the GXS data from data file
!

Use PORTLIB
Use MSFLIB
Implicit None

Character*(*) FileDirectory,FileName,SearchTag
Character(Len=1) Acomplished
Character(Len=150) Text

Integer FreeFile1, UsableNumGroups, jj, nFile

Real GXS(UsableNumGroups)

Logical LOpen, LStatus, WriteFlag

!1) Change directory to "FileDirectory":
!-------------------------------
LStatus = CHANGEDIRQQ(FileDirectory)

!2) make an inquiry about "FileName"
!-----------------------------------
Inquire(File= FileName,Opened=LOpen, Number=FreeFile1)

!3)Check inquiry results
!-----------------------
If(.Not. LOpen) Then
!File Not open --> open it for read
Call GetUnit(FreeFile1)
Open(Unit=FreeFile1,FILE= FileName,ACTION='Read')

Else
!file is open --> Rewind it
Rewind(FreeFile1)

End If

!4) Search file for the begining of needed data:
!-----------------------------------------------
Acomplished = 'N'

LoopStartData: Do While(Acomplished .Eq. 'N')

Read(FreeFile1,'(A150)') Text

If(Trim(Text) .Eq. Trim(SearchTag)) Then

!Skip a line
Read(FreeFile1,'(A150)') Text
!Start reading from the next line on

Do jj = 1, UsableNumGroups, 1

Read(FreeFile1,231) GXS(jj)
231 Format(66X, E13.5)

End Do

Acomplished = 'Y'

End If

End Do LoopStartData

!5) Close the file:
!------------------
Close(FreeFile1)

!6) Write to progress file:
!--------------------------
If(writeFlag) Then

Write(nFile,*)"Air cross section data:"

!==================================
!To check if data was read right:-
!----------------------------------
!DO jj =1, UsableNumGroups, 1
! write(nFile,231) GXS(jj)
!END DO

Write(nFile,*)""

End If

End Subroutine ReadGXSDataOnly
============================================================================
After I posted the question yesterday I seperated these parts from the entire program to trouble shoot it. It ran fine with now issues. The only difference between the two is the way I look for the data.
In the code that I have a problem with I use the following to look for the begining of the data

Material-Dry_Air|Probability [fraction (0-1)],Heat XS units (cGy - cm^2), GXS data units (1/cm)

The code I used yesterday I used "Material = Dry_Air" to look for the start of the data.

Both strings search seem to function correct. I get to the data but one well read it and the other will not. I use the same method to read more data therefore the program terminate at this point.
I am not sure at this point. any thought or do I need to post more of the code.


Mat M Al-Tamimi
mtamimi@netzero.net

 

[xwb]

66X ends up right in the middle of a number.

If you open up your file in notepad, switch off line wrap, click on the column where your data starts, look at the bottom of notepad, it will tell you which column it is in. On the the data you pasted, it is more like 62X and the format is E12.5

 

[gummibaer]

Material-Dry_Air|Probability [fraction (0-1)],Heat XS units (cGy - cm^2), GXS data units (1/cm)

The code I used yesterday I used "Material = Dry_Air" to look for the start of the data.

When I understand your code properly, then you drop one line after you have found Text to be equal to your SearchTag. Looking at your file I'd say you need two lines to drop if you looked for 'Material = Dry_Air', only one if you looked for 'Material-Dry_Air'. So if you exchange your SearchTag you should adjust the number of lines to drop too.

On the the data you pasted, it is more like 62X and the format is E12.5

No, you cannot tell. mtamimi choose not to include his code and data in code tags, so the blanks are not shown properly. html drops blanks and spoils the character count.

I would recommend list directed input and forget about the format altogether. If his data are in the 7th column the statement would be

...
integer iDummy
real rDummy
...
...
read (FreeFile1,*) iDummy, (rDummy, j = 1,5), GXS(jj)
...

Maybe you will have to get rid of the '|' delimiters first by preprocessing the inputfile and replace all '|' by ','. This can easyly be done in any texteditor by search and replace function.

And, at least while you debug your code, you should set writeFlag to .true. and remove the '!' from the output statements. Once you are sure you get what you want you can comment them off again.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[mtamimi]

Hi
First of all thank you guys for your input regrading my question. I do not Think I can thank you enough for your comments. I have been working for a while without help and a friend of mine suggest that I use your forum. I am not sure what is "code tags" or how to use it. And this is my first time asking a question in forum and I did not have my code at a url to point to it. So please for give me for my short comings.
In my format statment -66x - well position me right after the "|" in the data line. This is according to notepad++. Then two spaces then "1.96607E-03" before I get to the next "|". That is how I get E13.5. The funny thing is that format statement work fine if I do not have the floating point handler or the $Debug statment in code.
In monte carlo simulation the code generate a random number to make choices regarding a physical property. My code run fine when I loop through the process <= 200,000 (histories). The uncertainty in the answer well be large because of that. To reduce it I need to run it more. That is when I run into problems where the program terminates without giving a reason. That is why I added the floating point handler to trap that type of error. This caused the non ieee floating point exception erro as the fpe handler generated. and doe not allow the code to excute. Therefore, I do not think I know why or where the problem exist at high number of histories. I suspect that at large histories a random number gets generated that produces an error that well cause termination.
Is their a better way to trap that error than the way I added per FPS 4.0 recommendation?
I assumed that this error is FPE since its the one most post says it well terminate the program. Any thoughts??????? "Keep in mind that I really do not know the type of error generated."

Mat M. Al-Tamimi
mtamimi@netzero.net

 

[gummibaer]

I'm in doubt that your handler will do you any good. It needs a numerical errorcode to work - and if your system issues an error code it would very sure issue an errormessage on program termination.

So I take it your code works this way:
- you read in a big number of GXS-data
- controlled by a random number you select any of these GXS as some parameter for the next execution of a loop
- this works for quite some time and then your prog terminates without telling you why.

As a first step, before the loop starts computing data, write out your random number and the GXS together with all parameters that are modified between loop executions to your screen. Of course, your prog is now hurling figures at you, but you will see what the values are, when your prog quits execution. Check them all, if they are in the range as expected. If all looks okay, we will continue then.

BTW: Might be a good idea to have you rprog write a message to your screen, when the normal end of the execution is reached. Just to see that the exit your code takes really is an exception.


Norbert



The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[xwb]

Code tags are enclosing the code in

 and

 

[xwb]

Not sure how you would use your SIGFPE routine. I'll check my FPS4.0 and get back to you on this one. It is getting the error code that is the difficult bit. You could try the old fashioned err= for clues as to where the error might lie.

...
logical accomplished
...
accomplished = .false.
do jj = 1, UsableNumGroups, 1

    ! read the line into a buffer so we can print it
    read(FreeFile1, '(A150)') Text

    ! read the data from the buffer
    read(Text,231, err=999) GXS(jj)
231 Format(66X, E13.5)

end do
accomplished = .true.
999 if (.not. accomplished) then
   print *, 'Stopped after reading ', jj, ' lines in usablenumgroups'
   print *, Text
end if

 

[xwb]

In the hand_fpe, in the default case, print out what excnum you are getting.

If you look in MSFLP.F90 in the include directory, there are more codes listed like FPE$UNEMULATED, FPE$STACKUNDERFLOW etc.

If you are using the read (..err=) method, signalqq should not be called. In the error handler, (999 if statement) call GetStatusFPQQ(status) and print out the error status. These codes are listed in MSFLP.F90 after the FPE codes.

 

[mtamimi]

Hi
Thanks for your input.
1) In response to "gummibaer (Programmer)"
I did use the print statements as you mentioned earlier. I have so many variables that lookup and compute each time. And I need to run the same senario several millions of times. Writing to file (choose seperate files since one well be very large and impossible to open) was easy but generated a lot of files and slowed down the process & system. writing to the screen was better but occupied the PC with that program only. I was looking for a better way to find the problem.

2) In response to "xwb (Programmer)"
I do use the "ERR=" in other open statements in the program but not at that particular open statement. When I run my code i get the following:
In signal handler for SIG$FPE
signum = 8
exception = 138
Floating point exception: Non-IEEE type

Which mean that the case selected was the default. I do not know exactly what signum & exception are or what the default case is. According the FPS 4.0 help most these code are in the "MSFLIB.F90" but I could not find any thing as shown above. I do not have the "MSFLP.F90" in my include folder is it the same as "MSFLIB.F90" ???

I modified the error handler to include the Call GetStatusFPQQ(status), The program crashes and terminates. Is their a way to know why??

Function hand_fpe (signum, excnum)
!
!Purpose: Exception handler routine hand_fpe
!

!MS$ATTRIBUTES C :: hand_fpe
USE MSFLIB

Integer(2)  signum, excnum, status

	Call GetStatusFPQQ(status)

    WRITE(*,*) 'In signal handler for SIG$FPE'
    WRITE(*,*) 'signum = ', signum
    WRITE(*,*) 'exception = ', excnum
	WRITE(*,*) 'status  = ', status

	!Select Error Case
    Select Case(excnum)

      Case(FPE$INVALID)
        STOP ' Floating point exception: Invalid number'

	  Case(FPE$DENORMAL)
        STOP ' Floating point exception: Denormalized number'

	  Case(FPE$ZERODIVIDE)
        STOP ' Floating point exception: Zero divide'

	  Case(FPE$OVERFLOW)
        STOP ' Floating point exception: Overflow'

	  Case(FPE$UNDERFLOW)
        STOP ' Floating point exception: Underflow'

	  Case(FPE$INEXACT)
        STOP ' Floating point exception: Inexact precision'

	  Case default
        STOP ' Floating point exception: Non-IEEE type'

	End Select

    hand_fpe = 1

End Function hand_fpe

Thanks in advance for all the help.

 

[gummibaer]

Maybe I should explain the strategy I would use in your case.

(1) Track the values that cause the problem.
Since your prog performs quite a number of loops before it stops execution I would estimate a problem in the numerics. In my environment for instance, my code just stops executing when out of range numericals occur, like arcsin or arccos of arguments exceeding + 1.0 or - 1.0.

(2) Check data processing
Once the set of parameters that causes your problem is isolated, just have these data selected as the first by setting the appropriate random number nanually and printout intermediate values - or use an online debugger if you have one available.

I fear, there is no other way, tedious as it may be, to find out which data cause your problems.

Just while typing it occurs to me:
Might be a good idea to check your code, if you use any functions that allow a certain range of parameters only, like asin() or acos() etc. and check the arguments before calling this function. Something like

...
if (arg .lt. (-1.0) .or. arg .gt. 1.0) then
    write (*,*) 'got this bugger, random number is ',iRandom
    write (*,*) 'Arg = ',arg
!   print out anything that might help you to identify where this value came from
    return
endif
res = asin (arg)
...

If you do not have any problem here, I would say you will have to go the hard way as above and find out the values that trip your program.

Norbert



The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[xwb]

Signum = 8 = SIG$FPE. Basically, you told it to use the error handler when it got a signum of SIG$FPE. Just in case several signal handlers use the same routine, it tells which signal caused the handler to trigger.

138 is FPE$STACKOVERFLOW.

As I said earlier, do not call signalqq. If you let the err= handle it, it might just tell you which is the offending line.

It is possible that you are reading past end of file: do you have a end= in your read statements?

Debugging using FPS4.
[ol 1]
[li]Open your program in FPS4.[/li]
[li]Select the first executable statement, press F9. You will see a red dot on the left edge.[/li]
[li]Run the program using F5.[/li]
[li]It will stop where the red dot is with a yellow arrow.[/li]
[li]Select Debug from the menu bar. Select Exceptions... Scroll down the list - change the exception handling to Stop Always for all the float exceptions (C000008D..C0000093][/li]
[li]Click OK[/li]
[li]Press F5 to continue your program and wait until it hits the breakpoint[/li]
[li]This will tell you where the fault occured. If you use the read into a text buffer then read the numbers from the text buffer technique specified above, and examine the text buffer, it will tell you the contents of the line causing the problem.[/li]
[/ol]
It is possible that the error may not be where you think it is. You will get floating point overflow for real numbers > 1E38. If you are expecting numbers > 1E38, change to double precision.