From: Wolfgang Bangerth Date: Wed, 20 Jul 2005 00:34:59 +0000 (+0000) Subject: Duplicate code for inversion of matrices to avoid dependency to liblac X-Git-Tag: v8.0.0~13386 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d8dc00be8dfe471ead1562165eed36ab4b27d315;p=dealii.git 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 --- 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 -#include + 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 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 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