]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
FullMatrix::invert inverts all square matrices
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Jun 2000 20:04:58 +0000 (20:04 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Jun 2000 20:04:58 +0000 (20:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@3106 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/precondition_block.templates.h

index b7e9af29980c040bf408f51af1c80f24b3130439..dd9813f903ae1be8dcc05f5b21ebc0ef267e443e 100644 (file)
@@ -365,11 +365,17 @@ class FullMatrix : public Subscriptor
     number norm2 () const;
     
                                     /**
-                                     * A=Inverse(A). Inversion of this by
-                                     * Gauss-Jordan-algorithm. Note that this
-                                     * is a rather expensive operation, so
-                                     * you may not want to use it for
-                                     * larger matrices if not necessary.
+                                     * A=Inverse(A). Inversion of
+                                     * this by Gauss-Jordan-algorithm
+                                     * with partial pivoting.  This
+                                     * process is well-behaved for
+                                     * positive definite matrices,
+                                     * but be aware of round-off
+                                     * errors in the indefinite case.
+                                     *
+                                     * The numerical effort to invert
+                                     * an @p{n x n} matrix is of the
+                                     * order @p{n³}.
                                      */
     void gauss_jordan ();
 
@@ -384,28 +390,16 @@ class FullMatrix : public Subscriptor
     double determinant () const;
 
                                     /**
-                                     * Assign the inverse of the given
-                                     * matrix to @p{*this}. This function is
-                                     * only implemented (hardcoded) for
-                                     * square matrices of dimension one,
-                                     * two, three and four, since the
-                                     * amount of code needed grows quickly.
-                                     * The implementation does not use
-                                     * an elimination method like the
-                                     * Gauss-Jordan one, but rather sets
-                                     * the element directly; their values
-                                     * are precomputed symbolically using
-                                     * Maple. This way, we can avoid the
-                                     * overhead of loops and local variables
-                                     * but the number of lines of code
-                                     * grows rapidly.
-                                     *
-                                     * For all other sizes than the ones given
-                                     * above, an exception of type
-                                     * @p{ExcNotImplemented(dim_range)} is
-                                     * thrown, which you can catch and use
-                                     * some other method to invert the matrix,
-                                     * e.g. the @p{gauss_jordan} function.
+                                     * Assign the inverse of the
+                                     * given matrix to
+                                     * @p{*this}. This function is
+                                     * hardcoded for square matrices
+                                     * of dimension one to four,
+                                     * since the amount of code
+                                     * needed grows quickly.  For
+                                     * larger matrices, the method
+                                     * @ref{gauss_jordan} is invoked
+                                     * implicitly.
                                      */
     void invert (const FullMatrix<number> &M);
 
index 9a2f13c4939299bedd4902fa35f9f37437d714ab..d423afa1c3c58de26de9a6c3b3746c5dc7c55019 100644 (file)
@@ -1276,7 +1276,8 @@ FullMatrix<number>::invert (const FullMatrix<number> &M)
       }
 
       default:
-           AssertThrow (false, ExcNotImplemented(dim_range));
+           *this = M;
+           gauss_jordan();
     };    
 };
 
@@ -1437,5 +1438,4 @@ FullMatrix<number>::least_squares(Vector<number2>& dst, Vector<number2>& src)
   return sqrt(sum);
 }
 
-
 #endif
index b71ed8f3e1d001a68f3e22108a657bfd356fd5a1..22a34040841151681285531c3c1af3d0a420b557 100644 (file)
@@ -15,6 +15,7 @@
 
 
 #include <base/exceptions.h>
+#include <base/logstream.h>
 #include <lac/precondition_block.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
@@ -112,6 +113,7 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
 
   if (same_diagonal)
     {
+      deallog << "PreconditionBlock uses only one diagonal block" << endl;
                                       // Invert only the first block
                                       // This is a copy of the code in the
                                       // 'else' part, stripped of the outer loop
@@ -121,15 +123,7 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
        for (unsigned int column_cell=0; column_cell<blocksize; ++column_cell)
          M_cell(row_cell,column_cell)=M.el(row_cell,column_cell);
 
-      if (blocksize <=4)
-       {  
-         _inverse[0].invert(M_cell);
-       }
-      else
-       {
-         M_cell.gauss_jordan();
-         _inverse[0]=M_cell;
-       }
+      _inverse[0].invert(M_cell);
     }
   else
     {
@@ -161,19 +155,8 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
            for (unsigned int column_cell=0, column=cell*blocksize;
                 column_cell<blocksize; ++column_cell, ++column)
              M_cell(row_cell,column_cell)=M.el(row,column);
-//      try
-//     {
-         if (blocksize <=4)
-           {  
-             _inverse[cell].invert(M_cell);
-           }
-         else
-           {
-             M_cell.gauss_jordan();
-             _inverse[cell]=M_cell;
-           }
-//      }
-//      catch (ExcNotImplemented &)
+
+         _inverse[cell].invert(M_cell);
        }
     }
 }

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.