From 6071e508e889fc582b574ca71233841c3b4f9833 Mon Sep 17 00:00:00 2001
From: Benjamin Brands <benjamin.brands@fau.de>
Date: Thu, 12 Apr 2018 20:46:41 +0200
Subject: [PATCH] remedy for exceptions from within p_syevr/_syevr

---
 include/deal.II/lac/lapack_templates.h    | 32 +++++++++++++++++++++++
 include/deal.II/lac/scalapack.templates.h | 32 +++++++++++++++++++++++
 2 files changed, 64 insertions(+)

diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h
index e3c11a310e..8606ce83eb 100644
--- 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 a21ed45d23..05dfbfb8ea 100644
--- 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
-- 
2.39.5