From 3d713fc1a7e81e43fc8a5270350103cba0552232 Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 22 Aug 2001 11:53:44 +0000
Subject: [PATCH] Check pivot size relative to other matrix entries.

git-svn-id: https://svn.dealii.org/trunk@4904 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/lac/include/lac/full_matrix.h         |  5 +-
 .../lac/include/lac/full_matrix.templates.h   | 54 +++++++++++++------
 2 files changed, 41 insertions(+), 18 deletions(-)

diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h
index dc78f223ec..0f7eb948fd 100644
--- a/deal.II/lac/include/lac/full_matrix.h
+++ b/deal.II/lac/include/lac/full_matrix.h
@@ -613,7 +613,10 @@ class FullMatrix : public vector2d<number>
 				     /**
 				      * Exception
 				      */
-    DeclException0 (ExcNotRegular);
+    DeclException1 (ExcNotRegular,
+		    double,
+		    << "The maximal pivot is " << arg1
+		    << ", which is below the threshold. The matrix may be singular.");
 				     /**
 				      * Exception
 				      */
diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h
index 7a7e4894e8..b2b9b6e6b6 100644
--- a/deal.II/lac/include/lac/full_matrix.templates.h
+++ b/deal.II/lac/include/lac/full_matrix.templates.h
@@ -1340,17 +1340,33 @@ FullMatrix<number>::gauss_jordan()
   
 				   // Gauss-Jordan-Algorithmus
 				   // cf. Stoer I (4th Edition) p. 153
-  std::vector<unsigned int> p(n());
-
-  for (unsigned int i=0; i<n(); ++i)
+  const unsigned int N = n();
+
+				   // first get an estimate of the
+				   // size of the elements of this
+				   // matrix, for later checks whether
+				   // the pivot element is large
+				   // enough, or whether we have to
+				   // fear that the matrix is not
+				   // regular
+  double diagonal_sum = 0;
+  for (unsigned int i=0; i<N; ++i)
+    diagonal_sum += std::fabs(el(i,i));
+  const double typical_diagonal_element = diagonal_sum/N;
+  
+  std::vector<unsigned int> p(N);
+  for (unsigned int i=0; i<N; ++i)
     p[i] = i;
 
-  for (unsigned int j=0; j<n(); ++j)
+  for (unsigned int j=0; j<N; ++j)
     {
-				       // pivotsearch
+				       // pivot search: search that
+				       // part of the line on and
+				       // right of the diagonal for
+				       // the largest element
       number       max = std::fabs(el(j,j));
       unsigned int r   = j;
-      for (unsigned int i=j+1; i<n(); ++i)
+      for (unsigned int i=j+1; i<N; ++i)
 	{
 	  if (std::fabs(el(i,j)) > max)
 	    {
@@ -1358,11 +1374,15 @@ FullMatrix<number>::gauss_jordan()
 	      r = i;
 	    }
 	}
-      Assert(max>1.e-16, ExcNotRegular());
-				       // rowinterchange
+				       // check whether the pivot is
+				       // too small
+      Assert(max > 1.e-16*typical_diagonal_element,
+	     ExcNotRegular(max));
+      
+				       // row interchange
       if (r>j)
 	{
-	  for (unsigned int k=0; k<n(); ++k)
+	  for (unsigned int k=0; k<N; ++k)
 	    std::swap (el(j,k), el(r,k));
 
 	  std::swap (p[j], p[r]);
@@ -1371,29 +1391,29 @@ FullMatrix<number>::gauss_jordan()
 				       // transformation
       const number hr = 1./el(j,j);
       el(j,j) = hr;
-      for (unsigned int k=0; k<n(); ++k)
+      for (unsigned int k=0; k<N; ++k)
 	{
 	  if (k==j) continue;
-	  for (unsigned int i=0; i<n(); ++i)
+	  for (unsigned int i=0; i<N; ++i)
 	    {
 	      if (i==j) continue;
 	      el(i,k) -= el(i,j)*el(j,k)*hr;
 	    }
 	}
-      for (unsigned int i=0; i<n(); ++i)
+      for (unsigned int i=0; i<N; ++i)
 	{
 	  el(i,j) *= hr;
 	  el(j,i) *= -hr;
 	}
       el(j,j) = hr;
     }
-				   // columninterchange
-  std::vector<number> hv(n());
-  for (unsigned int i=0; i<n(); ++i)
+				   // column interchange
+  std::vector<number> hv(N);
+  for (unsigned int i=0; i<N; ++i)
     {
-      for (unsigned int k=0; k<n(); ++k)
+      for (unsigned int k=0; k<N; ++k)
 	hv[p[k]] = el(i,k);
-      for (unsigned int k=0; k<n(); ++k)
+      for (unsigned int k=0; k<N; ++k)
 	el(i,k) = hv[k];
     }
 }
-- 
2.39.5