]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The implementation of the ILU decomposition was rather ineffecient. It has been repla...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 29 Feb 2008 23:03:09 +0000 (23:03 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 29 Feb 2008 23:03:09 +0000 (23:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@15821 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.templates.h

index 0fcbdb44b3051931149c46e7776ea6d4c33f8a23..210c49a359dadd087715cc0f68848753628b69bc 100644 (file)
@@ -214,6 +214,15 @@ an integer id of the current thread.
 <h3>lac</h3>
 
 <ol>
+<li> Fixed: The implementation of SparseILU::decompose was rather
+inefficient in that it accessed random elements of the matrix in its
+inner loop. It has been replaced by the algorithm given in the book
+by Yves Saad: "Iterative methods for sparse linear systems", second
+edition, in section 10.3.2.
+<br>
+(WB 2008/2/29)
+</li>
+
 <li> Fixed: The implementation of SparseILU::vmult very needlessly called
 SparseMatrix::vmult just to throw away the (nonsensical) result away
 immediately. This is now fixed.
index 77b85ecdfae5df9eb32d1c87e58b4eab50e3ea51..28b85652903c69c7d51e1facf22f93a901c94996 100644 (file)
@@ -31,21 +31,9 @@ DEAL_II_NAMESPACE_OPEN
  * in the decomposition that do not fit into the sparsity structure of this
  * object are discarded.
  *
- * The algorithm used by this class is as follows (indices run from @p 0
- * to @p N-1):
- * @verbatim
- * copy original matrix into a[i,j]
- * 
- * for i=1..N-1
- *   a[i-1,i-1] = a[i-1,i-1]^{-1}
- *
- *   for k=0..i-1
- *     a[i,k] = a[i,k] * a[k,k]
- *
- *     for j=k+1..N-1
- *      if (a[i,j] exists & a[k,j] exists)
- *        a[i,j] -= a[i,k] * a[k,j]
- * @endverbatim
+ * The algorithm used by this class is essentially a copy of the
+ * algorithm given in the book Y. Saad: "Iterative methods for sparse
+ * linear systems", second edition, in section 10.3.2.
  *
  * 
  * <h3>Usage and state management</h3>
@@ -58,8 +46,7 @@ DEAL_II_NAMESPACE_OPEN
  * @<double@></tt>; others can be generated in application programs (see the
  * section on @ref Instantiations in the manual).
  * 
- * @author Wolfgang Bangerth, 1999, based on a similar implementation
- * by Malte Braack; unified interface: Ralf Hartmann
+ * @author Wolfgang Bangerth, 2008; unified interface: Ralf Hartmann
  */
 template <typename number>
 class SparseILU : public SparseLUDecomposition<number>
index fb3d6722abe18bf2b185bd844306cd38be45fcbc..c7cd70031bd1876a85161c193b06fdc620b3ea7d 100644 (file)
@@ -20,6 +20,7 @@
 #include <algorithm>
 #include <cmath>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 template <typename number>
@@ -66,113 +67,99 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
   if (strengthen_diagonal>0)
     this->strengthen_diagonal_impl();
 
-  const SparsityPattern             &sparsity = this->get_sparsity_pattern();
-  const std::size_t  * const rowstart_indices = sparsity.get_rowstart_indices();
-  const unsigned int * const column_numbers   = sparsity.get_column_numbers();
-  
-/*
-  PSEUDO-ALGORITHM
-  (indices=0..N-1)
-  
-  for i=1..N-1
-    a[i-1,i-1] = a[i-1,i-1]^{-1}
-
-    for k=0..i-1
-      a[i,k] = a[i,k] * a[k,k]
-
-      for j=k+1..N-1
-        if (a[i,j] exists & a[k,j] exists)
-          a[i,j] -= a[i,k] * a[k,j]
-*/
+                                  // in the following, we implement
+                                  // algorithm 10.4 in the book by
+                                  // Saad by translating in essence
+                                  // the algorithm given at the end
+                                  // of section 10.3.2, using the
+                                  // names of variables used there
+  const SparsityPattern     &sparsity = this->get_sparsity_pattern();
+  const std::size_t  * const ia       = sparsity.get_rowstart_indices();
+  const unsigned int * const ja       = sparsity.get_column_numbers();
 
+  number * luval = &this->global_entry (0);
 
-                                  // i := row
-  for (unsigned int row=1; row<this->m(); ++row)
+  const unsigned int N = this->m();
+  
+  std::vector<unsigned int> iw (N, numbers::invalid_unsigned_int);
+  
+  for (unsigned int k=0; k<N; ++k)
     {
-                                      // invert diagonal element of the
-                                      // previous row. this is a hack,
-                                      // which is possible since this
-                                      // element is not needed any more
-                                      // in the process of decomposition
-                                      // and since it makes the backward
-                                      // step when applying the decomposition
-                                      // significantly faster
-      AssertThrow((this->global_entry(rowstart_indices[row-1]) != 0),
-                 ExcDivideByZero());
+      const unsigned int j1 = ia[k],
+                        j2 = ia[k+1]-1;
+      
+      for (unsigned int j=j1; j<=j2; ++j)
+       iw[ja[j]] = j;
+
+                                      // the algorithm in the book
+                                      // works on the elements of row
+                                      // k left of the
+                                      // diagonal. however, since we
+                                      // store the diagonal element
+                                      // at the first position, start
+                                      // at the element after the
+                                      // diagonal and run as long as
+                                      // we don't walk into the right
+                                      // half
+      unsigned int j = j1+1;
+
+      label_150:
       
-      this->global_entry (rowstart_indices[row-1])
-       = 1./this->global_entry (rowstart_indices[row-1]);
-
-                                      // let k run over all lower-left
-                                      // elements of row i; skip
-                                      // diagonal element at start
-      const unsigned int * first_of_row
-       = &column_numbers[rowstart_indices[row]+1];
-      const unsigned int * first_after_diagonal = this->prebuilt_lower_bound[row];
-
-                                      // k := *col_ptr
-      for (const unsigned int * col_ptr = first_of_row;
-           col_ptr!=first_after_diagonal; ++col_ptr)
-       {
-         const unsigned int global_index_ik = col_ptr-column_numbers;
-         this->global_entry(global_index_ik) *= this->diag_element(*col_ptr);
-
-                                          // now do the inner loop over
-                                          // j. note that we need to do
-                                          // it in the right order, i.e.
-                                          // taking into account that the
-                                          // columns are sorted within each
-                                          // row correctly, but excluding
-                                          // the main diagonal entry
-                                           //
-                                           // the explicit use of operator()
-                                           // works around a bug in some gcc
-                                           // versions (see PR 18803)
-         const int global_index_ki = sparsity.operator()(*col_ptr,row);
-
-         if (global_index_ki != -1)
-           this->diag_element(row) -= this->global_entry(global_index_ik) *
-                                      this->global_entry(global_index_ki);
-
-         for (const unsigned int * j = col_ptr+1;
-              j<&column_numbers[rowstart_indices[row+1]];
-              ++j)
-           {
-                                              // get the locations of
-                                              // entries ij and kj in
-                                              // the matrix. note that k<i, k<j
-             
-//TODO:[WB] make code faster by using the following comment          
-                                              // note: this inner loop could
-                                              // be made considerably faster
-                                              // if we consulted the row
-                                              // with number *col_ptr,
-                                              // instead of always asking
-                                              // sparsity(*col_ptr,*j),
-                                              // since we traverse this
-                                              // row linearly. I just didn't
-                                              // have the time to figure out
-                                              // the details.
-                                               //
-                                               // the explicit use of
-                                               // operator() works around a
-                                               // bug in some gcc versions
-                                               // (see PR 18803)
-                     const int global_index_ij = j - &column_numbers[0],
-                       global_index_kj = sparsity.operator()(*col_ptr,*j);
-             if ((global_index_ij != -1) &&
-                 (global_index_kj != -1))
-               this->global_entry(global_index_ij) -= this->global_entry(global_index_ik) *
-                                                      this->global_entry(global_index_kj);
-           };
-       };
-    };
-
-                                  // Here the very last diagonal
-                                  // element still has to be inverted
-                                  // because the for-loop doesn't do
-                                  // it...
-  this->diag_element(this->m()-1) = 1./this->diag_element(this->m()-1);
+      unsigned int jrow = ja[j];
+      if (jrow >= k)
+       goto label_200;
+
+                                      // actual computations:
+      {
+       number t1 = luval[j] * luval[ia[jrow]];
+       luval[j] = t1;
+
+                                        // jj runs from just right of
+                                        // the diagonal to the end of
+                                        // the row
+       unsigned int jj = ia[jrow]+1;
+       while (ja[jj] < jrow)
+         ++jj;
+       for (; jj<ia[jrow+1]; ++jj)
+         {
+           const unsigned int jw = iw[ja[jj]];
+           if (jw != numbers::invalid_unsigned_int)
+             luval[jw] -= t1 * luval[jj];
+         }
+
+       ++j;
+       if (j<=j2)
+         goto label_150;
+      }
+
+      label_200:
+
+                                      // in the book there is an
+                                      // assertion that we have hit
+                                      // the diagonal element,
+                                      // i.e. that jrow==k. however,
+                                      // we store the diagonal
+                                      // element at the front, so
+                                      // jrow must actually be larger
+                                      // than k or j is already in
+                                      // the next row
+      Assert ((jrow > k) || (j==ia[k+1]), ExcInternalError());
+
+                                      // now we have to deal with the
+                                      // diagonal element. in the
+                                      // book it is located at
+                                      // position 'j', but here we
+                                      // use the convention of
+                                      // storing the diagonal element
+                                      // first, so instead of j we
+                                      // use uptr[k]=ia[k]
+      Assert (luval[ia[k]] != 0, ExcInternalError());
+      
+      luval[ia[k]] = 1./luval[ia[k]];
+           
+      for (unsigned int j=j1; j<=j2; ++j)
+       iw[ja[j]] = numbers::invalid_unsigned_int;
+    }
 }
 
 

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.