Actual source code: fp.c

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

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

 17: #include <petsc/private/petscimpl.h>
 18: #include <signal.h>
 19: #if defined(PETSC_HAVE_XMMINTRIN_H)
 20:   #include <xmmintrin.h>
 21: #endif

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

 30: /*@
 31:   PetscFPTrapPush - push a floating point trapping mode, restored using `PetscFPTrapPop()`

 33:   Not Collective

 35:   Input Parameter:
 36: . trap - `PETSC_FP_TRAP_ON` or `PETSC_FP_TRAP_OFF` or any of the values passable to `PetscSetFPTrap()`

 38:   Level: advanced

 40:   Notes:
 41:   This only changes the trapping if the new mode is different than the current mode.

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

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

 53: .seealso: `PetscFPTrapPop()`, `PetscSetFPTrap()`, `PetscDetermineInitialFPTrap()`
 54: @*/
 55: PetscErrorCode PetscFPTrapPush(PetscFPTrap trap)
 56: {
 57:   struct PetscFPTrapLink *link;

 59:   PetscFunctionBegin;
 60:   PetscCall(PetscNew(&link));
 61: #if defined(PETSC_HAVE_THREADSAFETY) && defined(PETSC_HAVE_OPENMP)
 62:   PetscPragmaOMP(critical)
 63: #endif
 64:   {
 65:     link->trapmode = _trapmode;
 66:     link->next     = _trapstack;
 67:     _trapstack     = link;
 68:   }
 69:   if (trap != _trapmode) PetscCall(PetscSetFPTrap(trap));
 70:   PetscFunctionReturn(PETSC_SUCCESS);
 71: }

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

 76:   Not Collective

 78:   Level: advanced

 80: .seealso: `PetscFPTrapPush()`, `PetscSetFPTrap()`, `PetscDetermineInitialFPTrap()`
 81: @*/
 82: PetscErrorCode PetscFPTrapPop(void)
 83: {
 84:   struct PetscFPTrapLink *link;

 86:   PetscFunctionBegin;
 87:   if (_trapstack->trapmode != _trapmode) PetscCall(PetscSetFPTrap(_trapstack->trapmode));
 88: #if defined(PETSC_HAVE_THREADSAFETY) && defined(PETSC_HAVE_OPENMP)
 89:   PetscPragmaOMP(critical)
 90: #endif
 91:   {
 92:     link       = _trapstack;
 93:     _trapstack = _trapstack->next;
 94:   }
 95:   PetscCall(PetscFree(link));
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: /*--------------------------------------- ---------------------------------------------------*/
100: #if defined(PETSC_HAVE_SUN4_STYLE_FPTRAP)
101:   #include <floatingpoint.h>

103: PETSC_EXTERN PetscErrorCode ieee_flags(char *, char *, char *, char **);
104: PETSC_EXTERN PetscErrorCode ieee_handler(char *, char *, sigfpe_handler_type(int, int, struct sigcontext *, char *));

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

120: /* this function gets called if a trap has occurred and been caught */
121: sigfpe_handler_type PetscDefaultFPTrap(int sig, int code, struct sigcontext *scp, char *addr)
122: {
123:   int            err_ind = -1;
124:   PetscErrorCode ierr;

126:   PetscFunctionBegin;
127:   for (int j = 0; error_codes[j].code_no; j++) {
128:     if (error_codes[j].code_no == code) err_ind = j;
129:   }

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

134:   ierr = PetscError(PETSC_COMM_SELF, PETSC_ERR_FP, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
135:   (void)ierr;
136:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

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

144:    Not Collective

146:    Input Parameter:
147: .  flag - values are
148: .vb
149:     PETSC_FP_TRAP_OFF   - do not trap any exceptions
150:     PETSC_FP_TRAP_ON - all exceptions that are possible on the system except underflow
151:     PETSC_FP_TRAP_INDIV - integer divide by zero
152:     PETSC_FP_TRAP_FLTOPERR - improper argument to function, for example with real numbers, the square root of a negative number
153:     PETSC_FP_TRAP_FLTOVF - overflow
154:     PETSC_FP_TRAP_FLTUND - underflow - not trapped by default on most systems
155:     PETSC_FP_TRAP_FLTDIV - floating point divide by zero
156:     PETSC_FP_TRAP_FLTINEX - inexact floating point result
157: .ve

159:    Options Database Key:
160: .  -fp_trap <off,on> - turn on or off trapping of floating point exceptions

162:    Level: advanced

164:    Notes:
165:    Preferred usage is `PetscFPTrapPush()` and `PetscFPTrapPop()` instead of this routine

167:    Currently only `PETSC_FP_TRAP_OFF` and `PETSC_FP_TRAP_ON` are handled. All others are treated as `PETSC_FP_TRAP_ON`.

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

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

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

185: .seealso: `PetscFPTrapPush()`, `PetscFPTrapPop()`, `PetscDetermineInitialFPTrap()`
186: @*/
187: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
188: {
189:   char *out;

191:   PetscFunctionBegin;
192:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
193:   (void)ieee_flags("clear", "exception", "all", &out);
194:   if (flag == PETSC_FP_TRAP_ON) {
195:     /*
196:       To trap more fp exceptions, including underflow, change the line below to
197:       if (ieee_handler("set","all",PetscDefaultFPTrap)) {
198:     */
199:     if (ieee_handler("set", "common", PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
200:   } else if (ieee_handler("clear", "common", PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

202:   _trapmode = flag;
203:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_SUN4_STYLE_FPTRAP\n"));
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

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

210:    Not Collective

212:    Note:
213:    Currently only supported on Linux and macOS. Checks if divide by zero is enable and if so declares that trapping is on.

215:    Level: advanced

217: .seealso: `PetscFPTrapPush()`, `PetscFPTrapPop()`, `PetscDetermineInitialFPTrap()`
218: @*/
219: PetscErrorCode PetscDetermineInitialFPTrap(void)
220: {
221:   PetscFunctionBegin;
222:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: /* -------------------------------------------------------------------------------------------*/
227: #elif defined(PETSC_HAVE_SOLARIS_STYLE_FPTRAP)
228:   #include <sunmath.h>
229:   #include <floatingpoint.h>
230:   #include <siginfo.h>
231:   #include <ucontext.h>

233: static struct {
234:   int   code_no;
235:   char *name;
236: } error_codes[] = {
237:   {FPE_FLTINV, "invalid floating point operand"},
238:   {FPE_FLTRES, "inexact floating point result" },
239:   {FPE_FLTDIV, "division-by-zero"              },
240:   {FPE_FLTUND, "floating point underflow"      },
241:   {FPE_FLTOVF, "floating point overflow"       },
242:   {0,          "unknown error"                 }
243: };
244:   #define SIGPC(scp) (scp->si_addr)

246: void PetscDefaultFPTrap(int sig, siginfo_t *scp, ucontext_t *uap)
247: {
248:   int            err_ind = -1, code = scp->si_code;
249:   PetscErrorCode ierr;

251:   PetscFunctionBegin;
252:   for (int j = 0; error_codes[j].code_no; j++) {
253:     if (error_codes[j].code_no == code) err_ind = j;
254:   }

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

259:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
260:   (void)ierr;
261:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
262: }

264: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
265: {
266:   char *out;

268:   PetscFunctionBegin;
269:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
270:   (void)ieee_flags("clear", "exception", "all", &out);
271:   if (flag == PETSC_FP_TRAP_ON) {
272:     if (ieee_handler("set", "common", (sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floating point handler\n");
273:   } else {
274:     if (ieee_handler("clear", "common", (sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");
275:   }
276:   _trapmode = flag;
277:   PetscCall(PetscInfo(NULL,"Using PETSC_HAVE_SOLARIS_STYLE_FPTRAP\n");
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

281: PetscErrorCode PetscDetermineInitialFPTrap(void)
282: {
283:   PetscFunctionBegin;
284:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /* ------------------------------------------------------------------------------------------*/
289: #elif defined(PETSC_HAVE_IRIX_STYLE_FPTRAP)
290:   #include <sigfpe.h>
291: static struct {
292:   int   code_no;
293:   char *name;
294: } error_codes[] = {
295:   {_INVALID, "IEEE operand error"      },
296:   {_OVERFL,  "floating point overflow" },
297:   {_UNDERFL, "floating point underflow"},
298:   {_DIVZERO, "floating point divide"   },
299:   {0,        "unknown error"           }
300: };
301: void PetscDefaultFPTrap(unsigned exception[], int val[])
302: {
303:   int            err_ind = -1, code = exception[0];
304:   PetscErrorCode ierr;

306:   PetscFunctionBegin;
307:   for (int j = 0; error_codes[j].code_no; j++) {
308:     if (error_codes[j].code_no == code) err_ind = j;
309:   }
310:   if (err_ind >= 0) ierr = (*PetscErrorPrintf)("*** %s occurred ***\n", error_codes[err_ind].name);
311:   else ierr = (*PetscErrorPrintf)("*** floating point error 0x%x occurred ***\n", code);

313:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
314:   (void)ierr;
315:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
316: }

318: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
319: {
320:   PetscFunctionBegin;
321:   if (flag == PETSC_FP_TRAP_ON) handle_sigfpes(_ON, , _EN_UNDERFL | _EN_OVERFL | _EN_DIVZERO | _EN_INVALID, PetscDefaultFPTrap, _ABORT_ON_ERROR, 0);
322:   else handle_sigfpes(_OFF, _EN_UNDERFL | _EN_OVERFL | _EN_DIVZERO | _EN_INVALID, 0, _ABORT_ON_ERROR, 0);
323:   _trapmode = flag;
324:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_IRIX_STYLE_FPTRAP\n"));
325:   PetscFunctionReturn(PETSC_SUCCESS);
326: }

328: PetscErrorCode PetscDetermineInitialFPTrap(void)
329: {
330:   PetscFunctionBegin;
331:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: /*----------------------------------------------- --------------------------------------------*/
336: #elif defined(PETSC_HAVE_RS6000_STYLE_FPTRAP)
337: /* In "fast" mode, floating point traps are imprecise and ignored.
338:    This is the reason for the fptrap(FP_TRAP_SYNC) call */
339: struct sigcontext;
340:   #include <fpxcp.h>
341:   #include <fptrap.h>
342:   #define FPE_FLTOPERR_TRAP (fptrap_t)(0x20000000)
343:   #define FPE_FLTOVF_TRAP   (fptrap_t)(0x10000000)
344:   #define FPE_FLTUND_TRAP   (fptrap_t)(0x08000000)
345:   #define FPE_FLTDIV_TRAP   (fptrap_t)(0x04000000)
346:   #define FPE_FLTINEX_TRAP  (fptrap_t)(0x02000000)

348: static struct {
349:   int   code_no;
350:   char *name;
351: } error_codes[] = {
352:   {FPE_FLTOPERR_TRAP, "IEEE operand error"           },
353:   {FPE_FLTOVF_TRAP,   "floating point overflow"      },
354:   {FPE_FLTUND_TRAP,   "floating point underflow"     },
355:   {FPE_FLTDIV_TRAP,   "floating point divide"        },
356:   {FPE_FLTINEX_TRAP,  "inexact floating point result"},
357:   < {0,                 "unknown error"                }
358: };
359:   #define SIGPC(scp)        (0) /* Info MIGHT be in scp->sc_jmpbuf.jmp_context.iar */
360: /*
361:    For some reason, scp->sc_jmpbuf does not work on the RS6000, even though
362:    it looks like it should from the include definitions.  It is probably
363:    some strange interaction with the "POSIX_SOURCE" that we require.
364: */

366: void PetscDefaultFPTrap(int sig, int code, struct sigcontext *scp)
367: {
368:   int            err_ind, j;
369:   fp_ctx_t       flt_context;
370:   PetscErrorCode ierr;

372:   PetscFunctionBegin;
373:   fp_sh_trap_info(scp, &flt_context);

375:   err_ind = -1;
376:   for (j = 0; error_codes[j].code_no; j++) {
377:     if (error_codes[j].code_no == flt_context.trap) err_ind = j;
378:   }

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

383:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
384:   (void)ierr;
385:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
386: }

388: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
389: {
390:   PetscFunctionBegin;
391:   if (flag == PETSC_FP_TRAP_ON) {
392:     signal(SIGFPE, (void (*)(int))PetscDefaultFPTrap);
393:     fp_trap(FP_TRAP_SYNC);
394:     /* fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW); */
395:     fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW);
396:   } else {
397:     signal(SIGFPE, SIG_DFL);
398:     fp_disable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW);
399:     fp_trap(FP_TRAP_OFF);
400:   }
401:   _trapmode = flag;
402:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_RS6000_STYLE_FPTRAP\n"));
403:   PetscFunctionReturn(PETSC_SUCCESS);
404: }

406: PetscErrorCode PetscDetermineInitialFPTrap(void)
407: {
408:   PetscFunctionBegin;
409:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: /* ------------------------------------------------------------*/
414: #elif defined(PETSC_HAVE_WINDOWS_COMPILERS)
415:   #include <float.h>
416: void PetscDefaultFPTrap(int sig)
417: {
418:   PetscErrorCode ierr;

420:   PetscFunctionBegin;
421:   ierr = (*PetscErrorPrintf)("*** floating point error occurred ***\n");
422:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
423:   (void)ierr;
424:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
425: }

427: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
428: {
429:   unsigned int cw;

431:   PetscFunctionBegin;
432:   if (flag == PETSC_FP_TRAP_ON) {
433:     /* cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW | _EM_UNDERFLOW; */
434:     cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW;
435:     PetscCheck(SIG_ERR != signal(SIGFPE, PetscDefaultFPTrap), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't set floating point handler");
436:   } else {
437:     cw = 0;
438:     PetscCheck(SIG_ERR != signal(SIGFPE, SIG_DFL), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't clear floating point handler");
439:   }
440:   (void)_controlfp(0, cw);
441:   _trapmode = flag;
442:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_WINDOWS_COMPILERS FPTRAP\n"));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: PetscErrorCode PetscDetermineInitialFPTrap(void)
447: {
448:   PetscFunctionBegin;
449:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: /* ------------------------------------------------------------*/
454: #elif defined(PETSC_HAVE_FENV_H) && !defined(__cplusplus)
455:   /*
456:    C99 style floating point environment.

458:    Note that C99 merely specifies how to save, restore, and clear the floating
459:    point environment as well as defining an enumeration of exception codes.  In
460:    particular, C99 does not specify how to make floating point exceptions raise
461:    a signal.  Glibc offers this capability through FE_NOMASK_ENV (or with finer
462:    granularity, feenableexcept()), xmmintrin.h offers _MM_SET_EXCEPTION_MASK().
463: */
464:   #include <fenv.h>
465:   #if defined(PETSC_HAVE_FE_VALUES)
466: typedef struct {
467:   int         code;
468:   const char *name;
469: } FPNode;
470: static const FPNode error_codes[] = {
471:   {FE_DIVBYZERO, "divide by zero"                                 },
472:   {FE_INEXACT,   "inexact floating point result"                  },
473:   {FE_INVALID,   "invalid floating point arguments (domain error)"},
474:   {FE_OVERFLOW,  "floating point overflow"                        },
475:   {FE_UNDERFLOW, "floating point underflow"                       },
476:   {0,            "unknown error"                                  }
477: };
478:   #endif

480: void PetscDefaultFPTrap(int sig)
481: {
482:   #if defined(PETSC_HAVE_FE_VALUES)
483:   const FPNode  *node;
484:   int            code;
485:   PetscBool      matched = PETSC_FALSE;
486:   #endif
487:   PetscErrorCode ierr;

489:   PetscFunctionBegin;
490:   /* Note: While it is possible for the exception state to be preserved by the
491:    * kernel, this seems to be rare which makes the following flag testing almost
492:    * useless.  But on a system where the flags can be preserved, it would provide
493:    * more detail.
494:    */
495:   #if defined(PETSC_HAVE_FE_VALUES)
496:   code = fetestexcept(FE_ALL_EXCEPT);
497:   for (node = &error_codes[0]; node->code; node++) {
498:     if (code & node->code) {
499:       matched = PETSC_TRUE;
500:       ierr    = (*PetscErrorPrintf)("*** floating point error \"%s\" occurred ***\n", node->name);
501:       code &= ~node->code; /* Unset this flag since it has been processed */
502:     }
503:   }
504:   if (!matched || code) { /* If any remaining flags are set, or we didn't process any flags */
505:     ierr = (*PetscErrorPrintf)("*** unknown floating point error occurred ***\n");
506:     ierr = (*PetscErrorPrintf)("The specific exception can be determined by running in a debugger.  When the\n");
507:     ierr = (*PetscErrorPrintf)("debugger traps the signal, the exception can be found with fetestexcept(0x%x)\n", FE_ALL_EXCEPT);
508:     ierr = (*PetscErrorPrintf)("where the result is a bitwise OR of the following flags:\n");
509:     ierr = (*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);
510:   }
511:   #else
512:   ierr = (*PetscErrorPrintf)("*** unknown floating point error occurred ***\n");
513:   #endif

515:   ierr = (*PetscErrorPrintf)("Try option -start_in_debugger\n");
516:   #if PetscDefined(USE_DEBUG)
517:     #if !PetscDefined(HAVE_THREADSAFETY)
518:   ierr = (*PetscErrorPrintf)("likely location of problem given in stack below\n");
519:   ierr = (*PetscErrorPrintf)("---------------------  Stack Frames ------------------------------------\n");
520:   ierr = PetscStackView(PETSC_STDOUT);
521:     #endif
522:   #else
523:   ierr = (*PetscErrorPrintf)("configure using --with-debugging=yes, recompile, link, and run \n");
524:   ierr = (*PetscErrorPrintf)("with -start_in_debugger to get more information on the crash.\n");
525:   #endif
526:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_INITIAL, "trapped floating point error");
527:   (void)ierr;
528:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
529: }

531: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
532: {
533:   PetscFunctionBegin;
534:   if (flag == PETSC_FP_TRAP_ON) {
535:     /* Clear any flags that are currently set so that activating trapping will not immediately call the signal handler. */
536:     PetscCheck(!feclearexcept(FE_ALL_EXCEPT), PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot clear floating point exception flags");
537:   #if defined(FE_NOMASK_ENV) && defined(PETSC_HAVE_FE_VALUES)
538:     /* Could use fesetenv(FE_NOMASK_ENV), but that causes spurious exceptions (like gettimeofday() -> PetscLogDouble). */
539:     /* PetscCheck(feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) != -1,PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot activate floating point exceptions"); */
540:     /* Doesn't work on AArch64 targets. There's a known hardware limitation. Need to detect hardware at configure time? */
541:     PetscCheck(feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW) != -1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot activate floating point exceptions");
542:     PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP with FE_NOMASK_ENV\n"));
543:   #elif defined PETSC_HAVE_XMMINTRIN_H
544:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO);
545:     /* _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW); */
546:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
547:     _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
548:     PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP with PETSC_HAVE_XMMINTRIN_H\n"));
549:   #else
550:     /* C99 does not provide a way to modify the environment so there is no portable way to activate trapping. */
551:     PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_FENV_H FPTRAP\n"));
552:   #endif
553:     PetscCheck(SIG_ERR != signal(SIGFPE, PetscDefaultFPTrap), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't set floating point handler");
554:   } else {
555:     PetscCheck(!fesetenv(FE_DFL_ENV), PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot disable floating point exceptions");
556:     /* can use _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_UNDERFLOW); if PETSC_HAVE_XMMINTRIN_H exists */
557:     PetscCheck(SIG_ERR != signal(SIGFPE, SIG_DFL), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't clear floating point handler");
558:   }
559:   _trapmode = flag;
560:   PetscFunctionReturn(PETSC_SUCCESS);
561: }

563: PetscErrorCode PetscDetermineInitialFPTrap(void)
564: {
565:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
566:   unsigned int flags;
567:   #endif

569:   PetscFunctionBegin;
570:   #if defined(FE_NOMASK_ENV)
571:   flags = fegetexcept();
572:   if (flags & FE_DIVBYZERO) {
573:   #elif defined PETSC_HAVE_XMMINTRIN_H
574:   flags = _MM_GET_EXCEPTION_MASK();
575:   if (!(flags & _MM_MASK_DIV_ZERO)) {
576:   #else
577:   PetscCall(PetscInfo(NULL, "Floating point trapping unknown, assuming off\n"));
578:   PetscFunctionReturn(PETSC_SUCCESS);
579:   #endif
580:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
581:     _trapmode = PETSC_FP_TRAP_ON;
582:     PetscCall(PetscInfo(NULL, "Floating point trapping is on by default %d\n", flags));
583:   } else {
584:     _trapmode = PETSC_FP_TRAP_OFF;
585:     PetscCall(PetscInfo(NULL, "Floating point trapping is off by default %d\n", flags));
586:   }
587:   PetscFunctionReturn(PETSC_SUCCESS);
588:   #endif
589: }

591: /* ------------------------------------------------------------*/
592: #elif defined(PETSC_HAVE_IEEEFP_H)
593:   #include <ieeefp.h>
594: void PetscDefaultFPTrap(int sig)
595: {
596:   PetscErrorCode ierr;

598:   PetscFunctionBegin;
599:   ierr = (*PetscErrorPrintf)("*** floating point error occurred ***\n");
600:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
601:   (void)ierr;
602:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
603: }

605: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
606: {
607:   PetscFunctionBegin;
608:   if (flag == PETSC_FP_TRAP_ON) {
609:   #if defined(PETSC_HAVE_FPRESETSTICKY)
610:     fpresetsticky(fpgetsticky());
611:   #elif defined(PETSC_HAVE_FPSETSTICKY)
612:     fpsetsticky(fpgetsticky());
613:   #endif
614:     fpsetmask(FP_X_INV | FP_X_DZ | FP_X_OFL | FP_X_OFL);
615:     PetscCheck(SIG_ERR != signal(SIGFPE, PetscDefaultFPTrap), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't set floating point handler");
616:   } else {
617:   #if defined(PETSC_HAVE_FPRESETSTICKY)
618:     fpresetsticky(fpgetsticky());
619:   #elif defined(PETSC_HAVE_FPSETSTICKY)
620:     fpsetsticky(fpgetsticky());
621:   #endif
622:     fpsetmask(0);
623:     PetscCheck(SIG_ERR != signal(SIGFPE, SIG_DFL), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't clear floating point handler");
624:   }
625:   _trapmode = flag;
626:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_IEEEFP_H FPTRAP\n"));
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: PetscErrorCode PetscDetermineInitialFPTrap(void)
631: {
632:   PetscFunctionBegin;
633:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
634:   PetscFunctionReturn(PETSC_SUCCESS);
635: }

637: /* -------------------------Default -----------------------------------*/
638: #else

640: static void PetscDefaultFPTrap(int sig)
641: {
642:   PetscErrorCode ierr;

644:   PetscFunctionBegin;
645:   ierr = (*PetscErrorPrintf)("*** floating point error occurred ***\n");
646:   ierr = PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
647:   (void)ierr;
648:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
649: }

651: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
652: {
653:   PetscFunctionBegin;
654:   if (flag == PETSC_FP_TRAP_ON) {
655:     if (SIG_ERR == signal(SIGFPE, PetscDefaultFPTrap)) PetscCall((*PetscErrorPrintf)("Can't set floatingpoint handler\n"));
656:   } else if (SIG_ERR == signal(SIGFPE, SIG_DFL)) PetscCall((*PetscErrorPrintf)("Can't clear floatingpoint handler\n"));

658:   _trapmode = flag;
659:   PetscCall(PetscInfo(NULL, "Using default FPTRAP\n"));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: PetscErrorCode PetscDetermineInitialFPTrap(void)
664: {
665:   PetscFunctionBegin;
666:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }
669: #endif