]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Invert full matrices using Lapack functions if these are detected. Fixed a constructo...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 13:25:13 +0000 (13:25 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 13:25:13 +0000 (13:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@17641 0785d39b-7218-0410-832d-ea1e28bc413d

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_templates.h
deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/source/full_matrix.cc
deal.II/lac/source/lapack_full_matrix.cc
deal.II/lac/source/trilinos_block_vector.cc

index dd89409b5cefac811dabc3bb7519dacd6f75715e..ffbebca87ca238d6e044cd58135cbaf2a16878b2 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -906,7 +906,15 @@ class FullMatrix : public Table<2,number>
                                      * well-behaved for positive
                                      * definite matrices, but be
                                      * aware of round-off errors in
-                                     * the indefinite case.
+                                     * the indefinite case. 
+                                     *
+                                     * In case deal.II was configured with
+                                     * LAPACK, the functions Xgetrf and
+                                     * Xgetri build an LU factorization and
+                                     * invert the matrix upon that
+                                     * factorization, providing best
+                                     * performance up to matrices with a
+                                     * few hundreds rows and columns.
                                      *
                                      * The numerical effort to invert
                                      * an <tt>n x n</tt> matrix is of the
@@ -1145,6 +1153,26 @@ class FullMatrix : public Table<2,number>
                                     //@}
     
     friend class Accessor;
+
+  private:
+
+#ifdef HAVE_LIBLAPACK
+                                     /**
+                                     * The vector storing the
+                                     * permutations applied for
+                                     * pivoting in the
+                                     * LU-factorization.
+                                     */
+    std::vector<int> ipiv;
+    
+                                    /**
+                                     * Workspace for calculating the
+                                     * inverse matrix from an LU
+                                     * factorization.
+                                     */
+    std::vector<number> inv_work;
+#endif
+
 };
 
 /**@}*/
index 73af2bc4263f4a678dee3e89028a3e47c1b3e307..08f678ce9a8f0a2b33dfafd8a12bf623fa5a38e9 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index 4576f23625838fde6a636146c71b8e4f90aed27b..e84ba90d6ee2c62c78dec7a93b3429481998f146 100644 (file)
@@ -383,7 +383,7 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      * inverse matrix from an LU
                                      * factorization.
                                      */
-    std::vector<number> iwork;
+    std::vector<number> inv_work;
     
                                     /**
                                      * Real parts of
index a44da9470e5e1e55680f7878deddd42dea7efa84..89fcebf77fc62fdbc89dc1d32344dfafa3c6201a 100644 (file)
@@ -51,9 +51,9 @@ 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);
+             int* ipiv, double* inv_work, const int* lwork, int* info);
 void sgetri_ (const int* n, float* A, const int* lda, 
-             int* ipiv, float* iwork, const int* lwork, int* info);
+             int* ipiv, float* inv_work, 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,
@@ -218,32 +218,32 @@ getrf (const int*, const int*, float*, const int*, int*, int*)
 #endif
 
 
-#ifdef HAVE_DGETRI_
+#ifdef HAVE_DGETRF_
 inline void
-getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info)
+getri (const int* n, double* A, const int* lda, int* ipiv, double* inv_work, const int* lwork, int* info)
 {
-  dgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+  dgetri_ (n,A,lda,ipiv,inv_work,lwork,info);
 }
 #else
 inline void
 getri (const int*, double*, const int*, int*, double*, const int*, int*)
 {
-  LAPACKSupport::ExcMissing("dgetrf");
+  LAPACKSupport::ExcMissing("dgetri");
 }
 #endif
 
 
-#ifdef HAVE_SGETRI_
+#ifdef HAVE_SGETRF_
 inline void
-getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info)
+getri (const int* n, float* A, const int* lda, int* ipiv, float* inv_work, const int* lwork, int* info)
 {
-  sgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+  sgetri_ (n,A,lda,ipiv,inv_work,lwork,info);
 }
 #else
 inline void
 getri (const int*, float*, const int*, int*, float*, const int*, int*)
 {
-  LAPACKSupport::ExcMissing("sgetrf");
+  LAPACKSupport::ExcMissing("sgetri");
 }
 #endif
 
index 34534dbb59ef9ce7a509625252065a4645c27168..ce144e0079453ff4dc72d9c1a3e8f7b946918b9a 100644 (file)
@@ -836,7 +836,11 @@ namespace TrilinosWrappers
                     :
                     BlockVectorBase<Vector > ()
   {
-    reinit(v);
+    this->components.resize (v.n_blocks());
+    this->block_indices = v.block_indices;
+    
+    for (unsigned int i=0; i<this->n_blocks(); ++i)
+      this->components[i] = v.components[i];
   }
 
 
index 0fa3339cdb1bc44d3e2a91da107bff519a3aee2d..e5c70d1090b5f991d9a69ff6cda667cd5e9744f2 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.
+#ifdef HAVE_LIBLAPACK
+
+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.
+  std::cout << " blas";
+  getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
+  std::cout << " 2 mal" << std::endl;
+
+  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"
 
 
@@ -50,4 +175,5 @@ TEMPL_OP_EQ(std::complex<float>,std::complex<long double>);
 
 #undef TEMPL_OP_EQ
 
+
 DEAL_II_NAMESPACE_CLOSE
index 20e8f44230a250d496ca8ec61c5088d8bca9de2a..386b85703feba12c68d1a32b52f2c79ccd1425e2 100644 (file)
@@ -197,8 +197,8 @@ LAPACKFullMatrix<number>::invert()
       Assert(info == 0, LACExceptions::ExcSingular());
     }
 
-  iwork.resize (mm);
-  getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info);
+  inv_work.resize (mm);
+  getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info);
 
   Assert(info >= 0, ExcInternalError());
   Assert(info == 0, LACExceptions::ExcSingular());
index 0d64f039474ce943958e86de32a0ee433e0bccb5..369b99ecc941586981d1428d01ee363de861d1d4 100644 (file)
@@ -231,9 +231,7 @@ namespace TrilinosWrappers
       components.resize(n_blocks());
   
     for (unsigned int i=0;i<n_blocks();++i)
-      components[i].reinit(v.block(i), true, false);
-      
-    collect_sizes();
+      components[i] = v.block(i);
   }
 
 

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.