From 76874a216004542243743a732045a0c98134d7a1 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 20 Feb 2016 18:41:57 -0500 Subject: [PATCH] Get rid of IterativeInverse in step-20. IterativeInverse was marked as deprecated in April 2015 in favor of LinearOperator. This puts this step in a strange place; we should not use deprecated features in tutorials, but at the same time we cannot use its replacement, which requires C++11 support. To resolve this, this step just provides its own replacement for a special case of IterativeInverse. --- examples/step-20/doc/intro.dox | 132 ++++++++++++-------------- examples/step-20/step-20.cc | 165 ++++++++++++++++----------------- 2 files changed, 139 insertions(+), 158 deletions(-) diff --git a/examples/step-20/doc/intro.dox b/examples/step-20/doc/intro.dox index 63ca74a0f7..9823e7097f 100644 --- a/examples/step-20/doc/intro.dox +++ b/examples/step-20/doc/intro.dox @@ -434,26 +434,25 @@ function vmult() that is all that the SolverCG class requires from an object representing a matrix. We can make our life a bit easier by also introducing an object that represents $M^{-1}$ and that has its own vmult() function that, if called, solves -the linear system with $M$; in fact, such a class already exists in -deal.II: this is accomplished by using the class -IterativeInverse. Using it, the class that implements the Schur -only needs to offer the vmult() +the linear system with $M$. Using this (which we will implement as the +InverseMatrix class in the body of the program), the class +that implements the Schur only needs to offer the vmult() function to perform a matrix-vector multiplication, using the algorithm above. Here are again the relevant parts of the code: @code -class SchurComplement +class SchurComplement : public Subscriptor { public: - SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv); + SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix > &inverse_mass); void vmult (Vector &dst, const Vector &src) const; private: const SmartPointer > system_matrix; - const SmartPointer > > m_inverse; + const SmartPointer > > inverse_mass; mutable Vector tmp1, tmp2; }; @@ -463,7 +462,7 @@ void SchurComplement::vmult (Vector &dst, const Vector &src) const { system_matrix->block(0,1).vmult (tmp1, src); // multiply with the top right block: B - m_inverse->vmult (tmp2, tmp1); // multiply with M^-1 + inverse_mass->vmult (tmp2, tmp1); // multiply with M^-1 system_matrix->block(1,0).vmult (dst, tmp2); // multiply with the bottom left block: B^T } @endcode @@ -490,37 +489,31 @@ matrix and the mass matrix, respectively: template void MixedLaplaceProblem::solve () { - PreconditionIdentity identity; - IterativeInverse > m_inverse; - m_inverse.initialize(system_matrix.block(0,0), identity); - m_inverse.solver.select("cg"); - static ReductionControl inner_control(1000, 0., 1.e-13); - m_inverse.solver.set_control(inner_control); - + InverseMatrix > inverse_mass (system_matrix.block(0,0)); Vector tmp (solution.block(0).size()); { + SchurComplement schur_complement (system_matrix, inverse_mass); Vector schur_rhs (solution.block(1).size()); - - m_inverse.vmult (tmp, system_rhs.block(0)); + inverse_mass.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); + SolverControl solver_control (solution.block(1).size(), + 1e-12*schur_rhs.l2_norm()); + SolverCG<> cg (solver_control); - cg.solve (SchurComplement(system_matrix, m_inverse), - solution.block(1), - schur_rhs, - PreconditionIdentity()); + PreconditionIdentity preconditioner; + cg.solve (schur_complement, solution.block(1), schur_rhs, + preconditioner); } + { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); - m_inverse.vmult (solution.block(0), tmp); + inverse_mass.vmult (solution.block(0), tmp); } } @endcode @@ -568,7 +561,7 @@ this looks almost as expensive as solving with $S$ right away. However, note that in the inner iteration, we do not have to calculate $M^{-1}$, but only the inverse of its diagonal, which is cheap. -The next step is to define a class that represents the approximate Schur +The next step is to define a class that represents this approximate Schur complement. This should look very much like the Schur complement class itself, except that it doesn't need the object representing $M^{-1}$ any more since we can compute the inverse of the diagonal of $M$ on the fly: @@ -576,26 +569,28 @@ since we can compute the inverse of the diagonal of $M$ on the fly: @code class ApproximateSchurComplement : public Subscriptor { - public: - ApproximateSchurComplement (const BlockSparseMatrix &A); +public: + ApproximateSchurComplement (const BlockSparseMatrix &A); - void vmult (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; - private: - const SmartPointer > system_matrix; +private: + const SmartPointer > system_matrix; - mutable Vector tmp1, tmp2; + mutable Vector tmp1, tmp2; }; -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); -} +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); + } @endcode Note how the vmult function differs in simply doing one Jacobi sweep @@ -606,18 +601,14 @@ the diagonal of $M$; in other words, the operation $({\textrm{diag}\ }M)^{-1}x$ on a vector $x$ is exactly what the function SparseMatrix::precondition_Jacobi above does.) -With all this, we already have the preconditioner: it should be the inverse of -the approximate Schur complement, i.e. we need code like this: +With all this, we nearly already have the preconditioner: it should be the +inverse of the approximate Schur complement. We implement this with +the InverseMatrix class: @code - ApproximateSchurComplement - approximate_schur_complement (system_matrix); - - IterativeInverse > - preconditioner; - preconditioner.initialize(approximate_schur_complement, identity); - preconditioner.solver.select("cg"); - preconditioner.solver.set_control(inner_control); + ApproximateSchurComplement approximate_schur (system_matrix); + InverseMatrix approximate_inverse + (approximate_schur); @endcode That's all! @@ -626,30 +617,21 @@ Taken together, the first block of our solve() function will then look like this: @code - Vector schur_rhs (solution.block(1).size()); - - m_inverse.vmult (tmp, system_rhs.block(0)); - 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); - - IterativeInverse > - preconditioner; - preconditioner.initialize(approximate_schur_complement, identity); - preconditioner.solver.select("cg"); - preconditioner.solver.set_control(inner_control); - - SolverControl solver_control (solution.block(1).size(), - 1e-12*schur_rhs.l2_norm()); - SolverCG<> cg (solver_control); - - cg.solve (schur_complement, solution.block(1), schur_rhs, - preconditioner); + SchurComplement schur_complement (system_matrix, inverse_mass); + Vector schur_rhs (solution.block(1).size()); + inverse_mass.vmult (tmp, system_rhs.block(0)); + system_matrix.block(1,0).vmult (schur_rhs, tmp); + schur_rhs -= system_rhs.block(1); + + SolverControl solver_control (solution.block(1).size(), + 1e-12*schur_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + ApproximateSchurComplement approximate_schur (system_matrix); + InverseMatrix approximate_inverse + (approximate_schur); + cg.solve (schur_complement, solution.block(1), schur_rhs, + approximate_inverse); @endcode Note how we pass the so-defined preconditioner to the solver working on the diff --git a/examples/step-20/step-20.cc b/examples/step-20/step-20.cc index 40260643a8..90fa8cc127 100644 --- a/examples/step-20/step-20.cc +++ b/examples/step-20/step-20.cc @@ -556,15 +556,62 @@ namespace Step20 // therefore not discuss the rationale for these classes here any more, but // rather only comment on implementational aspects. + // @sect4{The InverseMatrix class template} + + // There are a few places in this program where we will need either the + // action of the inverse of the mass matrix or the action of the inverse of + // the approximate Schur complement. Rather than explicitly calling + // SolverCG::solve every time that we need to solve such a system, we will + // wrap the action of either inverse in a simple class. The only things we + // would like to note are that this class is derived from + // Subscriptor and, as mentioned above, it stores a pointer to + // the underlying matrix with a SmartPointer object. + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix(const MatrixType &m); + + void vmult(Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + }; + + + template + InverseMatrix::InverseMatrix (const MatrixType &m) + : + matrix (&m) + {} + + + template + void InverseMatrix::vmult (Vector &dst, + const Vector &src) const + { + // To make the control flow simpler, we recreate both the ReductionControl + // and SolverCG objects every time this is called. This is not the most + // efficient choice because SolverCG instances allocate memory whenever + // they are created; this is just a tutorial so such inefficiencies are + // acceptable for the sake of exposition. + SolverControl solver_control (std::max(src.size(), static_cast (200)), + 1e-8*src.l2_norm()); + SolverCG<> cg (solver_control); + + dst = 0; + + cg.solve (*matrix, dst, src, PreconditionIdentity()); + } + // @sect4{The SchurComplement class template} // The next class is the Schur complement class. Its rationale has also been - // discussed in length in the introduction. The only things we would like to - // note is that the class, too, is derived from the Subscriptor - // class and that as mentioned above it stores pointers to the entire block - // matrix and the inverse of the mass matrix block using - // SmartPointer objects. + // discussed in length in the introduction. Like InverseMatrix, + // this class is derived from Subscriptor and stores SmartPointer s + // pointing to the system matrix and InverseMatrix wrapper. // // The vmult function requires two temporary vectors that we do // not want to re-allocate and free every time we call this function. Since @@ -583,25 +630,26 @@ namespace Step20 class SchurComplement : public Subscriptor { public: - SchurComplement (const BlockSparseMatrix &A, - const IterativeInverse > &Minv); + SchurComplement (const BlockSparseMatrix &A, + 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 IterativeInverse > &Minv) + SchurComplement + ::SchurComplement (const BlockSparseMatrix &A, + const InverseMatrix > &Minv) : system_matrix (&A), - m_inverse (&Minv), + inverse_mass (&Minv), tmp1 (A.block(0,0).m()), tmp2 (A.block(0,0).m()) {} @@ -619,29 +667,9 @@ namespace Step20 // @sect4{The ApproximateSchurComplement class template} // The third component of our solver and preconditioner system is the class - // that approximates the Schur complement so we can form a an InverseIterate - // object that approximates the inverse of the Schur complement. It follows - // the same pattern as the Schur complement class, with the only exception - // that we do not multiply with the inverse mass matrix in - // vmult, but rather just do a single Jacobi - // step. Consequently, the class also does not have to store a pointer to an - // inverse mass matrix object. - // - // We will later use this class as a template argument to the - // IterativeInverse class which will in turn want to use it as a - // template argument for the PointerMatrix class. The latter class - // has a function that requires us to also write a function that - // provides the product with the transpose of the matrix this object - // represents. As a consequence, in the code below, we also - // implement a Tvmult function here that represents the - // product of the transpose matrix with a vector. It is easy to see - // how this needs to be implemented here: since the matrix is - // symmetric, we can as well call vmult wherever the - // product with the transpose matrix is required. (Note, however, - // that even though we implement this function here, there will in - // fact not be any need for it as long as we use SolverCG as the - // solver since that solver does not ever call the function that - // provides this operation.) + // that approximates the Schur complement with the method described in the + // introduction. We will use this class to build a preconditioner for our + // system matrix. class ApproximateSchurComplement : public Subscriptor { public: @@ -649,8 +677,6 @@ namespace Step20 void vmult (Vector &dst, const Vector &src) const; - void Tvmult (Vector &dst, - const Vector &src) const; private: const SmartPointer > system_matrix; @@ -659,85 +685,58 @@ namespace Step20 }; - ApproximateSchurComplement::ApproximateSchurComplement (const BlockSparseMatrix &A) - : + + 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 + 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); } - - void ApproximateSchurComplement::Tvmult (Vector &dst, - const Vector &src) const - { - vmult (dst, src); - } - - - // @sect4{MixedLaplace::solve} // After all these preparations, we can finally write the function that // actually solves the linear problem. We will go through the two parts it // has that each solve one of the two equations, the first one for the // pressure (component 1 of the solution), then the velocities (component 0 - // of the solution). Both parts need an object representing the inverse mass - // matrix and an auxiliary vector, and we therefore declare these objects at - // the beginning of this function. + // of the solution). template void MixedLaplaceProblem::solve () { - ReductionControl inner_control(1000, 0., 1.e-13); - PreconditionIdentity identity; - IterativeInverse > m_inverse; - m_inverse.initialize(system_matrix.block(0,0), identity); - m_inverse.solver.select("cg"); - m_inverse.solver.set_control(inner_control); - + InverseMatrix > inverse_mass (system_matrix.block(0,0)); Vector tmp (solution.block(0).size()); // Now on to the first equation. The right hand side of it is $B^TM^{-1}F-G$, - // which is what we compute in the first few lines. We then declare the - // objects representing the Schur complement, its approximation, and the - // inverse of the approximation. Finally, we declare a solver object and - // hand off all these matrices and vectors to it to compute block 1 (the - // pressure) of the solution: + // which is what we compute in the first few lines: { + SchurComplement schur_complement (system_matrix, inverse_mass); Vector schur_rhs (solution.block(1).size()); - - m_inverse.vmult (tmp, system_rhs.block(0)); + inverse_mass.vmult (tmp, system_rhs.block(0)); 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); - - IterativeInverse > - preconditioner; - preconditioner.initialize(approximate_schur_complement, identity); - preconditioner.solver.select("cg"); - preconditioner.solver.set_control(inner_control); - - + // Now that we have the right hand side we can go ahead and solve for the + // pressure, using our approximation of the inverse as a preconditioner: SolverControl solver_control (solution.block(1).size(), 1e-12*schur_rhs.l2_norm()); - SolverCG<> cg (solver_control); + SolverCG<> cg (solver_control); + ApproximateSchurComplement approximate_schur (system_matrix); + InverseMatrix approximate_inverse + (approximate_schur); cg.solve (schur_complement, solution.block(1), schur_rhs, - preconditioner); + approximate_inverse); std::cout << solver_control.last_step() << " CG Schur complement iterations to obtain convergence." @@ -753,7 +752,7 @@ namespace Step20 tmp *= -1; tmp += system_rhs.block(0); - m_inverse.vmult (solution.block(0), tmp); + inverse_mass.vmult (solution.block(0), tmp); } } -- 2.39.5