Actual source code: fp.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  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: *
  6: */

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

 19: #include <petscsys.h>
 20: #include <signal.h>

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

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

 32:    Not Collective

 34:    Input Arguments:
 35: .    trap - PETSC_FP_TRAP_ON or PETSC_FP_TRAP_OFF

 37:    Level: advanced

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

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

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

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

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

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

 69:    Not Collective

 71:    Level: advanced

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

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

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

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

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

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

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

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

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

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

128:    Not Collective

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

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

136:    Level: advanced

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

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

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


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

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

176:   _trapmode = flag;
177:   return(0);
178: }

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

183:    Not Collective

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

188:    Level: advanced

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

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

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

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

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

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

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

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

236: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
237: {
238:   char *out;

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

252: PetscErrorCode  PetscDetermineInitialFPTrap(void)
253: {

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

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

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

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

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

297: PetscErrorCode  PetscDetermineInitialFPTrap(void)
298: {

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

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

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

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

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

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

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

341: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
342: {
343:   char *out;

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

357: PetscErrorCode  PetscDetermineInitialFPTrap(void)
358: {

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

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

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

394: void PetscDefaultFPTrap(int sig,int code,struct sigcontext *scp)
395: {
397:   int            err_ind,j;
398:   fp_ctx_t       flt_context;

401:   fp_sh_trap_info(scp,&flt_context);

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

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

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

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

440: PetscErrorCode  PetscDetermineInitialFPTrap(void)
441: {

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

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

460: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
461: {
462:   unsigned int cw;

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

477: PetscErrorCode  PetscDetermineInitialFPTrap(void)
478: {

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

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

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

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

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

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

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

579: PetscErrorCode  PetscDetermineInitialFPTrap(void)
580: {
581: #if defined(FE_NOMASK_ENV) || defined PETSC_HAVE_XMMINTRIN_H
582:   unsigned int   flags;
583: #endif

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

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

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

643: PetscErrorCode  PetscDetermineInitialFPTrap(void)
644: {

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

652: /* -------------------------Default -----------------------------------*/
653: #else

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

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

670:   _trapmode = on;
671:   return(0);
672: }

674: PetscErrorCode  PetscDetermineInitialFPTrap(void)
675: {

679:   PetscInfo(NULL,"Unable to determine initial floating point trapping. Assuming it is off\n");
680:   return(0);
681: }
682: #endif