From bcac4fcfd429086f0e8b8cf8f4466c66b585e206 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 12 Feb 2006 04:41:10 +0000 Subject: [PATCH] Introduce a preconditioner. git-svn-id: https://svn.dealii.org/trunk@12327 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-20/step-20.cc | 82 +++++++++++++++++++++++------ 1 file changed, 66 insertions(+), 16 deletions(-) diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 90e7a479e6..384368b3cb 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -739,57 +739,63 @@ void MixedLaplaceProblem::assemble_system () } + +template class InverseMatrix : public Subscriptor { public: - InverseMatrix (const SparseMatrix &m); + InverseMatrix (const Matrix &m); void vmult (Vector &dst, const Vector &src) const; private: - const SmartPointer > matrix; + const SmartPointer matrix; mutable GrowingVectorMemory<> vector_memory; }; -InverseMatrix::InverseMatrix (const SparseMatrix &m) +template +InverseMatrix::InverseMatrix (const Matrix &m) : matrix (&m) {} -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const +template +void InverseMatrix::vmult (Vector &dst, + const Vector &src) const { - SolverControl solver_control (matrix->m(), 1e-8*src.l2_norm()); - SolverCG<> cg (solver_control, vector_memory); + SolverControl solver_control (src.size(), 1e-8*src.l2_norm()); + SolverCG<> cg (solver_control, vector_memory); + dst = 0; + cg.solve (*matrix, dst, src, PreconditionIdentity()); } -class SchurComplement +class SchurComplement : public Subscriptor { public: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix &Minv); + const InverseMatrix > &Minv); void vmult (Vector &dst, const Vector &src) const; private: const SmartPointer > system_matrix; - const SmartPointer m_inverse; + const SmartPointer > > m_inverse; mutable Vector tmp1, tmp2; }; SchurComplement::SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix &Minv) + const InverseMatrix > &Minv) : system_matrix (&A), m_inverse (&Minv), @@ -807,10 +813,47 @@ void SchurComplement::vmult (Vector &dst, } + +class ApproximateSchurComplement : public Subscriptor +{ + public: + ApproximateSchurComplement (const BlockSparseMatrix &A); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + + mutable Vector tmp1, tmp2; +}; + + +ApproximateSchurComplement::ApproximateSchurComplement (const BlockSparseMatrix &A) + : + system_matrix (&A), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void ApproximateSchurComplement::vmult (Vector &dst, + const Vector &src) const +{ + system_matrix->block(0,1).vmult (tmp1, src); + system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); +} + + + + + template void MixedLaplaceProblem::solve () { - const InverseMatrix m_inverse (system_matrix.block(0,0)); + const InverseMatrix > + m_inverse (system_matrix.block(0,0)); Vector tmp (solution.block(0).size()); { @@ -820,14 +863,21 @@ void MixedLaplaceProblem::solve () system_matrix.block(1,0).vmult (schur_rhs, tmp); schur_rhs -= system_rhs.block(1); + SchurComplement + schur_complement (system_matrix, m_inverse); + + ApproximateSchurComplement + approximate_schur_complement (system_matrix); + + InverseMatrix + preconditioner (approximate_schur_complement); + SolverControl solver_control (system_matrix.block(0,0).m(), 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - cg.solve (SchurComplement(system_matrix, m_inverse), - solution.block(1), - schur_rhs, - PreconditionIdentity()); + cg.solve (schur_complement, solution.block(1), schur_rhs, + preconditioner); std::cout << " " << solver_control.last_step() << " CG Schur complement iterations needed to obtain convergence." -- 2.39.5