Actual source code: fp.c


  2: /*
  3:    IEEE error handler for all machines. Since each OS has
  4:    enough slight differences we have completely separate codes for each one.
  5: */

  7: /*
  8:   This feature test macro provides FE_NOMASK_ENV on GNU.  It must be defined
  9:   at the top of the file because other headers may pull in fenv.h even when
 10:   not strictly necessary.  Strictly speaking, we could include ONLY petscconf.h,
 11:   check PETSC_HAVE_FENV_H, and only define _GNU_SOURCE in that case, but such
 12:   shenanigans ought to be unnecessary.
 13: */
 14: #if !defined(_GNU_SOURCE)
 15: #define _GNU_SOURCE
 16: #endif

 18: #include <petsc/private/petscimpl.h>
 19: #include <signal.h>

 21: struct PetscFPTrapLink {
 22:   PetscFPTrap            trapmode;
 23:   struct PetscFPTrapLink *next;
 24: };
 25: static PetscFPTrap            _trapmode = PETSC_FP_TRAP_OFF; /* Current trapping mode; see PetscDetermineInitialFPTrap() */
 26: static struct PetscFPTrapLink *_trapstack;                   /* Any pushed states of _trapmode */

 28: /*@
 29:    PetscFPTrapPush - push a floating point trapping mode, restored using PetscFPTrapPop()

 31:    Not Collective

 33:    Input Parameter:
 34: .    trap - PETSC_FP_TRAP_ON or PETSC_FP_TRAP_OFF

 36:    Level: advanced

 38:    Notes:
 39:      This only changes the trapping if the new mode is different than the current mode.

 41:      This routine is called to turn off trapping for certain LAPACK routines that assume that dividing
 42:      by zero is acceptable. In particular the routine ieeeck().

 44:      Most systems by default have all trapping turned off, but certain Fortran compilers have
 45:      link flags that turn on trapping before the program begins.
 46: $       gfortran -ffpe-trap=invalid,zero,overflow,underflow,denormal
 47: $       ifort -fpe0

 49: .seealso: PetscFPTrapPop(), PetscSetFPTrap(), PetscDetermineInitialFPTrap()
 50: @*/
 51: PetscErrorCode PetscFPTrapPush(PetscFPTrap trap)
 52: {
 53:   PetscErrorCode         ierr;
 54:   struct PetscFPTrapLink *link;

 57:   PetscNew(&link);
 58:   link->trapmode = _trapmode;
 59:   link->next     = _trapstack;
 60:   _trapstack     = link;
 61:   if (trap != _trapmode) {PetscSetFPTrap(trap);}
 62:   return(0);
 63: }

 65: /*@
 66:    PetscFPTrapPop - push a floating point trapping mode, to be restored using PetscFPTrapPop()

 68:    Not Collective

 70:    Level: advanced

 72: .seealso: PetscFPTrapPush(), PetscSetFPTrap(), PetscDetermineInitialFPTrap()
 73: @*/
 74: PetscErrorCode PetscFPTrapPop(void)
 75: {
 76:   PetscErrorCode         ierr;
 77:   struct PetscFPTrapLink *link;

 80:   if (_trapstack->trapmode != _trapmode) {PetscSetFPTrap(_trapstack->trapmode);}
 81:   link       = _trapstack;
 82:   _trapstack = _trapstack->next;
 83:   PetscFree(link);
 84:   return(0);
 85: }

 87: /*--------------------------------------- ---------------------------------------------------*/
 88: #if defined(PETSC_HAVE_SUN4_STYLE_FPTRAP)
 89: #include <floatingpoint.h>

 91: PETSC_EXTERN PetscErrorCode ieee_flags(char*,char*,char*,char**);
 92: PETSC_EXTERN PetscErrorCode ieee_handler(char*,char*,sigfpe_handler_type(int,int,struct sigcontext*,char*));

 94: static struct { int code_no; char *name; } error_codes[] = {
 95:   { FPE_INTDIV_TRAP    ,"integer divide" },
 96:   { FPE_FLTOPERR_TRAP  ,"IEEE operand error" },
 97:   { FPE_FLTOVF_TRAP    ,"floating point overflow" },
 98:   { FPE_FLTUND_TRAP    ,"floating point underflow" },
 99:   { FPE_FLTDIV_TRAP    ,"floating pointing divide" },
100:   { FPE_FLTINEX_TRAP   ,"inexact floating point result" },
101:   { 0                  ,"unknown error" }
102: };
103: #define SIGPC(scp) (scp->sc_pc)

