From 2c6d04580e81eee8bc7f2e70be4fab474df7d67b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 11 Feb 2006 22:22:07 +0000 Subject: [PATCH] Restructure solver somewhat. git-svn-id: https://svn.dealii.org/trunk@12322 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-20/step-20.cc | 101 +++++++++++++++------------- 1 file changed, 54 insertions(+), 47 deletions(-) diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index da7758fdc8..6ae8e9aab0 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -739,12 +739,50 @@ void MixedLaplaceProblem::assemble_system () } +class InverseMatrix : public Subscriptor +{ + public: + InverseMatrix (const SparseMatrix &matrix); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > mass_matrix; + + mutable GrowingVectorMemory<> vector_memory; +}; + + +InverseMatrix::InverseMatrix (const SparseMatrix &matrix) + : + mass_matrix (&matrix) +{} + + +void InverseMatrix::vmult (Vector &dst, + const Vector &src) const +{ + SolverControl solver_control (mass_matrix->m(), 1e-8*src.l2_norm()); + SolverCG<> cg (solver_control, vector_memory); + + cg.solve (*mass_matrix, dst, src, PreconditionIdentity()); + + std::cout << " " << solver_control.last_step() + << " inner iterations needed to obtain convergence." + << std::endl; +} + + + class SchurComplement { public: - SchurComplement (const BlockSparseMatrix &A) + SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix &Minv) : - A (A), + system_matrix (&A), + m_inverse (&Minv), tmp1 (A.block(0,0).m()), tmp2 (A.block(0,0).m()) {} @@ -752,26 +790,14 @@ class SchurComplement void vmult (Vector &dst, const Vector &src) const { - A.block(0,1).vmult (tmp1, src); - - SolverControl solver_control (tmp1.size(), - 1e-8*tmp1.l2_norm()); - SolverCG<> cg (solver_control, vector_memory); - - cg.solve (A.block(0,0), tmp2, tmp1, - PreconditionIdentity()); - - std::cout << " " << solver_control.last_step() - << " inner iterations needed to obtain convergence." - << std::endl; - - A.block(1,0).vmult (dst, tmp2); + system_matrix->block(0,1).vmult (tmp1, src); + m_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); } private: - const BlockSparseMatrix &A; - - mutable GrowingVectorMemory<> vector_memory; + const SmartPointer > system_matrix; + const SmartPointer m_inverse; mutable Vector tmp1, tmp2; }; @@ -781,31 +807,22 @@ class SchurComplement template void MixedLaplaceProblem::solve () { + InverseMatrix m_inverse (system_matrix.block(0,0)); + Vector tmp (solution.block(0).size()); + { Vector schur_rhs (solution.block(1).size()); - { - Vector tmp (solution.block(0).size()); - - SolverControl solver_control (system_matrix.block(0,0).m(), - 1e-6*system_rhs.l2_norm()); - SolverCG<> cg (solver_control); - cg.solve (system_matrix.block(0,0), tmp, - system_rhs.block(0), PreconditionIdentity()); - - std::cout << " " << solver_control.last_step() - << " CG mass matrix iterations needed to obtain convergence." - << std::endl; - - system_matrix.block(1,0).vmult (schur_rhs, tmp); - schur_rhs -= system_rhs.block(1); - } + m_inverse.vmult (tmp, system_rhs.block(0)); + system_matrix.block(1,0).vmult (schur_rhs, tmp); + schur_rhs -= system_rhs.block(1); SolverControl solver_control (system_matrix.block(0,0).m(), 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - cg.solve (SchurComplement(system_matrix), solution.block(1), + cg.solve (SchurComplement(system_matrix, m_inverse), + solution.block(1), schur_rhs, PreconditionIdentity()); @@ -814,21 +831,11 @@ void MixedLaplaceProblem::solve () << std::endl; } { - Vector tmp (solution.block(0).size()); system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); - SolverControl solver_control (system_matrix.block(0,0).m(), - 1e-6*tmp.l2_norm()); - SolverCG<> cg (solver_control); - - cg.solve (system_matrix.block(0,0), solution.block(0), - tmp, PreconditionIdentity()); - - std::cout << " " << solver_control.last_step() - << " CG mass matrix iterations needed to obtain convergence." - << std::endl; + m_inverse.vmult (solution.block(0), tmp); } } -- 2.39.5