From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 8 Feb 2001 14:39:49 +0000 (+0000)
Subject: Fix SparseILU by a patch by Oliver Rheinbach <or@winfos.com>. Add a new test to check... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73027ad0e1d977f9791f141533fb2929d3e2414a;p=dealii-svn.git

Fix SparseILU by a patch by Oliver Rheinbach <or@winfos.com>. Add a new test to check for exact inversion in case we allow the ILU to use all elements of the matrix. Also check that the residual is small in case we use a rather sparse pattern for the ILU (the same one in fact as for the original five-point matrix.


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 f21b2852eb..1e48070b64 100644
--- 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 15f1d3990e..b568c769da 100644
--- 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 753f3b4093..c77bd18ac5 100644
--- 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)-