105: sigfpe_handler_type PetscDefaultFPTrap(int sig,int code,struct sigcontext *scp,char *addr)
106: {
108:   int            err_ind = -1,j;

111:   for (j = 0; error_codes[j].code_no; j++) {
112:     if (error_codes[j].code_no == code) err_ind = j;
113:   }

115:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred at pc=%X ***\n",error_codes[err_ind].name,SIGPC(scp));
116:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred at pc=%X ***\n",code,SIGPC(scp));

118:   PetscError(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
119:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
120:   return(0);
121: }

123: /*@
124:    PetscSetFPTrap - Enables traps/exceptions on common floating point errors.
125:                     This option may not work on certain systems.

127:    Not Collective

129:    Input Parameters:
130: .  flag - PETSC_FP_TRAP_ON, PETSC_FP_TRAP_OFF.

132:    Options Database Keys:
133: .  -fp_trap - Activates floating point trapping

135:    Level: advanced

137:    Description:
138:    On systems that support it, when called with PETSC_FP_TRAP_ON this routine causes floating point
139:    underflow, overflow, divide-by-zero, and invalid-operand (e.g., a NaN) to
140:    cause a message to be printed and the program to exit.

142:    Note:
143:    On many common systems, the floating
144:    point exception state is not preserved from the location where the trap
145:    occurred through to the signal handler.  In this case, the signal handler
146:    will just say that an unknown floating point exception occurred and which
147:    function it occurred in.  If you run with -fp_trap in a debugger, it will
148:    break on the line where the error occurred.  On systems that support C99
149:    floating point exception handling You can check which
150:    exception occurred using fetestexcept(FE_ALL_EXCEPT).  See fenv.h
151:    (usually at /usr/include/bits/fenv.h) for the enum values on your system.

153:    Caution:
154:    On certain machines, in particular the IBM PowerPC, floating point
155:    trapping may be VERY slow!

157: .seealso: PetscFPTrapPush(), PetscFPTrapPop(), PetscDetermineInitialFPTrap()
158: @*/
159: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
160: {
161:   char *out;

164:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
165:   (void) ieee_flags("clear","exception","all",&out);
166:   if (flag == PETSC_FP_TRAP_ON) {
167:     /*
168:       To trap more fp exceptions, including underflow, change the line below to
169:       if (ieee_handler("set","all",PetscDefaultFPTrap)) {
170:     */
171:     if (ieee_handler("set","common",PetscDefaultFPTrap))        (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
172:   } else if (ieee_handler("clear","common",PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

174:   _trapmode = flag;
175:   return(0);
176: }

178: /*@
179:    PetscDetermineInitialFPTrap - Attempts to determine the floating point trapping that exists when PetscInitialize() is called

181:    Not Collective

183:    Notes:
184:       Currently only supported on Linux and MacOS. Checks if divide by zero is enable and if so declares that trapping is on.

186:    Level: advanced

188: .seealso: PetscFPTrapPush(), PetscFPTrapPop(), PetscDetermineInitialFPTrap()
189: @*/
190: PetscErrorCode  PetscDetermineInitialFPTrap(void)
191: {

195:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
196:   return(0);
197: }

199: /* -------------------------------------------------------------------------------------------*/
200: #elif defined(PETSC_HAVE_SOLARIS_STYLE_FPTRAP)
201: #include <sunmath.h>
202: #include <floatingpoint.h>
203: #include <siginfo.h>
204: #include <ucontext.h>

206: static struct { int code_no; char *name; } error_codes[] = {
207:   { FPE_FLTINV,"invalid floating point operand"},
208:   { FPE_FLTRES,"inexact floating point result"},
209:   { FPE_FLTDIV,"division-by-zero"},
210:   { FPE_FLTUND,"floating point underflow"},
211:   { FPE_FLTOVF,"floating point overflow"},
212:   { 0,         "unknown error"}
213: };
214: #define SIGPC(scp) (scp->si_addr)

216: void PetscDefaultFPTrap(int sig,siginfo_t *scp,ucontext_t *uap)
217: {
218:   int            err_ind,j,code = scp->si_code;

222:   err_ind = -1;
223:   for (j = 0; error_codes[j].code_no; j++) {
224:     if (error_codes[j].code_no == code) err_ind = j;
225:   }

227:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred at pc=%X ***\n",error_codes[err_ind].name,SIGPC(scp));
228:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred at pc=%X ***\n",code,SIGPC(scp));

230:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
231:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
232: }

234: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
235: {
236:   char *out;

239:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
240:   (void) ieee_flags("clear","exception","all",&out);
241:   if (flag == PETSC_FP_TRAP_ON) {
242:     if (ieee_handler("set","common",(sigfpe_handler_type)PetscDefaultFPTrap))        (*PetscErrorPrintf)("Can't set floating point handler\n");
243:   } else {
244:     if (ieee_handler("clear","common",(sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");
245:   }
246:   _trapmode = flag;
247:   return(0);
248: }

250: PetscErrorCode  PetscDetermineInitialFPTrap(void)
251: {

255:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
256:   return(0);
257: }

259: /* ------------------------------------------------------------------------------------------*/
260: #elif defined(PETSC_HAVE_IRIX_STYLE_FPTRAP)
261: #include <sigfpe.h>
262: static struct { int code_no; char *name; } error_codes[] = {
263:   { _INVALID   ,"IEEE operand error" },
264:   { _OVERFL    ,"floating point overflow" },
265:   { _UNDERFL   ,"floating point underflow" },
266:   { _DIVZERO   ,"floating point divide" },
267:   { 0          ,"unknown error" }
268: } ;
269: void PetscDefaultFPTrap(unsigned exception[],int val[])
270: {
271:   int err_ind,j,code;

274:   code    = exception[0];
275:   err_ind = -1;
276:   for (j = 0; error_codes[j].code_no; j++) {
277:     if (error_codes[j].code_no == code) err_ind = j;
278:   }
279:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred ***\n",error_codes[err_ind].name);
280:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred ***\n",code);

282:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
283:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
284: }

286: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
287: {
289:   if (flag == PETSC_FP_TRAP_ON) handle_sigfpes(_ON,,_EN_UNDERFL|_EN_OVERFL|_EN_DIVZERO|_EN_INVALID,PetscDefaultFPTrap,_ABORT_ON_ERROR,0);
290:   else                          handle_sigfpes(_OFF,_EN_UNDERFL|_EN_OVERFL|_EN_DIVZERO|_EN_INVALID,0,_ABORT_ON_ERROR,0);
291:   _trapmode = flag;
292:   return(0);
293: }

295: PetscErrorCode  PetscDetermineInitialFPTrap(void)
296: {

300:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
301:   return(0);
302: }

304: /* -------------------------------------------------------------------------------------------*/
305: #elif defined(PETSC_HAVE_SOLARIS_STYLE_FPTRAP)
306: #include <sunmath.h>
307: #include <floatingpoint.h>
308: #include <siginfo.h>
309: #include <ucontext.h>

311: static struct { int code_no; char *name; } error_codes[] = {
312:   { FPE_FLTINV,"invalid floating point operand"},
313:   { FPE_FLTRES,"inexact floating point result"},
314:   { FPE_FLTDIV,"division-by-zero"},
315:   { FPE_FLTUND,"floating point underflow"},
316:   { FPE_FLTOVF,"floating point overflow"},
317:   { 0,         "unknown error"}
318: };
319: #define SIGPC(scp) (scp->si_addr)

321: void PetscDefaultFPTrap(int sig,siginfo_t *scp,ucontext_t *uap)
322: {
323:   int            err_ind,j,code = scp->si_code;

327:   err_ind = -1;
328:   for (j = 0; error_codes[j].code_no; j++) {
329:     if (error_codes[j].code_no == code) err_ind = j;
330:   }

332:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred at pc=%X ***\n",error_codes[err_ind].name,SIGPC(scp));
333:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred at pc=%X ***\n",code,SIGPC(scp));

335:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
336:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
337: }

339: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
340: {
341:   char *out;

344:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
345:   (void) ieee_flags("clear","exception","all",&out);
346:   if (flag == PETSC_FP_TRAP_ON) {
347:     if (ieee_handler("set","common",(sigfpe_handler_type)PetscDefaultFPTrap))        (*PetscErrorPrintf)("Can't set floating point handler\n");
348:   } else {
349:     if (ieee_handler("clear","common",(sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");
350:   }
351:   _trapmode = flag;
352:   return(0);
353: }

355: PetscErrorCode  PetscDetermineInitialFPTrap(void)
356: {

360:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
361:   return(0);
362: }

364: /*----------------------------------------------- --------------------------------------------*/
365: #elif defined(PETSC_HAVE_RS6000_STYLE_FPTRAP)
366: /* In "fast" mode, floating point traps are imprecise and ignored.
367:    This is the reason for the fptrap(FP_TRAP_SYNC) call */
368: struct sigcontext;
369: #include <fpxcp.h>
370: #include <fptrap.h>
371: #define FPE_FLTOPERR_TRAP (fptrap_t)(0x20000000)
372: #define FPE_FLTOVF_TRAP   (fptrap_t)(0x10000000)
373: #define FPE_FLTUND_TRAP   (fptrap_t)(0x08000000)
374: #define FPE_FLTDIV_TRAP   (fptrap_t)(0x04000000)
375: #define FPE_FLTINEX_TRAP  (fptrap_t)(0x02000000)

377: static struct { int code_no; char *name; } error_codes[] = {
378:   {FPE_FLTOPERR_TRAP   ,"IEEE operand error" },
379:   { FPE_FLTOVF_TRAP    ,"floating point overflow" },
380:   { FPE_FLTUND_TRAP    ,"floating point underflow" },
381:   { FPE_FLTDIV_TRAP    ,"floating point divide" },
382:   { FPE_FLTINEX_TRAP   ,"inexact floating point result" },
383:   { 0                  ,"unknown error" }
384: } ;
385: #define SIGPC(scp) (0) /* Info MIGHT be in scp->sc_jmpbuf.jmp_context.iar */
386: /*
387:    For some reason, scp->sc_jmpbuf does not work on the RS6000, even though
388:    it looks like it should from the include definitions.  It is probably
389:    some strange interaction with the "POSIX_SOURCE" that we require.
390: */

392: void PetscDefaultFPTrap(int sig,int code,struct sigcontext *scp)
393: {
395:   int            err_ind,j;
396:   fp_ctx_t       flt_context;

399:   fp_sh_trap_info(scp,&flt_context);

401:   err_ind = -1;
402:   for (j = 0; error_codes[j].code_no; j++) {
403:     if (error_codes[j].code_no == flt_context.trap) err_ind = j;
404:   }

406:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred ***\n",error_codes[err_ind].name);
407:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred ***\n",flt_context.trap);

409:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
410:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
411: }

413: PetscErrorCode PetscSetFPTrap(PetscFPTrap on)
414: {
416:   if (on == PETSC_FP_TRAP_ON) {
417:     signal(SIGFPE,(void (*)(int))PetscDefaultFPTrap);
418:     fp_trap(FP_TRAP_SYNC);
419:     fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW);
420:     /* fp_enable(mask) for individual traps.  Values are:
421:        TRP_INVALID
422:        TRP_DIV_BY_ZERO
423:        TRP_OVERFLOW
424:        TRP_UNDERFLOW
425:        TRP_INEXACT
426:        Can OR then together.
427:        fp_enable_all(); for all traps.
428:     */
429:   } else {
430:     signal(SIGFPE,SIG_DFL);
431:     fp_disable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW);
432:     fp_trap(FP_TRAP_OFF);
433:   }
434:   _trapmode = on;
435:   return(0);
436: }

438: PetscErrorCode  PetscDetermineInitialFPTrap(void)
439: {

443:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
444:   return(0);
445: }

447: /* ------------------------------------------------------------*/
448: #elif defined(PETSC_HAVE_WINDOWS_COMPILERS)
449: #include <float.h>
450: void PetscDefaultFPTrap(int sig)
451: {
453:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
454:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
455:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
456: }

458: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
459: {
460:   unsigned int cw;

463:   if (on == PETSC_FP_TRAP_ON) {
464:     cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW | _EM_UNDERFLOW;
465:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't set floating point handler\n");
466:   } else {
467:     cw = 0;
468:     if (SIG_ERR == signal(SIGFPE,SIG_DFL)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't clear floating point handler\n");
469:   }
470:   (void)_controlfp(0, cw);
471:   _trapmode = on;
472:   return(0);
473: }

475: PetscErrorCode  PetscDetermineInitialFPTrap(void)
476: {

480:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
481:   return(0);
482: }

484: /* ------------------------------------------------------------*/
485: #elif defined(PETSC_HAVE_FENV_H) && !defined(__cplusplus)
486: /*
487:    C99 style floating point environment.

489:    Note that C99 merely specifies how to save, restore, and clear the floating
490:    point environment as well as defining an enumeration of exception codes.  In
491:    particular, C99 does not specify how to make floating point exceptions raise
492:    a signal.  Glibc offers this capability through FE_NOMASK_ENV (or with finer
493:    granularity, feenableexcept()), xmmintrin.h offers _MM_SET_EXCEPTION_MASK().
494: */
495: #include <fenv.h>
496: typedef struct {int code; const char *name;} FPNode;
497: static const FPNode error_codes[] = {
498:   {FE_DIVBYZERO,"divide by zero"},
499:   {FE_INEXACT,  "inexact floating point result"},
500:   {FE_INVALID,  "invalid floating point arguments (domain error)"},
501:   {FE_OVERFLOW, "floating point overflow"},
502:   {FE_UNDERFLOW,"floating point underflow"},
503:   {0           ,"unknown error"}
504: };

506: void PetscDefaultFPTrap(int sig)
507: {
508:   const FPNode *node;
509:   int          code;
510:   PetscBool    matched = PETSC_FALSE;

513:   /* Note: While it is possible for the exception state to be preserved by the
514:    * kernel, this seems to be rare which makes the following flag testing almost
515:    * useless.  But on a system where the flags can be preserved, it would provide
516:    * more detail.
517:    */
518:   code = fetestexcept(FE_ALL_EXCEPT);
519:   for (node=&error_codes[0]; node->code; node++) {
520:     if (code & node->code) {
521:       matched = PETSC_TRUE;
522:       (*PetscErrorPrintf)("*** floating point error \"%s\" occurred ***\n",node->name);
523:       code &= ~node->code; /* Unset this flag since it has been processed */
524:     }
525:   }
526:   if (!matched || code) { /* If any remaining flags are set, or we didn't process any flags */
527:     (*PetscErrorPrintf)("*** unknown floating point error occurred ***\n");
528:     (*PetscErrorPrintf)("The specific exception can be determined by running in a debugger.  When the\n");
529:     (*PetscErrorPrintf)("debugger traps the signal, the exception can be found with fetestexcept(0x%x)\n",FE_ALL_EXCEPT);
530:     (*PetscErrorPrintf)("where the result is a bitwise OR of the following flags:\n");
531:     (*PetscErrorPrintf)("FE_INVALID=0x%x FE_DIVBYZERO=0x%x FE_OVERFLOW=0x%x FE_UNDERFLOW=0x%x FE_INEXACT=0x%x\n",FE_INVALID,FE_DIVBYZERO,FE_OVERFLOW,FE_UNDERFLOW,FE_INEXACT);
532:   }

534:   (*PetscErrorPrintf)("Try option -start_in_debugger\n");
535: #if PetscDefined(USE_DEBUG)
536:   (*PetscErrorPrintf)("likely location of problem given in stack below\n");
537:   (*PetscErrorPrintf)("---------------------  Stack Frames ------------------------------------\n");
538:   PetscStackView(PETSC_STDOUT);
539: #else
540:   (*PetscErrorPrintf)("configure using --with-debugging=yes, recompile, link, and run \n");
541:   (*PetscErrorPrintf)("with -start_in_debugger to get more information on the crash.\n");
542: #endif
543:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_INITIAL,"trapped floating point error");
544:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
545: }

547: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
548: {
550:   if (on == PETSC_FP_TRAP_ON) {
551:     /* Clear any flags that are currently set so that activating trapping will not immediately call the signal handler. */
552:     if (feclearexcept(FE_ALL_EXCEPT)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot clear floating point exception flags\n");
553: #if defined(FE_NOMASK_ENV)
554:     /* Could use fesetenv(FE_NOMASK_ENV), but that causes spurious exceptions (like gettimeofday() -> PetscLogDouble). */
555:     if (feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot activate floating point exceptions\n");
556: #elif defined PETSC_HAVE_XMMINTRIN_H
557:    _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO);
558:    _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW);
559:    _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
560:    _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
561: #else
562:     /* C99 does not provide a way to modify the environment so there is no portable way to activate trapping. */
563: #endif
564:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't set floating point handler\n");
565:   } else {
566:     if (fesetenv(FE_DFL_ENV)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot disable floating point exceptions");
567:     /* can use _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_UNDERFLOW); if PETSC_HAVE_XMMINTRIN_H exists */
568:     if (SIG_ERR == signal(SIGFPE,SIG_DFL)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't clear floating point handler\n");
569:   }
570:   _trapmode = on;
571:   return(0);
572: }

574: PetscErrorCode  PetscDetermineInitialFPTrap(void)
575: {
576: #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
577:   unsigned int   flags;
578: #endif

582: #if defined(FE_NOMASK_ENV)
583:   flags = fegetexcept();
584:   if (flags & FE_DIVBYZERO) {
585: #elif defined PETSC_HAVE_XMMINTRIN_H
586:   flags = _MM_GET_EXCEPTION_MASK();
587:   if (!(flags & _MM_MASK_DIV_ZERO)) {
588: #else
589:   PetscInfo(NULL,"Floating point trapping unknown, assuming off\n");
590:   return(0);
591: #endif
592: #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
593:     _trapmode = PETSC_FP_TRAP_ON;
594:     PetscInfo1(NULL,"Floating point trapping is on by default %d\n",flags);
595:   } else {
596:     _trapmode = PETSC_FP_TRAP_OFF;
597:     PetscInfo1(NULL,"Floating point trapping is off by default %d\n",flags);
598:   }
599:   return(0);
600: #endif
601: }

603: /* ------------------------------------------------------------*/
604: #elif defined(PETSC_HAVE_IEEEFP_H)
605: #include <ieeefp.h>
606: void PetscDefaultFPTrap(int sig)
607: {
609:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
610:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
611:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
612: }

614: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
615: {
617:   if (on == PETSC_FP_TRAP_ON) {
618: #if defined(PETSC_HAVE_FPPRESETSTICKY)
619:     fpresetsticky(fpgetsticky());
620: #elif defined(PETSC_HAVE_FPSETSTICKY)
621:     fpsetsticky(fpgetsticky());
622: #endif
623:     fpsetmask(FP_X_INV | FP_X_DZ | FP_X_OFL |  FP_X_OFL);
624:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't set floating point handler\n");
625:   } else {
626: #if defined(PETSC_HAVE_FPPRESETSTICKY)
627:     fpresetsticky(fpgetsticky());
628: #elif defined(PETSC_HAVE_FPSETSTICKY)
629:     fpsetsticky(fpgetsticky());
630: #endif
631:     fpsetmask(0);
632:     if (SIG_ERR == signal(SIGFPE,SIG_DFL)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't clear floating point handler\n");
633:   }
634:   _trapmode = on;
635:   return(0);
636: }

638: PetscErrorCode  PetscDetermineInitialFPTrap(void)
639: {

643:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
644:   return(0);
645: }

647: /* -------------------------Default -----------------------------------*/
648: #else

650: void PetscDefaultFPTrap(int sig)
651: {
653:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
654:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
655:   PETSCABORT(MPI_COMM_WORLD,PETSC_ERR_FP);
656: }

658: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
659: {
661:   if (on == PETSC_FP_TRAP_ON) {
662:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
663:   } else if (SIG_ERR == signal(SIGFPE,SIG_DFL))       (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

665:   _trapmode = on;
666:   return(0);
667: }

669: PetscErrorCode  PetscDetermineInitialFPTrap(void)
670: {

674:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
675:   return(0);
676: }
677: #endif