From d479ea6a61d6b3683fdce9b16d6590175f71696e Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 20 Jul 2005 00:34:59 +0000
Subject: [PATCH] Duplicate code for inversion of matrices to avoid dependency
 to liblac

git-svn-id: https://svn.dealii.org/trunk@11180 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/base/source/symmetric_tensor.cc | 103 +++++++++++++++++++-----
 1 file changed, 85 insertions(+), 18 deletions(-)

diff --git a/deal.II/base/source/symmetric_tensor.cc b/deal.II/base/source/symmetric_tensor.cc
index 2ecfc291e8..47a686a485 100644
--- a/deal.II/base/source/symmetric_tensor.cc
+++ b/deal.II/base/source/symmetric_tensor.cc
@@ -13,44 +13,111 @@
 
 
 #include <base/symmetric_tensor.h>
-#include <lac/full_matrix.h>
+
 
 template <>
 SymmetricTensor<4,3>
 invert (const SymmetricTensor<4,3> &t)
 {
-  SymmetricTensor<4,3> tmp;
+  SymmetricTensor<4,3> tmp = t;
 
                                    // this function follows the exact same
                                    // scheme as the 2d case, except that
                                    // hardcoding the inverse of a 6x6 matrix
                                    // is pretty wasteful. instead, we use the
                                    // Gauss-Jordan algorithm implemented for
-                                   // FullMatrix
-  FullMatrix<double> m(6,6);
-  for (unsigned int i=0; i<6; ++i)
-    for (unsigned int j=0; j<6; ++j)
-      m(i,j) = t.data[i][j];
-  m.gauss_jordan ();
-
-                                   // copy back and scale rows and
-                                   // columns. the mult matrix here is diag[1,
-                                   // 1, 1, 1/2, 1/2, 1/2]
-  for (unsigned int i=0; i<3; ++i)
-    for (unsigned int j=0; j<3; ++j)
-      tmp.data[i][j] = m(i,j);
+                                   // FullMatrix; the following code is copied
+                                   // from there because using the FullMatrix
+                                   // class would introduce circular
+                                   // references between libbase and liblac
+  const unsigned int N = 6;
+
+				   // 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(tmp.data[i][i]);
+  const double typical_diagonal_element = diagonal_sum/N;
+  
+  unsigned int p[N];
+  for (unsigned int i=0; i<N; ++i)
+    p[i] = i;
+
+  for (unsigned int j=0; j<N; ++j)
+    {
+				       // pivot search: search that
+				       // part of the line on and
+				       // right of the diagonal for
+				       // the largest element
+      double       max = std::fabs(tmp.data[j][j]);
+      unsigned int r   = j;
+      for (unsigned int i=j+1; i<N; ++i)
+        if (std::fabs(tmp.data[i][j]) > max)
+          {
+            max = std::fabs(tmp.data[i][j]);
+            r = i;
+          }
+				       // check whether the pivot is
+				       // too small
+      Assert(max > 1.e-16*typical_diagonal_element,
+	     ExcMessage("This tensor seems to be noninvertible"));
+      
+				       // row interchange
+      if (r>j)
+	{
+	  for (unsigned int k=0; k<N; ++k)
+	    std::swap (tmp.data[j][k], tmp.data[r][k]);
+
+	  std::swap (p[j], p[r]);
+	}
+
+				       // transformation
+      const double hr = 1./tmp.data[j][j];
+      tmp.data[j][j] = hr;
+      for (unsigned int k=0; k<N; ++k)
+	{
+	  if (k==j) continue;
+	  for (unsigned int i=0; i<N; ++i)
+	    {
+	      if (i==j) continue;
+	      tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
+	    }
+	}
+      for (unsigned int i=0; i<N; ++i)
+	{
+	  tmp.data[i][j] *= hr;
+	  tmp.data[j][i] *= -hr;
+	}
+      tmp.data[j][j] = hr;
+    }
+				   // column interchange
+  double hv[N];
+  for (unsigned int i=0; i<N; ++i)
+    {
+      for (unsigned int k=0; k<N; ++k)
+	hv[p[k]] = tmp.data[i][k];
+      for (unsigned int k=0; k<N; ++k)
+	tmp.data[i][k] = hv[k];
+    }
 
+                                   // scale rows and columns. the mult matrix
+                                   // here is diag[1, 1, 1, 1/2, 1/2, 1/2]
   for (unsigned int i=3; i<6; ++i)
     for (unsigned int j=0; j<3; ++j)
-      tmp.data[i][j] = m(i,j) / 2;
+      tmp.data[i][j] /= 2;
 
   for (unsigned int i=0; i<3; ++i)
     for (unsigned int j=3; j<6; ++j)
-      tmp.data[i][j] = m(i,j) / 2;
+      tmp.data[i][j] /= 2;
 
   for (unsigned int i=3; i<6; ++i)
     for (unsigned int j=3; j<6; ++j)
-      tmp.data[i][j] = m(i,j) / 4;
+      tmp.data[i][j] /= 4;
   
   return tmp;
 }
-- 
2.39.5