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` or any of the values passable to `PetscSetFPTrap()`

 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:   struct PetscFPTrapLink *link;

 55:   PetscNew(&link);
 56: #if defined(PETSC_HAVE_THREADSAFETY) && defined(PETSC_HAVE_OPENMP)
 57:   #pragma omp critical
 58: #endif
 59:   {
 60:     link->trapmode = _trapmode;
 61:     link->next     = _trapstack;
 62:     _trapstack     = link;
 63:   }
 64:   if (trap != _trapmode) PetscSetFPTrap(trap);
 65:   return 0;
 66: }

 68: /*@
 69:    PetscFPTrapPop - push a floating point trapping mode, to be restored using `PetscFPTrapPop()`

 71:    Not Collective

 73:    Level: advanced

 75: .seealso: `PetscFPTrapPush()`, `PetscSetFPTrap()`, `PetscDetermineInitialFPTrap()`
 76: @*/
 77: PetscErrorCode PetscFPTrapPop(void)
 78: {
 79:   struct PetscFPTrapLink *link;

 81:   if (_trapstack->trapmode != _trapmode) PetscSetFPTrap(_trapstack->trapmode);
 82: #if defined(PETSC_HAVE_THREADSAFETY) && defined(PETSC_HAVE_OPENMP)
 83:   #pragma omp critical
 84: #endif
 85:   {
 86:     link       = _trapstack;
 87:     _trapstack = _trapstack->next;
 88:   }
 89:   PetscFree(link);
 90:   return 0;
 91: }

 93: /*--------------------------------------- ---------------------------------------------------*/
 94: #if defined(PETSC_HAVE_SUN4_STYLE_FPTRAP)
 95:   #include <floatingpoint.h>

 97: PETSC_EXTERN PetscErrorCode ieee_flags(char *, char *, char *, char **);
 98: PETSC_EXTERN PetscErrorCode ieee_handler(char *, char *, sigfpe_handler_type(int, int, struct sigcontext *, char *));

100: static struct {
101:   int   code_no;
102:   char *name;
103: } error_codes[] = {
104:   {FPE_INTDIV_TRAP,   "integer divide"               },
105:   {FPE_FLTOPERR_TRAP, "IEEE operand error"           },
106:   {FPE_FLTOVF_TRAP,   "floating point overflow"      },
107:   {FPE_FLTUND_TRAP,   "floating point underflow"     },
108:   {FPE_FLTDIV_TRAP,   "floating pointing divide"     },
109:   {FPE_FLTINEX_TRAP,  "inexact floating point result"},
110:   {0,                 "unknown error"                }
111: };
112:   #define SIGPC(scp) (scp->sc_pc)

114: /* this function gets called if a trap has occurred and been caught */
115: sigfpe_handler_type PetscDefaultFPTrap(int sig, int code, struct sigcontext *scp, char *addr)
116: {
117:   int err_ind = -1;

119:   for (int j = 0; error_codes[j].code_no; j++) {
120:     if (error_codes[j].code_no == code) err_ind = j;
121:   }

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

126:   (void)PetscError(PETSC_COMM_SELF, PETSC_ERR_FP, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
127:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
128:   return 0;
129: }

131: /*@
132:    PetscSetFPTrap - Enables traps/exceptions on common floating point errors. This option may not work on certain systems or only a
133:    subset of exceptions may be trapable.

135:    Not Collective

137:    Input Parameters:
138: .  flag - values are
139: .vb
140:     PETSC_FP_TRAP_OFF   - do not trap any exceptions
141:     PETSC_FP_TRAP_ON - all exceptions that are possible on the system except underflow
142:     PETSC_FP_TRAP_INDIV - integer divide by zero
143:     PETSC_FP_TRAP_FLTOPERR - improper argument to function, for example with real numbers, the square root of a negative number
144:     PETSC_FP_TRAP_FLTOVF - overflow
145:     PETSC_FP_TRAP_FLTUND - underflow - not trapped by default on most systems
146:     PETSC_FP_TRAP_FLTDIV - floating point divide by zero
147:     PETSC_FP_TRAP_FLTINEX - inexact floating point result
148: .ve

150:    Options Database Keys:
151: .  -fp_trap <off,on> - turn on or off trapping of floating point exceptions

153:    Level: advanced

155:    Notes:
156:    Currently only `PETSC_FP_TRAP_OFF` and `PETSC_FP_TRAP_ON` are handled. All others are treated as `PETSC_FP_TRAP_ON`.

158:    The values are bit values and may be |ed together in the function call

160:    On systems that support it this routine causes floating point
161:    overflow, divide-by-zero, and invalid-operand (e.g., a NaN), but not underflow, to
162:    cause a message to be printed and the program to exit.

164:    On many common systems, the floating
165:    point exception state is not preserved from the location where the trap
166:    occurred through to the signal handler.  In this case, the signal handler
167:    will just say that an unknown floating point exception occurred and which
168:    function it occurred in.  If you run with -fp_trap in a debugger, it will
169:    break on the line where the error occurred.  On systems that support C99
170:    floating point exception handling You can check which
171:    exception occurred using fetestexcept(FE_ALL_EXCEPT).  See fenv.h
172:    (usually at /usr/include/bits/fenv.h) for the enum values on your system.

174: .seealso: `PetscFPTrapPush()`, `PetscFPTrapPop()`, `PetscDetermineInitialFPTrap()`
175: @*/
176: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
177: {
178:   char *out;

180:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
181:   (void)ieee_flags("clear", "exception", "all", &out);
182:   if (flag == PETSC_FP_TRAP_ON) {
183:     /*
184:       To trap more fp exceptions, including underflow, change the line below to
185:       if (ieee_handler("set","all",PetscDefaultFPTrap)) {
186:     */
187:     if (ieee_handler("set", "common", PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
188:   } else if (ieee_handler("clear", "common", PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

190:   _trapmode = flag;
191:   PetscInfo(NULL, "Using PETSC_HAVE_SUN4_STYLE_FPTRAP\n");
192:   return 0;
193: }

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

198:    Not Collective

200:    Note:
201:       Currently only supported on Linux and MacOS. Checks if divide by zero is enable and if so declares that trapping is on.

203:    Level: advanced

205: .seealso: `PetscFPTrapPush()`, `PetscFPTrapPop()`, `PetscDetermineInitialFPTrap()`
206: @*/
207: PetscErrorCode PetscDetermineInitialFPTrap(void)
208: {
209:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
210:   return 0;
211: }

213: /* -------------------------------------------------------------------------------------------*/
214: #elif defined(PETSC_HAVE_SOLARIS_STYLE_FPTRAP)
215:   #include <sunmath.h>
216:   #include <floatingpoint.h>
217:   #include <siginfo.h>
218:   #include <ucontext.h>

220: static struct {
221:   int   code_no;
222:   char *name;
223: } error_codes[] = {
224:   {FPE_FLTINV, "invalid floating point operand"},
225:   {FPE_FLTRES, "inexact floating point result" },
226:   {FPE_FLTDIV, "division-by-zero"              },
227:   {FPE_FLTUND, "floating point underflow"      },
228:   {FPE_FLTOVF, "floating point overflow"       },
229:   {0,          "unknown error"                 }
230: };
231:   #define SIGPC(scp) (scp->si_addr)

233: void PetscDefaultFPTrap(int sig, siginfo_t *scp, ucontext_t *uap)
234: {
235:   int err_ind = -1, code = scp->si_code;

237:   for (int j = 0; error_codes[j].code_no; j++) {
238:     if (error_codes[j].code_no == code) err_ind = j;
239:   }

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

244:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
245:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
246: }

248: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
249: {
250:   char *out;

252:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
253:   (void)ieee_flags("clear", "exception", "all", &out);
254:   if (flag == PETSC_FP_TRAP_ON) {
255:     if (ieee_handler("set", "common", (sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floating point handler\n");
256:   } else {
257:     if (ieee_handler("clear", "common", (sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");
258:   }
259:   _trapmode = flag;
260:   PetscInfo(NULL,"Using PETSC_HAVE_SOLARIS_STYLE_FPTRAP\n";
261:   return 0;
262: }

264: PetscErrorCode PetscDetermineInitialFPTrap(void)
265: {
266:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
267:   return 0;
268: }

270: /* ------------------------------------------------------------------------------------------*/
271: #elif defined(PETSC_HAVE_IRIX_STYLE_FPTRAP)
272:   #include <sigfpe.h>
273: static struct {
274:   int   code_no;
275:   char *name;
276: } error_codes[] = {
277:   {_INVALID, "IEEE operand error"      },
278:   {_OVERFL,  "floating point overflow" },
279:   {_UNDERFL, "floating point underflow"},
280:   {_DIVZERO, "floating point divide"   },
281:   {0,        "unknown error"           }
282: };
283: void PetscDefaultFPTrap(unsigned exception[], int val[])
284: {
285:   int err_ind = -1, code = exception[0];

287:   for (int j = 0; error_codes[j].code_no; j++) {
288:     if (error_codes[j].code_no == code) err_ind = j;
289:   }
290:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred ***\n", error_codes[err_ind].name);
291:   else (*PetscErrorPrintf)("*** floating point error 0x%x occurred ***\n", code);

293:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
294:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
295: }

297: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
298: {
299:   if (flag == PETSC_FP_TRAP_ON) handle_sigfpes(_ON, , _EN_UNDERFL | _EN_OVERFL | _EN_DIVZERO | _EN_INVALID, PetscDefaultFPTrap, _ABORT_ON_ERROR, 0);
300:   else handle_sigfpes(_OFF, _EN_UNDERFL | _EN_OVERFL | _EN_DIVZERO | _EN_INVALID, 0, _ABORT_ON_ERROR, 0);
301:   _trapmode = flag;
302:   PetscInfo(NULL, "Using PETSC_HAVE_IRIX_STYLE_FPTRAP\n");
303:   return 0;
304: }

306: PetscErrorCode PetscDetermineInitialFPTrap(void)
307: {
308:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
309:   return 0;
310: }

312: /*----------------------------------------------- --------------------------------------------*/
313: #elif defined(PETSC_HAVE_RS6000_STYLE_FPTRAP)
314: /* In "fast" mode, floating point traps are imprecise and ignored.
315:    This is the reason for the fptrap(FP_TRAP_SYNC) call */
316: struct sigcontext;
317:   #include <fpxcp.h>
318:   #include <fptrap.h>
319:   #define FPE_FLTOPERR_TRAP (fptrap_t)(0x20000000)
320:   #define FPE_FLTOVF_TRAP   (fptrap_t)(0x10000000)
321:   #define FPE_FLTUND_TRAP   (fptrap_t)(0x08000000)
322:   #define FPE_FLTDIV_TRAP   (fptrap_t)(0x04000000)
323:   #define FPE_FLTINEX_TRAP  (fptrap_t)(0x02000000)

325: static struct {
326:   int   code_no;
327:   char *name;
328: } error_codes[] = {
329:   {FPE_FLTOPERR_TRAP, "IEEE operand error"           },
330:   {FPE_FLTOVF_TRAP,   "floating point overflow"      },
331:   {FPE_FLTUND_TRAP,   "floating point underflow"     },
332:   {FPE_FLTDIV_TRAP,   "floating point divide"        },
333:   {FPE_FLTINEX_TRAP,  "inexact floating point result"},
334:   < {0,                 "unknown error"                }
335: };
336:   #define SIGPC(scp)        (0) /* Info MIGHT be in scp->sc_jmpbuf.jmp_context.iar */
337: /*
338:    For some reason, scp->sc_jmpbuf does not work on the RS6000, even though
339:    it looks like it should from the include definitions.  It is probably
340:    some strange interaction with the "POSIX_SOURCE" that we require.
341: */

343: void PetscDefaultFPTrap(int sig, int code, struct sigcontext *scp)
344: {
345:   int      err_ind, j;
346:   fp_ctx_t flt_context;

348:   fp_sh_trap_info(scp, &flt_context);

350:   err_ind = -1;
351:   for (j = 0; error_codes[j].code_no; j++) {
352:     if (error_codes[j].code_no == flt_context.trap) err_ind = j;
353:   }

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

358:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
359:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
360: }

362: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
363: {
364:   if (flag == PETSC_FP_TRAP_ON) {
365:     signal(SIGFPE, (void (*)(int))PetscDefaultFPTrap);
366:     fp_trap(FP_TRAP_SYNC);
367:     /* fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW); */
368:     fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW);
369:   } else {
370:     signal(SIGFPE, SIG_DFL);
371:     fp_disable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW);
372:     fp_trap(FP_TRAP_OFF);
373:   }
374:   _trapmode = flag;
375:   PetscInfo(NULL, "Using PETSC_HAVE_RS6000_STYLE_FPTRAP\n");
376:   return 0;
377: }

379: PetscErrorCode PetscDetermineInitialFPTrap(void)
380: {
381:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
382:   return 0;
383: }

385: /* ------------------------------------------------------------*/
386: #elif defined(PETSC_HAVE_WINDOWS_COMPILERS)
387:   #include <float.h>
388: void PetscDefaultFPTrap(int sig)
389: {
390:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
391:   PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
392:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
393: }

395: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
396: {
397:   unsigned int cw;

399:   if (flag == PETSC_FP_TRAP_ON) {
400:     /* cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW | _EM_UNDERFLOW; */
401:     cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW;
403:   } else {
404:     cw = 0;
406:   }
407:   (void)_controlfp(0, cw);
408:   _trapmode = flag;
409:   PetscInfo(NULL, "Using PETSC_HAVE_WINDOWS_COMPILERS FPTRAP\n");
410:   return 0;
411: }

413: PetscErrorCode PetscDetermineInitialFPTrap(void)
414: {
415:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
416:   return 0;
417: }

419: /* ------------------------------------------------------------*/
420: #elif defined(PETSC_HAVE_FENV_H) && !defined(__cplusplus)
421:   /*
422:    C99 style floating point environment.

424:    Note that C99 merely specifies how to save, restore, and clear the floating
425:    point environment as well as defining an enumeration of exception codes.  In
426:    particular, C99 does not specify how to make floating point exceptions raise
427:    a signal.  Glibc offers this capability through FE_NOMASK_ENV (or with finer
428:    granularity, feenableexcept()), xmmintrin.h offers _MM_SET_EXCEPTION_MASK().
429: */
430:   #include <fenv.h>
431: typedef struct {
432:   int         code;
433:   const char *name;
434: } FPNode;
435: static const FPNode error_codes[] = {
436:   {FE_DIVBYZERO, "divide by zero"                                 },
437:   {FE_INEXACT,   "inexact floating point result"                  },
438:   {FE_INVALID,   "invalid floating point arguments (domain error)"},
439:   {FE_OVERFLOW,  "floating point overflow"                        },
440:   {FE_UNDERFLOW, "floating point underflow"                       },
441:   {0,            "unknown error"                                  }
442: };

444: void PetscDefaultFPTrap(int sig)
445: {
446:   const FPNode *node;
447:   int           code;
448:   PetscBool     matched = PETSC_FALSE;

450:   /* Note: While it is possible for the exception state to be preserved by the
451:    * kernel, this seems to be rare which makes the following flag testing almost
452:    * useless.  But on a system where the flags can be preserved, it would provide
453:    * more detail.
454:    */
455:   code = fetestexcept(FE_ALL_EXCEPT);
456:   for (node = &error_codes[0]; node->code; node++) {
457:     if (code & node->code) {
458:       matched = PETSC_TRUE;
459:       (*PetscErrorPrintf)("*** floating point error \"%s\" occurred ***\n", node->name);
460:       code &= ~node->code; /* Unset this flag since it has been processed */
461:     }
462:   }
463:   if (!matched || code) { /* If any remaining flags are set, or we didn't process any flags */
464:     (*PetscErrorPrintf)("*** unknown floating point error occurred ***\n");
465:     (*PetscErrorPrintf)("The specific exception can be determined by running in a debugger.  When the\n");
466:     (*PetscErrorPrintf)("debugger traps the signal, the exception can be found with fetestexcept(0x%x)\n", FE_ALL_EXCEPT);
467:     (*PetscErrorPrintf)("where the result is a bitwise OR of the following flags:\n");
468:     (*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);
469:   }

471:   (*PetscErrorPrintf)("Try option -start_in_debugger\n");
472:   #if PetscDefined(USE_DEBUG)
473:   (*PetscErrorPrintf)("likely location of problem given in stack below\n");
474:   (*PetscErrorPrintf)("---------------------  Stack Frames ------------------------------------\n");
475:   PetscStackView(PETSC_STDOUT);
476:   #else
477:   (*PetscErrorPrintf)("configure using --with-debugging=yes, recompile, link, and run \n");
478:   (*PetscErrorPrintf)("with -start_in_debugger to get more information on the crash.\n");
479:   #endif
480:   PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_INITIAL, "trapped floating point error");
481:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
482: }

484: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
485: {
486:   if (flag == PETSC_FP_TRAP_ON) {
487:     /* Clear any flags that are currently set so that activating trapping will not immediately call the signal handler. */
489:   #if defined(FE_NOMASK_ENV)
490:     /* Could use fesetenv(FE_NOMASK_ENV), but that causes spurious exceptions (like gettimeofday() -> PetscLogDouble). */
492:     /* Doesn't work on AArch64 targets. There's a known hardware limitation. Need to detect hardware at configure time? */
494:     PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP with FE_NOMASK_ENV\n");
495:   #elif defined PETSC_HAVE_XMMINTRIN_H
496:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO);
497:     /* _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW); */
498:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
499:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
500:     PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP with PETSC_HAVE_XMMINTRIN_H\n");
501:   #else
502:     /* C99 does not provide a way to modify the environment so there is no portable way to activate trapping. */
503:     PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP\n");
504:   #endif
506:   } else {
508:     /* can use _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_UNDERFLOW); if PETSC_HAVE_XMMINTRIN_H exists */
510:   }
511:   _trapmode = flag;
512:   return 0;
513: }

515: PetscErrorCode PetscDetermineInitialFPTrap(void)
516: {
517:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
518:   unsigned int flags;
519:   #endif

521:   #if defined(FE_NOMASK_ENV)
522:   flags = fegetexcept();
523:   if (flags & FE_DIVBYZERO) {
524:   #elif defined PETSC_HAVE_XMMINTRIN_H
525:   flags = _MM_GET_EXCEPTION_MASK();
526:   if (!(flags & _MM_MASK_DIV_ZERO)) {
527:   #else
528:   PetscInfo(NULL, "Floating point trapping unknown, assuming off\n");
529:   return 0;
530:   #endif
531:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
532:     _trapmode = PETSC_FP_TRAP_ON;
533:     PetscInfo(NULL, "Floating point trapping is on by default %d\n", flags);
534:   } else {
535:     _trapmode = PETSC_FP_TRAP_OFF;
536:     PetscInfo(NULL, "Floating point trapping is off by default %d\n", flags);
537:   }
538:   return 0;
539:   #endif
540: }

542: /* ------------------------------------------------------------*/
543: #elif defined(PETSC_HAVE_IEEEFP_H)
544:   #include <ieeefp.h>
545: void PetscDefaultFPTrap(int sig)
546: {
547:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
548:   PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
549:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
550: }

552: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
553: {
554:   if (flag == PETSC_FP_TRAP_ON) {
555:   #if defined(PETSC_HAVE_FPPRESETSTICKY)
556:     fpresetsticky(fpgetsticky());
557:   #elif defined(PETSC_HAVE_FPSETSTICKY)
558:     fpsetsticky(fpgetsticky());
559:   #endif
560:     fpsetmask(FP_X_INV | FP_X_DZ | FP_X_OFL | FP_X_OFL);
562:   } else {
563:   #if defined(PETSC_HAVE_FPPRESETSTICKY)
564:     fpresetsticky(fpgetsticky());
565:   #elif defined(PETSC_HAVE_FPSETSTICKY)
566:     fpsetsticky(fpgetsticky());
567:   #endif
568:     fpsetmask(0);
570:   }
571:   _trapmode = flag;
572:   PetscInfo(NULL, "Using PETSC_HAVE_IEEEFP_H FPTRAP\n");
573:   return 0;
574: }

576: PetscErrorCode PetscDetermineInitialFPTrap(void)
577: {
578:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
579:   return 0;
580: }

582: /* -------------------------Default -----------------------------------*/
583: #else

585: void PetscDefaultFPTrap(int sig)
586: {
587:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
588:   PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
589:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
590: }

592: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
593: {
594:   if (flag == PETSC_FP_TRAP_ON) {
595:     if (SIG_ERR == signal(SIGFPE, PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
596:   } else if (SIG_ERR == signal(SIGFPE, SIG_DFL)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

598:   _trapmode = flag;
599:   PetscInfo(NULL, "Using default FPTRAP\n");
600:   return 0;
601: }

603: PetscErrorCode PetscDetermineInitialFPTrap(void)
604: {
605:   PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n");
606:   return 0;
607: }
608: #endif