]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added the functionality to calculate the inverse in LAPACKFullMatrix class.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Nov 2008 18:49:50 +0000 (18:49 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Nov 2008 18:49:50 +0000 (18:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@17547 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/configure
deal.II/configure.in
deal.II/doc/news/changes.h
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/lapack_full_matrix.h
deal.II/lac/include/lac/lapack_support.h
deal.II/lac/include/lac/lapack_templates.h
deal.II/lac/source/full_matrix.inst.in
deal.II/lac/source/lapack_full_matrix.cc

index 33b7a60eaa78e8bb42fd7354800977068c715065..8adaaed7337216e48e7ebb0afb9be54063f43db2 100755 (executable)
@@ -15420,7 +15420,7 @@ done
 
 
 
-for ac_func in dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_
+for ac_func in dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_
 do
 as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
 { echo "$as_me:$LINENO: checking for $ac_func" >&5
index fd9b75872932157e426dbb0894db0fb9e7f1a350..03371f18f178c1005d9a7db4e5da7b4a30f25ba3 100644 (file)
@@ -542,7 +542,7 @@ if test "x$with_blas" != "x" -a "x$with_blas" != "xno" ; then
 fi
 
 AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_])
-AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_])
+AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_])
 
 dnl -------------------------------------------------------------
 dnl       set include paths of several libraries
index 5586bff1fe2da06b642c9a6a738449afb8e39a9f..844155aa66324a024848ea564bc4f842a0d904ba 100644 (file)
@@ -322,6 +322,16 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li>
+  <p>
+  New: The class LAPACKFullMatrix can now invert full matrices using 
+  the (optimized) LAPACK functions getrf and getri. The speedup over 
+  the FullMatrix::gauss_jordan() function is a factor of two for matrices
+  with 100 rows and columns, and grows with matrix size.
+  <br>
+  (Martin Kronbichler 2008/11/11)
+  </p>
+
   <li>
   <p>
   Fixed: The BlockMatrixBase::clear() function that is used by all other
index 976398c612f0b318bfdf2f884dcdc4d358a4cede..dd89409b5cefac811dabc3bb7519dacd6f75715e 100644 (file)
@@ -19,6 +19,7 @@
 #include <base/table.h>
 #include <lac/exceptions.h>
 #include <lac/identity_matrix.h>
+#include <lac/lapack_full_matrix.h>
 
 #include <vector>
 #include <iomanip>
@@ -330,6 +331,16 @@ class FullMatrix : public Table<2,number>
     FullMatrix<number> &
     operator = (const IdentityMatrix &id);
     
+                                    /**
+                                     * Assignment operator for a
+                                     * LapackFullMatrix. The calling matrix
+                                     * must be of the same size as the
+                                     * LAPACK matrix.
+                                     */
+    template <typename number2>
+    FullMatrix<number> &
+    operator = (const LAPACKFullMatrix<number2>&);
+    
                                      /**
                                      * Assignment from different
                                      * matrix classes. This
index aafc65b128056389dcd147dd128bec0800da20e9..73af2bc4263f4a678dee3e89028a3e47c1b3e307 100644 (file)
@@ -108,6 +108,22 @@ FullMatrix<number>::operator = (const IdentityMatrix &id)
 
 
 
+template <typename number>
+template <typename number2>
+FullMatrix<number> &
+FullMatrix<number>::operator = (const LAPACKFullMatrix<number2>& M)
+{
+  Assert (this->m() == M.n_rows(), ExcDimensionMismatch(this->m(), M.n_rows()));
+  Assert (this->n() == M.n_cols(), ExcDimensionMismatch(this->n(), M.n_rows()));
+  for (unsigned int i=0;i<this->m();++i)
+    for (unsigned int j=0;j<this->n();++j)
+      (*this)(i,j) = M(i,j);
+  
+  return *this;
+}
+
+
+
 template <typename number>
 bool
 FullMatrix<number>::all_zero () const
index 0f238cef55fc92ac78e00c7d598eb5d751ee1147..4576f23625838fde6a636146c71b8e4f90aed27b 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -231,6 +231,14 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      */
     void compute_lu_factorization ();
 
+                                    /**
+                                     * Invert the matrix by first computing
+                                     * an LU factorization with the LAPACK
+                                     * function Xgetrf and then building
+                                     * the actual inverse using Xgetri.
+                                     */
+    void invert ();
+
                                     /**
                                      * Solve the linear system with
                                      * right hand side given by
@@ -370,6 +378,13 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      */
     std::vector<int> ipiv;
     
