Make SparseILU::vmult faster.

git-svn-id: https://svn.dealii.org/trunk@18622 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 851560f9a63dce97b441c529549219cb026234ae..f25798f2b3c62cab6a57b1586e7114b84fce5df2 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -593,6 +593,16 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li>
+  <p>
+  Updated: The SparseILU::vmult kernel has been re-written using the same
+  data structures as in SparseMatrix::vmult, which reduces the count of
+  operations by one third and the execution times on typical problems by
+  ten to twenty precent.
+  <br>
+  (Martin Kronbichler 2009/04/16)
+  </p>
+
   <li>
   <p>
   New: There is now a new class VectorView<Number> that allows views of
@@ -610,7 +620,7 @@ inconvenience this causes.
   <li>
   <p>
   Updated: The local_to_global functions in ConstraintMatrix got smarter,
-  which imcreases the speed of sparsity pattern generation with
+  which increases the speed of sparsity pattern generation with
   ConstraintMatrix argument, and makes writing into sparse matrices using
   distribute_local_to_global faster, in case there are many constraints.
   <br>

diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h
index 9050076badc575f914b8f6d3e67e770fec7b67e4..42b0b536fa585980f78143c4ee1f22b08e758b75 100644 (file)
--- a/deal.II/lac/include/lac/sparse_decomposition.templates.h
+++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h
@@ -178,6 +178,18 @@ template <typename somenumber>
 void
 SparseLUDecomposition<number>::copy_from (const SparseMatrix<somenumber>& matrix)
 {
+                                  // check whether we use the same sparsity
+                                  // pattern as the input matrix
+  if (&this->get_sparsity_pattern() == &matrix.get_sparsity_pattern())
+    {
+      const somenumber * input_ptr = matrix.val;
+      number * this_ptr = this->val;
+      const number * const end_ptr = this_ptr + this->n_nonzero_elements();
+      for ( ; this_ptr != end_ptr; ++input_ptr, ++this_ptr)
+       *this_ptr = *input_ptr;
+      return;
+    }
+
                                    // preset the elements
   std::fill_n (&this->global_entry(0),
                this->n_nonzero_elements(),

diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h
index c7cd70031bd1876a85161c193b06fdc620b3ea7d..c19807c119b36ebbe13ba1125b81993bad516ee3 100644 (file)
--- a/deal.II/lac/include/lac/sparse_ilu.templates.h
+++ b/deal.II/lac/include/lac/sparse_ilu.templates.h
@@ -54,7 +54,6 @@ template <typename somenumber>
 void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                   const double strengthen_diagonal)
 {
-  SparseLUDecomposition<number>::decompose (matrix, strengthen_diagonal);
   Assert (matrix.m()==matrix.n(), ExcNotQuadratic ());
   Assert (this->m()==this->n(),   ExcNotQuadratic ());
   Assert (matrix.m()==this->m(),  ExcDimensionMismatch(matrix.m(), this->m()));
@@ -62,7 +61,7 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
   Assert (strengthen_diagonal>=0,
          ExcInvalidStrengthening (strengthen_diagonal));
 
-  this->copy_from (matrix);
+  SparseLUDecomposition<number>::decompose (matrix, strengthen_diagonal);
 
   if (strengthen_diagonal>0)
     this->strengthen_diagonal_impl();
@@ -77,7 +76,7 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
   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);
+  number * luval = this->SparseMatrix<number>::val;
 
   const unsigned int N = this->m();
   
@@ -178,6 +177,7 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
     = this->get_sparsity_pattern().get_rowstart_indices();
   const unsigned int * const column_numbers
     = this->get_sparsity_pattern().get_column_numbers();
+
                                   // solve LUx=b in two steps:
                                   // first Ly = b, then
                                   //       Ux = y
@@ -199,9 +199,13 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
                                       // find the position where the part
                                       // right of the diagonal starts
       const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row];
-      
-      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
-       dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
+
+      somenumber dst_row = 0;
+      const number * luval = this->SparseMatrix<number>::val + 
+                            (rowstart - column_numbers);
+      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col, ++luval)
+       dst_row += *luval * dst(*col);
+      dst(row) -= dst_row;
     }
 
                                   // now the backward solve. same
@@ -220,8 +224,13 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
                                       // right of the diagonal starts
       const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row];
       
-      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
-       dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
+      somenumber dst_row = 0;
+      const number * luval = this->SparseMatrix<number>::val + 
+                            (first_after_diagonal - column_numbers);
+      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col, ++luval)
+       dst_row += *luval * dst(*col);
+
+      dst(row) -= dst_row;
 
                                       // scale by the diagonal element.
                                       // note that the diagonal element
@@ -244,6 +253,7 @@ void SparseILU<number>::Tvmult (Vector<somenumber>       &dst,
     = this->get_sparsity_pattern().get_rowstart_indices();
   const unsigned int * const column_numbers
     = this->get_sparsity_pattern().get_column_numbers();
+
                                   // solve (LU)'x=b in two steps:
                                   // first U'y = b, then
                                   //       L'x = y
@@ -268,9 +278,12 @@ void SparseILU<number>::Tvmult (Vector<somenumber>       &dst,
                                       // find the position where the part
                                       // right of the diagonal starts
       const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row];
-      
-      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
-       tmp(*col) += this->global_entry (col-column_numbers) * dst(row);
+
+      const somenumber dst_row = dst (row);
+      const number * luval = this->SparseMatrix<number>::val + 
+                            (first_after_diagonal - column_numbers);
+      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col, ++luval)
+       tmp(*col) += *luval * dst_row;
     }
 
                                   // now the backward solve. same
@@ -292,9 +305,12 @@ void SparseILU<number>::Tvmult (Vector<somenumber>       &dst,
                                       // find the position where the part
                                       // right of the diagonal starts
       const unsigned int * const first_after_diagonal = this->prebuilt_lower_bound[row];
-      
-      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
-       tmp(*col) += this->global_entry (col-column_numbers) * dst(row);
+
+      const somenumber dst_row = dst (row);
+      const number * luval = this->SparseMatrix<number>::val + 
+                            (rowstart - column_numbers);
+      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col, ++luval)
+       tmp(*col) += *luval * dst_row;
     }
 }