From 3bdd80eadce586c8ceeaf250c3f63d0ec7bc9fde Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 12 Feb 2006 19:15:14 +0000 Subject: [PATCH] More documentation. git-svn-id: https://svn.dealii.org/trunk@12334 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-20/step-20.cc | 150 +++++++++++++++++++++++++++- 1 file changed, 147 insertions(+), 3 deletions(-) diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 384368b3cb..8821bab96e 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -739,7 +739,138 @@ void MixedLaplaceProblem::assemble_system () } - + // @sect3{Linear solvers and preconditioners} + + // The linear solvers and + // preconditioners we use in this + // example have been discussed in + // significant detail already in the + // introduction. We will therefore + // not discuss the rationale for + // these classes here any more, but + // rather only comment on + // implementational aspects. + + // @sect4{The ``InverseMatrix'' class template} + + // The first component of our linear + // solver scheme was the creation of + // a class that acts like the inverse + // of a matrix, i.e. which has a + // ``vmult'' function that multiplies + // a vector with an inverse matrix by + // solving a linear system. + // + // While most of the code below + // should be obvious given the + // purpose of this class, two + // comments are in order. First, the + // class is derived from the + // ``Subscriptor'' class so that we + // can use the ``SmartPointer'' class + // with inverse matrix objects. The + // use of the ``Subscriptor'' class + // has been explained before in + // step-7 and step-20. The present + // class also sits on the receiving + // end of this + // ``Subscriptor''/``SmartPointer'' + // pair: it holds its pointer to the + // matrix it is supposed to be the + // inverse of through a + // ``SmartPointer'' to make sure that + // this matrix is not destroyed while + // we still have a pointer to it. + // + // Secondly, we realize that we will + // probably perform many + // matrix-vector products with + // inverse matrix objects. Now, every + // time we do so, we have to call the + // CG solver to solve a linear + // system. To work, the CG solver + // needs to allocate four temporary + // vectors that it will release again + // at the end of its operation. What + // this means is that through + // repeated calls to the ``vmult'' + // function of this class we have to + // allocate and release vectors over + // and over again. + // + // The natural question is then: + // Wouldn't it be nice if we could + // avoid this, and allocate vectors + // only once? In fact, deal.II offers + // a way to do exactly this. What all + // the linear solvers do is not to + // allocate memory using ``new'' and + // ``delete'', but rather to allocate + // them from an object derived from + // the ``VectorMemory'' class (see + // the module on Vector memory + // management in the API reference + // manual). By default, the linear + // solvers use a derived class + // ``PrimitiveVectorMemory'' that, + // ever time a vector is requested, + // allocates one using ``new'', and + // calls ``delete'' on it again once + // the solver returns it to the + // ``PrimitiveVectorMemory'' + // object. This is the appropriate + // thing to do if we do not + // anticipate that the vectors may be + // reused any time soon. + // + // On the other hand, for the present + // case, we would like to have a + // vector memory object that + // allocates vectors when asked by a + // linear solver, but when the linear + // solver returns the vectors, the + // vector memory object holds on to + // them for later requests by linear + // solvers. The + // ``GrowingVectorMemory'' class does + // exactly this: when asked by a + // linear solver for a vector, it + // first looks whether it has unused + // ones in its pool and if so offers + // this vector. If it doesn't, it + // simply grows its pool. Vectors are + // only returned to the C++ runtime + // memory system once the + // ``GrowingVectorMemory'' object is + // destroyed itself. + // + // What we therefore need to do is + // have the present matrix have an + // object of type + // ``GrowingVectorMemory'' as a + // member variable and use it + // whenever we create a linear solver + // object. There is a slight + // complication here: Since the + // ``vmult'' function is marked as + // ``const'' (it doesn't change the + // state of the object, after all, + // and simply operates on its + // arguments), it can only pass an + // unchanging vector memory object to + // the solvers. The solvers, however, + // do change the state of the vector + // memory object, even though this + // has no impact on the actual state + // of the inverse matrix object. The + // compiler would therefore flag any + // such attempt as an error, if we + // didn't make use of a rarely used + // feature of C++: we mark the + // variable as ``mutable''. What this + // does is to allow us to change a + // member variable even from a + // ``const'' member function. template class InverseMatrix : public Subscriptor { @@ -763,6 +894,19 @@ InverseMatrix::InverseMatrix (const Matrix &m) {} + // Here now is the function that + // implements multiplication with the + // inverse matrix by calling a CG + // solver. Note how we pass the + // vector memory object discussed + // above to the linear solver. Note + // also that we set the solution + // vector to zero before starting the + // solve, since we do not want to use + // the possible previous and unknown + // content of that variable as + // starting vector for the linear + // solve: template void InverseMatrix::vmult (Vector &dst, const Vector &src) const @@ -788,14 +932,14 @@ class SchurComplement : public Subscriptor 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), -- 2.39.5