Actual source code: fp.c

petsc-3.12.5 2020-03-29
Report Typos and Errors

  2: /*
  3: *   IEEE error handler for all machines. Since each machine 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 */
 27: static struct PetscFPTrapLink *_trapstack;                   /* Any pushed states of _trapmode */

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

 32:    Not Collective

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

 37:    Level: advanced

 39: .seealso: PetscFPTrapPop(), PetscSetFPTrap()
 40: @*/
 41: PetscErrorCode PetscFPTrapPush(PetscFPTrap trap)
 42: {
 43:   PetscErrorCode         ierr;
 44:   struct PetscFPTrapLink *link;

 47:   PetscNew(&link);
 48:   link->trapmode = _trapmode;
 49:   link->next     = _trapstack;
 50:   _trapstack     = link;
 51:   if (trap != _trapmode) {PetscSetFPTrap(trap);}
 52:   return(0);
 53: }

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

 58:    Not Collective

 60:    Level: advanced

 62: .seealso: PetscFPTrapPush(), PetscSetFPTrap()
 63: @*/
 64: PetscErrorCode PetscFPTrapPop(void)
 65: {
 66:   PetscErrorCode         ierr;
 67:   struct PetscFPTrapLink *link;

 70:   if (_trapstack->trapmode != _trapmode) {PetscSetFPTrap(_trapstack->trapmode);}
 71:   link       = _trapstack;
 72:   _trapstack = _trapstack->next;
 73:   PetscFree(link);
 74:   return(0);
 75: }

 77: /*--------------------------------------- ---------------------------------------------------*/
 78: #if defined(PETSC_HAVE_SUN4_STYLE_FPTRAP)
 79: #include <floatingpoint.h>

 81: PETSC_EXTERN PetscErrorCode ieee_flags(char*,char*,char*,char**);
 82: PETSC_EXTERN PetscErrorCode ieee_handler(char*,char*,sigfpe_handler_type(int,int,struct sigcontext*,char*));

 84: static struct { int code_no; char *name; } error_codes[] = {
 85:   { FPE_INTDIV_TRAP    ,"integer divide" },
 86:   { FPE_FLTOPERR_TRAP  ,"IEEE operand error" },
 87:   { FPE_FLTOVF_TRAP    ,"floating point overflow" },
 88:   { FPE_FLTUND_TRAP    ,"floating point underflow" },
 89:   { FPE_FLTDIV_TRAP    ,"floating pointing divide" },
 90:   { FPE_FLTINEX_TRAP   ,"inexact floating point result" },
 91:   { 0                  ,"unknown error" }
 92: };
 93: #define SIGPC(scp) (scp->sc_pc)

 95: sigfpe_handler_type PetscDefaultFPTrap(int sig,int code,struct sigcontext *scp,char *addr)
 96: {
 98:   int            err_ind = -1,j;

101:   for (j = 0; error_codes[j].code_no; j++) {
102:     if (error_codes[j].code_no == code) err_ind = j;
103:   }

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

108:   PetscError(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
109:   MPI_Abort(PETSC_COMM_WORLD,0);
110:   return(0);
111: }

113: /*@
114:    PetscSetFPTrap - Enables traps/exceptions on common floating point errors.
115:                     This option may not work on certain machines.

117:    Not Collective

119:    Input Parameters:
120: .  flag - PETSC_FP_TRAP_ON, PETSC_FP_TRAP_OFF.

122:    Options Database Keys:
123: .  -fp_trap - Activates floating point trapping

125:    Level: advanced

127:    Description:
128:    On systems that support it, this routine causes floating point
129:    overflow, divide-by-zero, and invalid-operand (e.g., a NaN) to
130:    cause a message to be printed and the program to exit.

132:    Note:
133:    On many common systems including x86 and x86-64 Linux, the floating
134:    point exception state is not preserved from the location where the trap
135:    occurred through to the signal handler.  In this case, the signal handler
136:    will just say that an unknown floating point exception occurred and which
137:    function it occurred in.  If you run with -fp_trap in a debugger, it will
138:    break on the line where the error occurred.  You can check which
139:    exception occurred using fetestexcept(FE_ALL_EXCEPT).  See fenv.h
140:    (usually at /usr/include/bits/fenv.h) for the enum values on your system.

142:    Caution:
143:    On certain machines, in particular the IBM rs6000, floating point
144:    trapping is VERY slow!


147: .seealso: PetscFPTrapPush(), PetscFPTrapPop()
148: @*/
149: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
150: {
151:   char *out;

154:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
155:   (void) ieee_flags("clear","exception","all",&out);
156:   if (flag == PETSC_FP_TRAP_ON) {
157:     /*
158:       To trap more fp exceptions, including underflow, change the line below to
159:       if (ieee_handler("set","all",PetscDefaultFPTrap)) {
160:     */
161:     if (ieee_handler("set","common",PetscDefaultFPTrap))        (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
162:   } else if (ieee_handler("clear","common",PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

164:   _trapmode = flag;
165:   return(0);
166: }

168: /* -------------------------------------------------------------------------------------------*/
169: #elif defined(PETSC_HAVE_SOLARIS_STYLE_FPTRAP)
170: #include <sunmath.h>
171: #include <floatingpoint.h>
172: #include <siginfo.h>
173: #include <ucontext.h>

175: static struct { int code_no; char *name; } error_codes[] = {
176:   { FPE_FLTINV,"invalid floating point operand"},
177:   { FPE_FLTRES,"inexact floating point result"},
178:   { FPE_FLTDIV,"division-by-zero"},
179:   { FPE_FLTUND,"floating point underflow"},
180:   { FPE_FLTOVF,"floating point overflow"},
181:   { 0,         "unknown error"}
182: };
183: #define SIGPC(scp) (scp->si_addr)

185: void PetscDefaultFPTrap(int sig,siginfo_t *scp,ucontext_t *uap)
186: {
187:   int            err_ind,j,code = scp->si_code;

191:   err_ind = -1;
192:   for (j = 0; error_codes[j].code_no; j++) {
193:     if (error_codes[j].code_no == code) err_ind = j;
194:   }

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

199:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
200:   MPI_Abort(PETSC_COMM_WORLD,0);
201: }

203: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
204: {
205:   char *out;

208:   /* Clear accumulated exceptions.  Used to suppress meaningless messages from f77 programs */
209:   (void) ieee_flags("clear","exception","all",&out);
210:   if (flag == PETSC_FP_TRAP_ON) {
211:     if (ieee_handler("set","common",(sigfpe_handler_type)PetscDefaultFPTrap))        (*PetscErrorPrintf)("Can't set floating point handler\n");
212:   } else if (ieee_handler("clear","common",(sigfpe_handler_type)PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");
213:   _trapmode = flag;
214:   return(0);
215: }

217: /* ------------------------------------------------------------------------------------------*/

219: #elif defined(PETSC_HAVE_IRIX_STYLE_FPTRAP)
220: #include <sigfpe.h>
221: static struct { int code_no; char *name; } error_codes[] = {
222:   { _INVALID   ,"IEEE operand error" },
223:   { _OVERFL    ,"floating point overflow" },
224:   { _UNDERFL   ,"floating point underflow" },
225:   { _DIVZERO   ,"floating point divide" },
226:   { 0          ,"unknown error" }
227: } ;
228: void PetscDefaultFPTrap(unsigned exception[],int val[])
229: {
230:   int err_ind,j,code;

233:   code    = exception[0];
234:   err_ind = -1;
235:   for (j = 0; error_codes[j].code_no; j++) {
236:     if (error_codes[j].code_no == code) err_ind = j;
237:   }
238:   if (err_ind >= 0) (*PetscErrorPrintf)("*** %s occurred ***\n",error_codes[err_ind].name);
239:   else              (*PetscErrorPrintf)("*** floating point error 0x%x occurred ***\n",code);

241:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
242:   MPI_Abort(PETSC_COMM_WORLD,0);
243: }

245: PetscErrorCode PetscSetFPTrap(PetscFPTrap flag)
246: {
248:   if (flag == PETSC_FP_TRAP_ON) handle_sigfpes(_ON,_EN_OVERFL|_EN_DIVZERO|_EN_INVALID,PetscDefaultFPTrap,_ABORT_ON_ERROR,0);
249:   else                          handle_sigfpes(_OFF,_EN_OVERFL|_EN_DIVZERO|_EN_INVALID,0,_ABORT_ON_ERROR,0);

251:   _trapmode = flag;
252:   return(0);
253: }
254: /*----------------------------------------------- --------------------------------------------*/
255: /* In "fast" mode, floating point traps are imprecise and ignored.
256:    This is the reason for the fptrap(FP_TRAP_SYNC) call */
257: #elif defined(PETSC_HAVE_RS6000_STYLE_FPTRAP)
258: struct sigcontext;
259: #include <fpxcp.h>
260: #include <fptrap.h>
261: #define FPE_FLTOPERR_TRAP (fptrap_t)(0x20000000)
262: #define FPE_FLTOVF_TRAP   (fptrap_t)(0x10000000)
263: #define FPE_FLTUND_TRAP   (fptrap_t)(0x08000000)
264: #define FPE_FLTDIV_TRAP   (fptrap_t)(0x04000000)
265: #define FPE_FLTINEX_TRAP  (fptrap_t)(0x02000000)

267: static struct { int code_no; char *name; } error_codes[] = {
268:   {FPE_FLTOPERR_TRAP   ,"IEEE operand error" },
269:   { FPE_FLTOVF_TRAP    ,"floating point overflow" },
270:   { FPE_FLTUND_TRAP    ,"floating point underflow" },
271:   { FPE_FLTDIV_TRAP    ,"floating point divide" },
272:   { FPE_FLTINEX_TRAP   ,"inexact floating point result" },
273:   { 0                  ,"unknown error" }
274: } ;
275: #define SIGPC(scp) (0) /* Info MIGHT be in scp->sc_jmpbuf.jmp_context.iar */
276: /*
277:    For some reason, scp->sc_jmpbuf does not work on the RS6000, even though
278:    it looks like it should from the include definitions.  It is probably
279:    some strange interaction with the "POSIX_SOURCE" that we require.
280: */

282: void PetscDefaultFPTrap(int sig,int code,struct sigcontext *scp)
283: {
285:   int            err_ind,j;
286:   fp_ctx_t       flt_context;

289:   fp_sh_trap_info(scp,&flt_context);

291:   err_ind = -1;
292:   for (j = 0; error_codes[j].code_no; j++) {
293:     if (error_codes[j].code_no == flt_context.trap) err_ind = j;
294:   }

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

299:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
300:   MPI_Abort(PETSC_COMM_WORLD,0);
301: }

303: PetscErrorCode PetscSetFPTrap(PetscFPTrap on)
304: {
306:   if (on == PETSC_FP_TRAP_ON) {
307:     signal(SIGFPE,(void (*)(int))PetscDefaultFPTrap);
308:     fp_trap(FP_TRAP_SYNC);
309:     fp_enable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW);
310:     /* fp_enable(mask) for individual traps.  Values are:
311:        TRP_INVALID
312:        TRP_DIV_BY_ZERO
313:        TRP_OVERFLOW
314:        TRP_UNDERFLOW
315:        TRP_INEXACT
316:        Can OR then together.
317:        fp_enable_all(); for all traps.
318:     */
319:   } else {
320:     signal(SIGFPE,SIG_DFL);
321:     fp_disable(TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW);
322:     fp_trap(FP_TRAP_OFF);
323:   }
324:   _trapmode = on;
325:   return(0);
326: }

328: #elif defined(PETSC_HAVE_FENV_H) && !defined(__cplusplus)
329: /*
330:    C99 style floating point environment.

332:    Note that C99 merely specifies how to save, restore, and clear the floating
333:    point environment as well as defining an enumeration of exception codes.  In
334:    particular, C99 does not specify how to make floating point exceptions raise
335:    a signal.  Glibc offers this capability through FE_NOMASK_ENV (or with finer
336:    granularity, feenableexcept()), xmmintrin.h offers _MM_SET_EXCEPTION_MASK().
337: */
338: #include <fenv.h>
339: typedef struct {int code; const char *name;} FPNode;
340: static const FPNode error_codes[] = {
341:   {FE_DIVBYZERO,"divide by zero"},
342:   {FE_INEXACT,  "inexact floating point result"},
343:   {FE_INVALID,  "invalid floating point arguments (domain error)"},
344:   {FE_OVERFLOW, "floating point overflow"},
345:   {FE_UNDERFLOW,"floating point underflow"},
346:   {0           ,"unknown error"}
347: };

349: void PetscDefaultFPTrap(int sig)
350: {
351:   const FPNode *node;
352:   int          code;
353:   PetscBool    matched = PETSC_FALSE;

356:   /* Note: While it is possible for the exception state to be preserved by the
357:    * kernel, this seems to be rare which makes the following flag testing almost
358:    * useless.  But on a system where the flags can be preserved, it would provide
359:    * more detail.
360:    */
361:   code = fetestexcept(FE_ALL_EXCEPT);
362:   for (node=&error_codes[0]; node->code; node++) {
363:     if (code & node->code) {
364:       matched = PETSC_TRUE;
365:       (*PetscErrorPrintf)("*** floating point error \"%s\" occurred ***\n",node->name);
366:       code &= ~node->code; /* Unset this flag since it has been processed */
367:     }
368:   }
369:   if (!matched || code) { /* If any remaining flags are set, or we didn't process any flags */
370:     (*PetscErrorPrintf)("*** unknown floating point error occurred ***\n");
371:     (*PetscErrorPrintf)("The specific exception can be determined by running in a debugger.  When the\n");
372:     (*PetscErrorPrintf)("debugger traps the signal, the exception can be found with fetestexcept(0x%x)\n",FE_ALL_EXCEPT);
373:     (*PetscErrorPrintf)("where the result is a bitwise OR of the following flags:\n");
374:     (*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);
375:   }

377:   (*PetscErrorPrintf)("Try option -start_in_debugger\n");
378: #if defined(PETSC_USE_DEBUG)
379:   if (!PetscStackActive()) (*PetscErrorPrintf)("  or try option -log_stack\n");
380:   else {
381:     (*PetscErrorPrintf)("likely location of problem given in stack below\n");
382:     (*PetscErrorPrintf)("---------------------  Stack Frames ------------------------------------\n");
383:     PetscStackView(PETSC_STDOUT);
384:   }
385: #endif
386: #if !defined(PETSC_USE_DEBUG)
387:   (*PetscErrorPrintf)("configure using --with-debugging=yes, recompile, link, and run \n");
388:   (*PetscErrorPrintf)("with -start_in_debugger to get more information on the crash.\n");
389: #endif
390:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_INITIAL,"trapped floating point error");
391:   MPI_Abort(PETSC_COMM_WORLD,0);
392: }

394: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
395: {
397:   if (on == PETSC_FP_TRAP_ON) {
398:     /* Clear any flags that are currently set so that activating trapping will not immediately call the signal handler. */
399:     if (feclearexcept(FE_ALL_EXCEPT)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot clear floating point exception flags\n");
400: #if defined FE_NOMASK_ENV
401:     /* We could use fesetenv(FE_NOMASK_ENV), but that causes spurious exceptions (like gettimeofday() -> PetscLogDouble). */
402:     if (feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot activate floating point exceptions\n");
403: #elif defined PETSC_HAVE_XMMINTRIN_H
404:     _MM_SET_EXCEPTION_MASK(_MM_MASK_INEXACT | _MM_MASK_UNDERFLOW);
405: #else
406:     /* C99 does not provide a way to modify the environment so there is no portable way to activate trapping. */
407: #endif
408:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't set floating point handler\n");
409:   } else {
410:     if (fesetenv(FE_DFL_ENV)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot disable floating point exceptions");
411:     if (SIG_ERR == signal(SIGFPE,SIG_DFL)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can't clear floating point handler\n");
412:   }
413:   _trapmode = on;
414:   return(0);
415: }

417: /* -------------------------Default -----------------------------------*/
418: #else

420: void PetscDefaultFPTrap(int sig)
421: {
423:   (*PetscErrorPrintf)("*** floating point error occurred ***\n");
424:   PetscError(PETSC_COMM_SELF,0,"User provided function","Unknown file",PETSC_ERR_FP,PETSC_ERROR_REPEAT,"floating point error");
425:   MPI_Abort(PETSC_COMM_WORLD,0);
426: }

428: PetscErrorCode  PetscSetFPTrap(PetscFPTrap on)
429: {
431:   if (on == PETSC_FP_TRAP_ON) {
432:     if (SIG_ERR == signal(SIGFPE,PetscDefaultFPTrap)) (*PetscErrorPrintf)("Can't set floatingpoint handler\n");
433:   } else if (SIG_ERR == signal(SIGFPE,SIG_DFL))       (*PetscErrorPrintf)("Can't clear floatingpoint handler\n");

435:   _trapmode = on;
436:   return(0);
437: }
438: #endif