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: .vb
 47:        gfortran -ffpe-trap=invalid,zero,overflow,underflow,denormal
 48:        ifort -fpe0
 49: .ve

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

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

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

 74:    Not Collective

 76:    Level: advanced

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

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

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

101: PETSC_EXTERN PetscErrorCode ieee_flags(char *, char *, char *, char **);
102: PETSC_EXTERN PetscErrorCode ieee_handler(char *, char *, sigfpe_handler_type(int, int, struct sigcontext *, char *));

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

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

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

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

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

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

142:    Not Collective

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

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

160:    Level: advanced

162:    Notes:
163:    Preferred usage is `PetscFPTrapPush()` and `PetscFPTrapPop()` instead of this routine

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

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

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

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

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

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

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

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

208:    Not Collective

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

213:    Level: advanced

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

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

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

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

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

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

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

262: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
263: {
264:   char *out;

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

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

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

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

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

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

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

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

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

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

370:   PetscFunctionBegin;
371:   fp_sh_trap_info(scp, &flt_context);

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

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

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

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

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

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

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

425: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
426: {
427:   unsigned int cw;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

635: /* -------------------------Default -----------------------------------*/
636: #else

638: void PetscDefaultFPTrap(int sig)
639: {
640:   PetscErrorCode ierr;

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

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

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

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