remedy for exceptions from within p_syevr/_syevr

diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h
index e3c11a310e06508c8ebf8b681653bb28ad5cd473..8606ce83ebc5f0da62b6a97a7132afca5fd0ce38 100644 (file)
--- a/include/deal.II/lac/lapack_templates.h
+++ b/include/deal.II/lac/lapack_templates.h
@@ -20,6 +20,10 @@
 #include <deal.II/base/config.h>
 #include <deal.II/lac/lapack_support.h>
 
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+#include <cfenv>
+#endif
+
 extern "C"
 {
 // vector update of the form y += alpha*x with a scalar, x,y vectors
@@ -1401,7 +1405,21 @@ syevr(const char *jobz,
       types::blas_int *liwork,
       types::blas_int *info)
 {
+  /*
+   * Netlib and Atlas Lapack perform floating point tests (e.g. divide-by-zero) within the call to dsyevr
+   * causing floating point exceptions to be thrown (at least in debug mode). Therefore, we wrap the calls
+   * to dsyevr into the following code to supress the exception.
+   */
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fenv_t fp_exceptions;
+  feholdexcept(&fp_exceptions);
+#endif
+
   dsyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info);
+
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fesetenv(&fp_exceptions);
+#endif
 }
 
 inline void
@@ -1427,7 +1445,21 @@ syevr(const char *jobz,
       types::blas_int *liwork,
       types::blas_int *info)
 {
+  /*
+   * Netlib and Atlas Lapack perform floating point tests (e.g. divide-by-zero) within the call to ssyevr
+   * causing floating point exceptions to be thrown (at least in debug mode). Therefore, we wrap the calls
+   * to ssyevr into the following code to supress the exception.
+   */
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fenv_t fp_exceptions;
+  feholdexcept(&fp_exceptions);
+#endif
+
   ssyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info);
+
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fesetenv(&fp_exceptions);
+#endif
 }
 #endif
 

diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h
index a21ed45d239091ec2e84d824d339f57fb7d549ae..05dfbfb8eaeb851ab616d87b7512d0a78b389681 100644 (file)
--- a/include/deal.II/lac/scalapack.templates.h
+++ b/include/deal.II/lac/scalapack.templates.h
@@ -24,6 +24,10 @@
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/mpi.templates.h>
 
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+#include <cfenv>
+#endif
+
 // useful examples:
 // https://stackoverflow.com/questions/14147705/cholesky-decomposition-scalapack-error/14203864
 // http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=139   // second post by Julien Langou
@@ -1902,7 +1906,21 @@ inline void psyevr(const char *jobz,
                    int *liwork,
                    int *info)
 {
+  /*
+   * Netlib ScaLAPACK performs floating point tests (e.g. divide-by-zero) within the call to pdsyevr
+   * causing floating point exceptions to be thrown (at least in debug mode). Therefore, we wrap the calls
+   * to pdsyevr into the following code to supress the exception.
+   */
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fenv_t fp_exceptions;
+  feholdexcept(&fp_exceptions);
+#endif
+
   pdsyevr_(jobz,range,uplo,n,A,IA,JA,DESCA,VL,VU,IL,IU,m,nz,w,Z,IZ,JZ,DESCZ,work,lwork,iwork,liwork,info);
+
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fesetenv(&fp_exceptions);
+#endif
 }
 
 inline void psyevr(const char *jobz,
@@ -1930,7 +1948,21 @@ inline void psyevr(const char *jobz,
                    int *liwork,
                    int *info)
 {
+  /*
+   * Netlib ScaLAPACK performs floating point tests (e.g. divide-by-zero) within the call to pssyevr
+   * causing floating point exceptions to be thrown (at least in debug mode). Therefore, we wrap the calls
+   * to pssyevr into the following code to supress the exception.
+   */
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fenv_t fp_exceptions;
+  feholdexcept(&fp_exceptions);
+#endif
+
   pssyevr_(jobz,range,uplo,n,A,IA,JA,DESCA,VL,VU,IL,IU,m,nz,w,Z,IZ,JZ,DESCZ,work,lwork,iwork,liwork,info);
+
+#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
+  fesetenv(&fp_exceptions);
+#endif
 }
 
 #endif // DEAL_II_WITH_SCALAPACK