]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Correctly implement gauss_jordan using LAPACK.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 22:04:18 +0000 (22:04 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Nov 2008 22:04:18 +0000 (22:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@17655 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 7a27ce59dc849ab4f8a5908016e5029a68d749a0..bf4aaa34df7d98878b7d84fd966fb29bbac27f6a 100644 (file)
@@ -1153,26 +1153,6 @@ class FullMatrix : public Table<2,number>
                                     //@}
     
     friend class Accessor;
-
-  private:
-
-#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
-                                     /**
-                                     * 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 643bc966bc597e96bbac8a82ce4cf3018c23e365..a6282b63f42830160b44259483422d09a88de2dc 100644 (file)
@@ -1235,6 +1235,43 @@ FullMatrix<number>::print_formatted (
 }
 
 
+namespace internal
+{
+                                  // a namespace into which we import
+                                  // the global definitions of the
+                                  // LAPACK functions getrf for
+                                  // various data types and add
+                                  // general templates for all other
+                                  // types. this allows to call the
+                                  // getrf function for all data
+                                  // types but generated an exception
+                                  // if something is called for a
+                                  // data type not supported by
+                                  // LAPACK.
+  namespace getrf_switch
+  {
+    using ::getrf;
+
+    template <typename T>
+    void
+    getrf (const int*, const int*, T*, const int*, int*, int*)
+    {
+      Assert (false, LAPACKSupport::ExcMissing("dgetr for this data type"));
+    }
+    
+
+    using ::getri;
+
+    template <typename T>
+    void
+    getri (const int*, T*, const int*, int*, T*, const int*, int*)
+    {
+      Assert (false, LAPACKSupport::ExcMissing("dgetri for this data type"));
+    }    
+  }
+}
+
+
 
 template <typename number>
 void
@@ -1272,21 +1309,25 @@ FullMatrix<number>::gauss_jordan ()
                                       // A^{-1}.
 
       const int nn = this->n();
-      ipiv.resize(nn);
+
+                                      // workspace for permutations
+      std::vector<int> ipiv(nn);
       int info;
 
                                       // Use the LAPACK function getrf for 
                                       // calculating the LU factorization.
-      getrf(&nn, &nn, this->data(), &nn, &ipiv[0], &info);
+      internal::getrf_switch::getrf(&nn, &nn, this->val, &nn, &ipiv[0], &info);
 
       Assert(info >= 0, ExcInternalError());
       Assert(info == 0, LACExceptions::ExcSingular());
 
-      inv_work.resize (nn);
+                                      // scratch array
+      std::vector<number> inv_work (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);
+      internal::getrf_switch::getri(&nn, this->val, &nn, &ipiv[0], &inv_work[0], &nn, &info);
 
       Assert(info >= 0, ExcInternalError());
       Assert(info == 0, LACExceptions::ExcSingular());

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.