]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge the new inversion code into the general template.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 21:29:51 +0000 (21:29 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 21:29:51 +0000 (21:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@17650 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/source/full_matrix.cc

index 08f678ce9a8f0a2b33dfafd8a12bf623fa5a38e9..643bc966bc597e96bbac8a82ce4cf3018c23e365 100644 (file)
@@ -20,6 +20,7 @@
 #include <base/template_constraints.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
+#include <lac/lapack_templates.h>
 
 #include <vector>
 #include <cmath>
@@ -1241,9 +1242,64 @@ FullMatrix<number>::gauss_jordan ()
 {
   Assert (!this->empty(), ExcEmptyMatrix());  
   Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
-  
-                                  // Gauss-Jordan-Algorithmus
-                                  // cf. Stoer I (4th Edition) p. 153
+
+                                  // see if we can use Lapack
+                                  // algorithms for this and if the
+                                  // type for 'number' works for us:
+#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && \
+    defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
+  if (types_are_equal<number,double>::value
+      ||
+      types_are_equal<number,float>::value)
+    {
+                                      // In case we have the LAPACK functions 
+                                      // getrf and getri detected at configure, 
+                                      // we use these algorithms for inversion 
+                                      // since they provide better performance 
+                                      // than the deal.II native functions. 
+                                      //
+                                      // Note that BLAS/LAPACK stores matrix 
+                                      // elements column-wise (i.e., all values in 
+                                      // one column, then all in the next, etc.), 
+                                      // whereas the FullMatrix stores them 
+                                      // row-wise.
+                                      // We ignore that difference, and give our
+                                      // row-wise data to LAPACK,
+                                      // let LAPACK build the inverse of the
+                                      // transpose matrix, and read the result as
+                                      // if it were row-wise again. In other words,
+                                      // we just got ((A^T)^{-1})^T, which is
+                                      // A^{-1}.
+
+      const int nn = this->n();
+      ipiv.resize(nn);
+      int info;
+
+                                      // Use the LAPACK function getrf for 
+                                      // calculating the LU factorization.
+      getrf(&nn, &nn, this->data(), &nn, &ipiv[0], &info);
+
+      Assert(info >= 0, ExcInternalError());
+      Assert(info == 0, LACExceptions::ExcSingular());
+
+      inv_work.resize (nn);
+                                      // Use the LAPACK function getri for
+                                      // calculating the actual inverse using
+                                      // the LU factorization.
+      getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
+
+      Assert(info >= 0, ExcInternalError());
+      Assert(info == 0, LACExceptions::ExcSingular());
+
+      return;
+    }
+
+#endif
+
+                                  // otherwise do it by hand. use the
+                                  // Gauss-Jordan-Algorithmus from
+                                  // Stoer & Bulirsch I (4th Edition)
+                                  // p. 153
   const unsigned int N = n();
 
                                   // first get an estimate of the
index bf7d74b8b4630d9e724bad16ef1c7e7e9f689ad0..0fa3339cdb1bc44d3e2a91da107bff519a3aee2d 100644 (file)
 
 
 #include <lac/full_matrix.templates.h>
-#include <lac/lapack_templates.h>
 #include <base/logstream.h>
 
 DEAL_II_NAMESPACE_OPEN
 
-
-                                   // Need to explicitly state the Lapack
-                                   // inversion since it only works with
-                                   // floats and doubles in case LAPACK was
-                                   // detected by configure.
-#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
-
-template <>
-void
-FullMatrix<float>::gauss_jordan ()
-{
-  Assert (!this->empty(), ExcEmptyMatrix());  
-  Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
-  
-                                   // In case we have the LAPACK functions 
-                                   // getrf and getri detected at configure, 
-                                   // we use these algorithms for inversion 
-                                   // since they provide better performance 
-                                   // than the deal.II native functions. 
-                                   //
-                                   // Note that BLAS/LAPACK stores matrix 
-                                   // elements column-wise (i.e., all values in 
-                                   // one column, then all in the next, etc.), 
-                                   // whereas the FullMatrix stores them 
-                                   // row-wise.
-                                   // We ignore that difference, and give our
-                                   // row-wise data to LAPACK,
-                                   // let LAPACK build the inverse of the
-                                   // transpose matrix, and read the result as
-                                   // if it were row-wise again. In other words,
-                                   // we just got ((A^T)^{-1})^T, which is
-                                   // A^{-1}.
-
-  const int nn = this->n();
-  float* values = const_cast<float*> (this->data());
-  ipiv.resize(nn);
-  int info;
-
-                                   // Use the LAPACK function getrf for 
-                                   // calculating the LU factorization.
-  getrf(&nn, &nn, values, &nn, &ipiv[0], &info);
-
-  Assert(info >= 0, ExcInternalError());
-  Assert(info == 0, LACExceptions::ExcSingular());
-
-  inv_work.resize (nn);
-                                   // Use the LAPACK function getri for
-                                   // calculating the actual inverse using
-                                   // the LU factorization.
-  getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
-
-  Assert(info >= 0, ExcInternalError());
-  Assert(info == 0, LACExceptions::ExcSingular());
-}
-
-template <>
-void
-FullMatrix<double>::gauss_jordan ()
-{
-  Assert (!this->empty(), ExcEmptyMatrix());  
-  Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
-  
-                                   // In case we have the LAPACK functions 
-                                   // getrf and getri detected at configure, 
-                                   // we use these algorithms for inversion 
-                                   // since they provide better performance 
-                                   // than the deal.II native functions. 
-                                   //
-                                   // Note that BLAS/LAPACK stores matrix 
-                                   // elements column-wise (i.e., all values in 
-                                   // one column, then all in the next, etc.), 
-                                   // whereas the FullMatrix stores them 
-                                   // row-wise.
-                                   // We ignore that difference, and give our
-                                   // row-wise data to LAPACK,
-                                   // let LAPACK build the inverse of the
-                                   // transpose matrix, and read the result as
-                                   // if it were row-wise again. In other words,
-                                   // we just got ((A^T)^{-1})^T, which is
-                                   // A^{-1}.
-
-  const int nn = this->n();
-  double* values = const_cast<double*> (this->data());
-  ipiv.resize(nn);
-  int info;
-
-                                   // Use the LAPACK function getrf for 
-                                   // calculating the LU factorization.
-  getrf(&nn, &nn, values, &nn, &ipiv[0], &info);
-
-  Assert(info >= 0, ExcInternalError());
-  Assert(info == 0, LACExceptions::ExcSingular());
-
-  inv_work.resize (nn);
-                                   // Use the LAPACK function getri for
-                                   // calculating the actual inverse using
-                                   // the LU factorization.
-  getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
-
-  Assert(info >= 0, ExcInternalError());
-  Assert(info == 0, LACExceptions::ExcSingular());
-}
-
-                                   // ... and now the usual instantiations
-                                   // of gauss_jordan() and all the rest.
-template void FullMatrix<long double>::gauss_jordan ();
-template void FullMatrix<std::complex<float> >::gauss_jordan ();
-template void FullMatrix<std::complex<double> >::gauss_jordan ();
-template void FullMatrix<std::complex<long double> >::gauss_jordan ();
-
-#else
-
-template void FullMatrix<float>::gauss_jordan ();
-template void FullMatrix<double>::gauss_jordan ();
-template void FullMatrix<long double>::gauss_jordan ();
-template void FullMatrix<std::complex<float> >::gauss_jordan ();
-template void FullMatrix<std::complex<double> >::gauss_jordan ();
-template void FullMatrix<std::complex<long double> >::gauss_jordan ();
-
-#endif
-
-
 #include "full_matrix.inst"
 
 
@@ -173,5 +50,4 @@ TEMPL_OP_EQ(std::complex<float>,std::complex<long double>);
 
 #undef TEMPL_OP_EQ
 
-
 DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.