Fix SparseILU by a patch by Oliver Rheinbach <or@winfos.com>. Add a new test to check...

git-svn-id: https://svn.dealii.org/branches/Branch-3-1@3873 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html b/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html
index f21b2852ebf726681fa0eef8d0613d7b743065d5..1e48070b648cf2ec74ba698fb7f66cb2d0b4dfc5 100644 (file)
--- a/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html
+++ b/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html
@@ -279,6 +279,12 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       Fixed: <code class="class">SparseILU</code> had two bugs which
+       are now fixed.
+       <br>
+       (<a href="mailto:or@winfos.com">Oliver Rheinbach</a>, 2001/02/02)
+       </p>
 
   <li> <p>
        Improved: block classes have a variable number of blocks now,

diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h
index 15f1d3990e0768fccb2a547916415dc92821c4d5..b568c769da0a8ef5acff3dca1256d25282f1766d 100644 (file)
--- a/deal.II/lac/include/lac/sparse_ilu.h
+++ b/deal.II/lac/include/lac/sparse_ilu.h
@@ -164,16 +164,20 @@ class SparseILU : protected SparseMatrix<number>
     void reinit (const SparsityPattern &sparsity);
 
                                     /**
-                                     * Perform the incomplete LU factorization
-                                     * of the given matrix.
+                                     * Perform the incomplete LU
+                                     * factorization of the given
+                                     * matrix.
                                      *
-                                     * Note that the sparsity structures of
-                                     * the decomposition and the matrix passed
-                                     * to this function need not be equal,
-                                     * but that the pattern used by this
-                                     * matrix needs to contain all elements
-                                     * used by the matrix to be decomposed.
-                                     * Fill-in is thus allowed.
+                                     * Note that the sparsity
+                                     * structures of the
+                                     * decomposition and the matrix
+                                     * passed to this function need
+                                     * not be equal, but that the
+                                     * pattern used by this matrix
+                                     * needs to contain all elements
+                                     * used by the matrix to be
+                                     * decomposed.  Fill-in is thus
+                                     * allowed.
                                      */
     template <typename somenumber>
     void decompose (const SparseMatrix<somenumber> &matrix,
@@ -189,7 +193,9 @@ class SparseILU : protected SparseMatrix<number>
                              const Vector<somenumber> &src) const;
 
                                     /**
-                                     * Same as @p{apply_decomposition}, format for LAC.
+                                     * Same as
+                                     * @p{apply_decomposition},
+                                     * format for LAC.
                                      */
     template <typename somenumber>
     void vmult (Vector<somenumber>       &dst,

diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h
index 753f3b40934aca9c0b233ec906182a90574bf680..c77bd18ac5751d617591d79d6d5f79162396a565 100644 (file)
--- a/deal.II/lac/include/lac/sparse_ilu.templates.h
+++ b/deal.II/lac/include/lac/sparse_ilu.templates.h
@@ -176,7 +176,12 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                           // columns are sorted within each
                                           // row correctly, but excluding
                                           // the main diagonal entry
-         bool left_of_diagonal = true;
+         const int global_index_ki = sparsity(*col_ptr,row);
+
+         if (global_index_ki != -1)
+           diag_element(row) -= global_entry(global_index_ik) *
+                                global_entry(global_index_ki);
+
          for (const unsigned int * j = col_ptr+1;
               j<&column_numbers[rowstart_indices[row+1]];
               ++j)
@@ -191,26 +196,6 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                               // row linearly. I just didn't
                                               // have the time to figure out
                                               // the details.
-             
-                                              // check whether we have just
-                                              // traversed the diagonal
-                                              //
-                                              // note that j should never point
-                                              // to the diagonal itself!
-             if (left_of_diagonal && (*j > row))
-               {
-                 Assert (*j != row, ExcInternalError());
-                 
-                 left_of_diagonal = false;
-                                                  // a[i,i] -= a[i,k]*a[k,i]
-                 const int global_index_ki = sparsity(*col_ptr,row);
-
-                 if (global_index_ki != -1)
-                   diag_element(row) -= global_entry(global_index_ik) *
-                                        global_entry(global_index_ki);
-                 
-               };
-             
                      const int global_index_ij = j - &column_numbers[0],
                        global_index_kj = sparsity(*col_ptr,*j);
              if ((global_index_ij != -1) &&
@@ -221,6 +206,12 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
        };
     };
 
+                                  // Here the very last diagonal
+                                  // element still has to be inverted
+                                  // because the for-loop doesn't do
+                                  // it...
+ diag_element(m()-1)=1./diag_element(m()-1);
+
 /*
   OLD CODE, rather crude first implementation with an algorithm taken
   from 'W. Hackbusch, G. Wittum: Incomplete Decompositions (ILU)-