Xref: utzoo comp.unix.ultrix:7443 comp.lang.fortran:5563 Path: utzoo!utgpu!news-server.csri.toronto.edu!bonnie.concordia.ca!uunet!zaphod.mps.ohio-state.edu!hobbes.physics.uiowa.edu!news.iastate.edu!IASTATE.EDU!keinert From: keinert@IASTATE.EDU (Keinert Fritz) Newsgroups: comp.unix.ultrix,comp.lang.fortran Subject: Re: How to detect NaN's Message-ID: <1991May31.130807@IASTATE.EDU> Date: 31 May 91 18:08:07 GMT References: <1991May30.204332.16506@litwin.com> Sender: news@news.iastate.edu (USENET News System) Reply-To: keinert@IASTATE.EDU (Keinert Fritz) Organization: Iowa State University Lines: 293 In article <1991May30.204332.16506@litwin.com>, vlr@litwin.com (Vic Rice) writes: > I have encountered a problem on a DECStation 5000 using the Mips Fortran > compiler. Some operation I am performing is generating a NaN which then > proceeds to propogate to other dependent variables. Since this is in > a loop, I am not exactly sure where it starts. > > How can I test a variable for the presence of a NaN or Infinity ??? > -- > Dr. Victor L. Rice > Litwin Process Automation One of my biggest gripes with both the Sun and Mips Fortran compilers is that there is no easy way to turn on floating point exception handling. This results in exactly the problem Vic describes: at the end of the program, most of your output is infinity or NaN, and you have no clue where it started. Here is a simple floating point exception handler that I wrote (Mips Fortran version). Compile it with C and link with your Fortran (or C) program, like this: % cc -c fpx_mips.c % f77 prog.f fpx_mips.o -o prog or something like that. To use it, just insert a statement call standard_fpx() or call setcsr(whatever) as the first executable statement in your Fortran program. I welcome suggestions for improvement, but I am not prepared to answer a lot of questions or fix bugs in this code. Good Luck. A last remark for the knowledgeable: the reason I am trapping SIGILL is a known bug in Ultrix 4.1. Floating point errors cause SIGILL instead of SIGFPE. -- Fritz Keinert phone: (515) 294-5128 Department of Mathematics fax: (515) 294-5454 Iowa State University e-mail: keinert@iastate.edu Ames, IA 50011 -------------------------------------------------------------------------- /********************************************************************** * * * Floating point exception handler for MIPS machines * * Fritz Keinert (keinert@iastate.edu) * * April 5, 1991 * * * * building on routines by Howard Jespersen and Srini Uppugunduri * * * * As they say in NETLIB: Everything free comes with no guarantee. * * * * This file contains the following routines: * * * * int getcsr() - returns the contents of the control/status * * register (CSR) of the MIPS floating point * * unit (FPU). The register contents are * * described in /usr/include/mips/fpu.h * * * * int getcsr_() - same routine, callable from Fortran * * * * int setcsr(csr) - sets CSR to new value, returns old value; * * also installs fpx_handler * * * * int setcsr_(csr) - same routine, callable from Fortran * * * * void fpx_handler(...) - floating point exception handler. Prints * * an error message and dies. * * * * int standard_fpx() - same as setcsr(3584), to set up the * * most frequently useful behavior. * * * * * * int standard_fpx_() - same routine, callable from Fortran * * * * * ********************************************************************** * * * The floating point exception handling works like this: * * * * First, you have to tell the floating point unit to enable * * interrupts, by setting the CSR accordingly. (I assume this * * stands for Control and Status Register). The correct value * * for CSR is found by adding * * * * one of the following, to define the rounding behavior: * * * * 0 - round to nearest * * 1 - zero * * 2 - + infinity * * 3 - - infinity * * * * plus as many as needed of the following, to turn on * * interrupts: * * * * 128 - cause interrupt on inexact result * * 256 - underflow * * 512 - overflow * * 1024 - division by zero * * 2048 - invalid operation * * * * The standard values are: * * * * 1. round to nearest, interrupt on overflow, division by zero * * and invalid operation. CSR = 3584. * * * * 2. round to nearest, interrupt on underflow, overflow, * * division by zero and invalid operation. CSR = 3840. * * * * Second, you have to install an interrupt handler to give you * * a message when an interrupt actually happens. Otherwise, the * * operating system will just ignore it, or kill your program with * * a highly informative message like "illegal operation". * * * * My interrupt handler just prints a message and dies. For more * * sophisticated behavior, like "print a message on underflow, set * * the result to zero, and continue", or a traceback of where this * * happened, you are on your own. * * * * Use the dbx debugger to find out where in your source code the * * problem is. Here is an example: program test dies with the * * message * * * * caught signal 4, code 0 * * floating point error: overflow * * at or near location 0x400448 * * * * Call up dbx, execute the first line of your program, set the * * program counter to the address reported by the interrupt routine, * * and ask dbx where that is: * * * * % dbx test * * dbx version 2.0 * * Type 'help' for help. * * reading symbolic information ... * * [using test.test] * * test: 6 print*,'round to nearest = 0' * * (dbx) stop at 6 * * [2] stop at "test.f":6 * * (dbx) run * * [2] stopped at [test.test:6 ,0x40020c] * * print*,'round to nearest = 0' * * (dbx) assign $pc=0x400448 * * 4195400 * * (dbx) where * * > 0 test.test() ["test.f":24, 0x400448] * * 1 main.main(0x7fffb9c4, 0x7fffb9cc, 0x0, 0x0, 0x0) * * ["main.c":35, 0x400f2c] * * (dbx) quit * * * * and we know: it was around line 24 in file test.f. * * * * This method may not work if the error is in a system subroutine. * * For example, if you try sqrt(-1.), this will tell you that the * * problem is inside a square root, but there may be many square * * roots in your program. * * * * If all else fails, just run your program inside dbx until it * * dies, and then type "where". This will get you an exact * * traceback of where you are (sqrt function, called from line * * ?? of subroutine ??, which was called from line ?? in main * * program). * * * * You may have to use the "-g" flag for compiling, to get the * * debugger to work properly. * * * **********************************************************************/ #include #include #include #include void fpx_handler(int sig, int code, struct sigcontext *scp) { union fpc_csr csr; csr.fc_word = scp->sc_fpc_csr; fprintf(stderr,"caught signal %d, code %d\n", sig, code); if (csr.fc_struct.ex_invalid) { fprintf(stderr,"floating point error: invalid operation\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (csr.fc_struct.ex_divide0) { fprintf(stderr,"floating point error: division by zero\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (csr.fc_struct.ex_underflow) { fprintf(stderr,"floating point error: underflow\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (csr.fc_struct.ex_overflow) { fprintf(stderr,"floating point error: overflow\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (csr.fc_struct.ex_inexact) { fprintf(stderr,"floating point error: inexact operation\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (csr.fc_struct.ex_unimplemented) { fprintf(stderr,"floating point error: unimplemented operation\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } if (code == BRK_DIVZERO) { fprintf(stderr,"integer error: division by zero\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); exit(1); } exit(1); fprintf(stderr,"floating point error: unknown exception code\n"); fprintf(stderr,"at or near location 0x%x\n",scp->sc_pc); fprintf(stderr,"fp status word 0x%x\n",csr); exit(1); } int getcsr() /* * return contents of fpu control and status register */ { return get_fpc_csr(); } int getcsr_() /* * Fortran wrapper for getcsr() */ { return getcsr(); } int setcsr(unsigned int csr) /* * set up floating point interrupt handler, turn on interrupts */ { signal(SIGTRAP, fpx_handler); signal(SIGFPE, fpx_handler); signal(SIGILL, fpx_handler); return set_fpc_csr(csr); } int setcsr_(unsigned int *csr) /* * Fortran wrapper for setcsr() */ { return setcsr(*csr); } int standard_fpx() /* * set up standard error handling: die on overflow, * division by zero, illegal operation, but continue * on underflow. Round to nearest. */ { return setcsr((unsigned int) 3584); } int standard_fpx_() /* * Fortran wrapper for standard_fpx() */ { return standard_fpx(); }