From 9dc82af995038720efcc2a83135406f85f0d5d9c Mon Sep 17 00:00:00 2001 From: oliver Date: Thu, 15 Sep 2005 10:09:08 +0000 Subject: [PATCH] Implemented the Tvmult method for the SparseILU preconditioner. git-svn-id: https://svn.dealii.org/trunk@11443 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_ilu.h | 15 ++++ .../lac/include/lac/sparse_ilu.templates.h | 73 ++++++++++++++++++- deal.II/lac/source/sparse_ilu.cc | 8 ++ 3 files changed, 94 insertions(+), 2 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index e8175cc612..6743edabe3 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -154,6 +154,21 @@ class SparseILU : public SparseLUDecomposition void vmult (Vector &dst, const Vector &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 + void Tvmult (Vector &dst, + const Vector &src) const; + + /** * Determine an estimate for the * memory consumption (in bytes) diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 27b5731cdd..67e1cd3325 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -220,8 +220,8 @@ void SparseILU::vmult (Vector &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::vmult (Vector &dst, } +template +template +void SparseILU::Tvmult (Vector &dst, + const Vector &src) const +{ + SparseLUDecomposition::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 tmp (N); + + dst = src; + for (unsigned int row=0; rowdiag_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 unsigned int diff --git a/deal.II/lac/source/sparse_ilu.cc b/deal.II/lac/source/sparse_ilu.cc index abe906ede8..24ecc3d86c 100644 --- a/deal.II/lac/source/sparse_ilu.cc +++ b/deal.II/lac/source/sparse_ilu.cc @@ -20,12 +20,16 @@ template void SparseILU::decompose (const SparseMatrix & const double); template void SparseILU::vmult (Vector &, const Vector &) const; +template void SparseILU::Tvmult (Vector &, + const Vector &) const; template void SparseILU::initialize (const SparseMatrix &, const AdditionalData data); template void SparseILU::decompose (const SparseMatrix &, const double); template void SparseILU::vmult (Vector &, const Vector &) const; +template void SparseILU::Tvmult (Vector &, + const Vector &) const; template class SparseILU; @@ -35,9 +39,13 @@ template void SparseILU::decompose (const SparseMatrix &, const double); template void SparseILU::vmult (Vector &, const Vector &) const; +template void SparseILU::Tvmult (Vector &, + const Vector &) const; template void SparseILU::initialize (const SparseMatrix &, const AdditionalData data); template void SparseILU::decompose (const SparseMatrix &, const double); template void SparseILU::vmult (Vector &, const Vector &) const; +template void SparseILU::Tvmult (Vector &, + const Vector &) const; -- 2.39.5