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;

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

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

    133:   (void)PetscError(PETSC_COMM_SELF, PETSC_ERR_FP, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    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;

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

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

    256:   PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    257:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    258: }

    260: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
    261: {
    262:   char *out;

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

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

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

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

    308:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    309:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    310: }

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

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

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

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

    360: void PetscDefaultFPTrap(int sig, int code, struct sigcontext *scp)
    361: {
    362:   int      err_ind, j;
    363:   fp_ctx_t flt_context;

    365:   PetscFunctionBegin;
    366:   fp_sh_trap_info(scp, &flt_context);

    368:   err_ind = -1;
    369:   for (j = 0; error_codes[j].code_no; j++) {
    370:     if (error_codes[j].code_no == flt_context.trap) err_ind = j;
    371:   }

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

    376:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    377:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    378: }

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

    398: PetscErrorCode PetscDetermineInitialFPTrap(void)
    399: {
    400:   PetscFunctionBegin;
    401:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
    402:   PetscFunctionReturn(PETSC_SUCCESS);
    403: }

    405: /* ------------------------------------------------------------*/
    406: #elif defined(PETSC_HAVE_WINDOWS_COMPILERS)
    407:   #include <float.h>
    408: void PetscDefaultFPTrap(int sig)
    409: {
    410:   PetscFunctionBegin;
    411:   (void)(*PetscErrorPrintf)("*** floating point error occurred ***\n");
    412:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    413:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    414: }

    416: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
    417: {
    418:   unsigned int cw;

    420:   PetscFunctionBegin;
    421:   if (flag == PETSC_FP_TRAP_ON) {
    422:     /* cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW | _EM_UNDERFLOW; */
    423:     cw = _EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW;
    424:     PetscCheck(SIG_ERR != signal(SIGFPE, PetscDefaultFPTrap), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't set floating point handler");
    425:   } else {
    426:     cw = 0;
    427:     PetscCheck(SIG_ERR != signal(SIGFPE, SIG_DFL), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't clear floating point handler");
    428:   }
    429:   (void)_controlfp(0, cw);
    430:   _trapmode = flag;
    431:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_WINDOWS_COMPILERS FPTRAP\n"));
    432:   PetscFunctionReturn(PETSC_SUCCESS);
    433: }

    435: PetscErrorCode PetscDetermineInitialFPTrap(void)
    436: {
    437:   PetscFunctionBegin;
    438:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
    439:   PetscFunctionReturn(PETSC_SUCCESS);
    440: }

    442: /* ------------------------------------------------------------*/
    443: #elif defined(PETSC_HAVE_FENV_H) && !defined(__cplusplus)
    444:   /*
    445:    C99 style floating point environment.

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

    469: void PetscDefaultFPTrap(int sig)
    470: {
    471:   #if defined(PETSC_HAVE_FE_VALUES)
    472:   const FPNode *node;
    473:   int           code;
    474:   PetscBool     matched = PETSC_FALSE;
    475:   #endif

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

    503:   (void)(*PetscErrorPrintf)("Try option -start_in_debugger\n");
    504:   #if PetscDefined(USE_DEBUG)
    505:     #if !PetscDefined(HAVE_THREADSAFETY)
    506:   (void)(*PetscErrorPrintf)("likely location of problem given in stack below\n");
    507:   (void)(*PetscErrorPrintf)("---------------------  Stack Frames ------------------------------------\n");
    508:   (void)PetscStackView(PETSC_STDOUT);
    509:     #endif
    510:   #else
    511:   (void)(*PetscErrorPrintf)("configure using --with-debugging=yes, recompile, link, and run \n");
    512:   (void)(*PetscErrorPrintf)("with -start_in_debugger to get more information on the crash.\n");
    513:   #endif
    514:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_INITIAL, "trapped floating point error");
    515:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    516: }

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

    550: PetscErrorCode PetscDetermineInitialFPTrap(void)
    551: {
    552:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
    553:   unsigned int flags;
    554:   #endif

    556:   PetscFunctionBegin;
    557:   #if defined(FE_NOMASK_ENV)
    558:   flags = fegetexcept();
    559:   if (flags & FE_DIVBYZERO) {
    560:   #elif defined PETSC_HAVE_XMMINTRIN_H
    561:   flags = _MM_GET_EXCEPTION_MASK();
    562:   if (!(flags & _MM_MASK_DIV_ZERO)) {
    563:   #else
    564:   PetscCall(PetscInfo(NULL, "Floating point trapping unknown, assuming off\n"));
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566:   #endif
    567:   #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
    568:     _trapmode = PETSC_FP_TRAP_ON;
    569:     PetscCall(PetscInfo(NULL, "Floating point trapping is on by default %d\n", flags));
    570:   } else {
    571:     _trapmode = PETSC_FP_TRAP_OFF;
    572:     PetscCall(PetscInfo(NULL, "Floating point trapping is off by default %d\n", flags));
    573:   }
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575:   #endif
    576: }

    578: /* ------------------------------------------------------------*/
    579: #elif defined(PETSC_HAVE_IEEEFP_H)
    580:   #include <ieeefp.h>
    581: void PetscDefaultFPTrap(int sig)
    582: {
    583:   PetscFunctionBegin;
    584:   (void)(*PetscErrorPrintf)("*** floating point error occurred ***\n");
    585:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    586:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    587: }

    589: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
    590: {
    591:   PetscFunctionBegin;
    592:   if (flag == PETSC_FP_TRAP_ON) {
    593:   #if defined(PETSC_HAVE_FPRESETSTICKY)
    594:     fpresetsticky(fpgetsticky());
    595:   #elif defined(PETSC_HAVE_FPSETSTICKY)
    596:     fpsetsticky(fpgetsticky());
    597:   #endif
    598:     fpsetmask(FP_X_INV | FP_X_DZ | FP_X_OFL | FP_X_OFL);
    599:     PetscCheck(SIG_ERR != signal(SIGFPE, PetscDefaultFPTrap), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't set floating point handler");
    600:   } else {
    601:   #if defined(PETSC_HAVE_FPRESETSTICKY)
    602:     fpresetsticky(fpgetsticky());
    603:   #elif defined(PETSC_HAVE_FPSETSTICKY)
    604:     fpsetsticky(fpgetsticky());
    605:   #endif
    606:     fpsetmask(0);
    607:     PetscCheck(SIG_ERR != signal(SIGFPE, SIG_DFL), PETSC_COMM_SELF, PETSC_ERR_LIB, "Can't clear floating point handler");
    608:   }
    609:   _trapmode = flag;
    610:   PetscCall(PetscInfo(NULL, "Using PETSC_HAVE_IEEEFP_H FPTRAP\n"));
    611:   PetscFunctionReturn(PETSC_SUCCESS);
    612: }

    614: PetscErrorCode PetscDetermineInitialFPTrap(void)
    615: {
    616:   PetscFunctionBegin;
    617:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
    618:   PetscFunctionReturn(PETSC_SUCCESS);
    619: }

    621: /* -------------------------Default -----------------------------------*/
    622: #else

    624: static void PetscDefaultFPTrap(int sig)
    625: {
    626:   PetscFunctionBegin;
    627:   (void)(*PetscErrorPrintf)("*** floating point error occurred ***\n");
    628:   (void)PetscError(PETSC_COMM_SELF, 0, NULL, NULL, PETSC_ERR_FP, PETSC_ERROR_REPEAT, "floating point error");
    629:   PETSCABORT(MPI_COMM_WORLD, PETSC_ERR_FP);
    630: }

    632: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
    633: {
    634:   PetscFunctionBegin;
    635:   if (flag == PETSC_FP_TRAP_ON) {
    636:     if (SIG_ERR == signal(SIGFPE, PetscDefaultFPTrap)) PetscCall((*PetscErrorPrintf)("Can't set floatingpoint handler\n"));
    637:   } else if (SIG_ERR == signal(SIGFPE, SIG_DFL)) PetscCall((*PetscErrorPrintf)("Can't clear floatingpoint handler\n"));

    639:   _trapmode = flag;
    640:   PetscCall(PetscInfo(NULL, "Using default FPTRAP\n"));
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }

    644: PetscErrorCode PetscDetermineInitialFPTrap(void)
    645: {
    646:   PetscFunctionBegin;
    647:   PetscCall(PetscInfo(NULL, "Unable to determine initial floating point trapping. Assuming it is off\n"));
    648:   PetscFunctionReturn(PETSC_SUCCESS);
    649: }
    650: #endif