+                                    /**
+                                     * Workspace for calculating the
+                                     * inverse matrix from an LU
+                                     * factorization.
+                                     */
+    std::vector<number> iwork;
+    
                                     /**
                                      * Real parts of
                                      * eigenvalues. Filled by
index dcfab53d447b28db0e0f99e8cbdab967160ec6d1..034eb9f4e67d1ff5dee9b3bd61bb1fe76aca66f7 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -34,6 +34,8 @@ namespace LAPACKSupport
   {
                                         /// Contents is actually a matrix.
        matrix,
+                                        /// Contents is the inverse of a matrix.
+       inverse_matrix,
                                         /// Contents is an LU decomposition.
        lu,
                                         /// Eigenvalue vector is filled
@@ -51,6 +53,8 @@ namespace LAPACKSupport
       {
        case matrix:
              return "matrix";
+       case inverse_matrix:
+             return "inverse matrix";
        case lu:
              return "lu decomposition";
        case eigenvalues:
index cf1463678533fa0917f4b2e0ffe7d39795407720..a44da9470e5e1e55680f7878deddd42dea7efa84 100644 (file)
@@ -5,7 +5,7 @@
 //    This file was automatically generated from blas.h.in
 //    See blastemplates in the deal.II contrib directory
 //
-//    Copyright (C) 2005, 2006, 2007 by the deal authors
+//    Copyright (C) 2005, 2006, 2007, 2008 by the deal authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -49,6 +49,11 @@ void dgetrf_ (const int* m, const int* n, double* A,
              const int* lda, int* ipiv, int* info);
 void sgetrf_ (const int* m, const int* n, float* A,
              const int* lda, int* ipiv, int* info);
+// Invert matrix from LU factorization
+void dgetri_ (const int* n, double* A, const int* lda, 
+             int* ipiv, double* iwork, const int* lwork, int* info);
+void sgetri_ (const int* n, float* A, const int* lda, 
+             int* ipiv, float* iwork, const int* lwork, int* info);
 // Apply forward/backward substitution to LU factorization
 void dgetrs_ (const char* trans, const int* n, const int* nrhs,
              const double* A, const int* lda, const int* ipiv,
@@ -213,6 +218,36 @@ getrf (const int*, const int*, float*, const int*, int*, int*)
 #endif
 
 
+#ifdef HAVE_DGETRI_
+inline void
+getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info)
+{
+  dgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+}
+#else
+inline void
+getri (const int*, double*, const int*, int*, double*, const int*, int*)
+{
+  LAPACKSupport::ExcMissing("dgetrf");
+}
+#endif
+
+
+#ifdef HAVE_SGETRI_
+inline void
+getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info)
+{
+  sgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+}
+#else
+inline void
+getri (const int*, float*, const int*, int*, float*, const int*, int*)
+{
+  LAPACKSupport::ExcMissing("sgetrf");
+}
+#endif
+
+
 #ifdef HAVE_DGETRS_
 inline void
 getrs (const char* trans, const int* n, const int* nrhs, const double* A, const int* lda, const int* ipiv, double* b, const int* ldb, int* info)
index e1e0e89a839812720cae8b0532b340a9982fa07a..4b527ba37d29be702b28d3958b9af1d49832c2a7 100644 (file)
@@ -24,6 +24,9 @@ for (S : REAL_SCALARS)
 
 for (S1, S2 : REAL_SCALARS)
   {
+    template
+      FullMatrix<S1>& FullMatrix<S1>::operator = (const LAPACKFullMatrix<S2>&);
+
     template
       void FullMatrix<S1>::fill<S2> (
        const FullMatrix<S2>&, unsigned, unsigned, unsigned, unsigned);
index 4a8a01ff1adb20341910831b715f1c39b36cc596..20e8f44230a250d496ca8ec61c5088d8bca9de2a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -175,6 +175,38 @@ LAPACKFullMatrix<number>::compute_lu_factorization()
 }
 
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::invert()
+{
+  Assert(state == matrix || state == lu, 
+        ExcState(state));
+  const int mm = this->n_rows();
+  const int nn = this->n_cols();
+  Assert (nn == mm, ExcNotQuadratic());
+
+  number* values = const_cast<number*> (this->data());
+  ipiv.resize(mm);
+  int info;
+
+  if (state == matrix)
+    {
+      getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
+
+      Assert(info >= 0, ExcInternalError());
+      Assert(info == 0, LACExceptions::ExcSingular());
+    }
+
+  iwork.resize (mm);
+  getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info);
+
+  Assert(info >= 0, ExcInternalError());
+  Assert(info == 0, LACExceptions::ExcSingular());
+  
+  state = inverse_matrix;
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>& v,

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.