]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implemented the Tvmult method for the SparseILU preconditioner.
authoroliver <oliver@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Sep 2005 10:09:08 +0000 (10:09 +0000)
committeroliver <oliver@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Sep 2005 10:09:08 +0000 (10:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@11443 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/source/sparse_ilu.cc

index e8175cc61294314f82b508023454b1ac4fdcd280..6743edabe3f6107e9d5dedf14f4d67fd80fda5e5 100644 (file)
@@ -154,6 +154,21 @@ class SparseILU : public SparseLUDecomposition<number>
     void vmult (Vector<somenumber>       &dst,
                const Vector<somenumber> &src) const;
 
+
+                                    /**
+                                     * Apply the transpose of the 
+                                     * incomplete decomposition,
+                                     * i.e. do one forward-backward step
+                                     * $dst=(LU)^{-T}src$.
+                                     *
+                                     * The @p initialize function
+                                     * needs to be called beforehand.
+                                     */
+    template <typename somenumber>
+    void Tvmult (Vector<somenumber>       &dst,
+                const Vector<somenumber> &src) const;
+
+
                                     /**
                                      * Determine an estimate for the
                                      * memory consumption (in bytes)
index 27b5731cdd3be648cedc49e016b2452cf078e298..67e1cd3325d9226445fd237b5c08f50154079e53 100644 (file)
@@ -220,8 +220,8 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
                                   // done.
                                   //
                                   // note that we need to scale now,
-                                  // since the diagonal is not zero
-                                  // now
+                                  // since the diagonal is not equal to
+                                   // one now
   for (int row=N-1; row>=0; --row)
     {
                                       // get end of this row
@@ -241,6 +241,75 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
 }
 
 
+template <typename number>
+template <typename somenumber>
+void SparseILU<number>::Tvmult (Vector<somenumber>       &dst,
+                               const Vector<somenumber> &src) const 
+{
+  SparseLUDecomposition<number>::Tvmult (dst, src);
+
+  Assert (dst.size() == src.size(), ExcDimensionMismatch(dst.size(), src.size()));
+  Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
+  
+  const unsigned int N=dst.size();
+  const unsigned int * const rowstart_indices
+    = 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
+                                  //
+                                  // first a forward solve. Due to the
+                                   // fact that the transpose of U'
+                                   // is not easily accessible, a
+                                   // temporary vector is required.
+  Vector<somenumber> tmp (N);
+
+  dst = src;
+  for (unsigned int row=0; row<N; ++row)
+    {
+      dst(row) -= tmp (row);
+                                      // scale by the diagonal element.
+                                      // note that the diagonal element
+                                      // was stored inverted
+      dst(row) *= this->diag_element(row);
+
+                                      // get end of this row
+      const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
+                                      // 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);
+    };
+
+                                  // now the backward solve. same
+                                  // procedure, but we need not set
+                                  // dst before, since this is already
+                                  // done.
+                                  //
+                                  // note that we no scaling is required
+                                   // now, since the diagonal is one
+                                  // now
+  tmp = 0;
+  for (int row=N-1; row>=0; --row)
+    {
+      dst(row) -= tmp (row);
+
+                                      // get start of this row. skip the
+                                      // diagonal element
+      const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
+                                      // 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);
+    };
+}
+
 
 template <typename number>
 unsigned int
index abe906ede84a63dc4eb84e682bb48a8ce63c03ca..24ecc3d86cbf9ef8742cd1226422138e7d9258c1 100644 (file)
@@ -20,12 +20,16 @@ template void SparseILU<double>::decompose<double> (const SparseMatrix<double> &
                                                    const double);
 template void SparseILU<double>::vmult <double> (Vector<double> &,
                                                  const Vector<double> &) const;
+template void SparseILU<double>::Tvmult <double> (Vector<double> &,
+                                                 const Vector<double> &) const;
 template void SparseILU<double>::initialize<float> (const SparseMatrix<float> &,
                                                    const AdditionalData data);
 template void SparseILU<double>::decompose<float> (const SparseMatrix<float> &,
                                                   const double);
 template void SparseILU<double>::vmult<float> (Vector<float> &,
                                                const Vector<float> &) const;
+template void SparseILU<double>::Tvmult<float> (Vector<float> &,
+                                               const Vector<float> &) const;
 
 
 template class SparseILU<float>;
@@ -35,9 +39,13 @@ template void SparseILU<float>::decompose<double> (const SparseMatrix<double> &,
                                                   const double);
 template void SparseILU<float>::vmult<double> (Vector<double> &,
                                                const Vector<double> &) const;
+template void SparseILU<float>::Tvmult<double> (Vector<double> &,
+                                               const Vector<double> &) const;
 template void SparseILU<float>::initialize<float> (const SparseMatrix<float> &,
                                                   const AdditionalData data);
 template void SparseILU<float>::decompose<float> (const SparseMatrix<float> &,
                                                  const double);
 template void SparseILU<float>::vmult<float> (Vector<float> &,
                                               const Vector<float> &) const;
+template void SparseILU<float>::Tvmult<float> (Vector<float> &,
+                                              const Vector<float> &) const;

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.