From: Guido Kanschat Date: Fri, 2 Jul 2010 20:30:37 +0000 (+0000) Subject: use library class instead of hand woven thingy X-Git-Tag: v8.0.0~5854 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e4bf9ac69e67debe833dc12fa39c9a5dcd2fc8d6;p=dealii.git use library class instead of hand woven thingy git-svn-id: https://svn.dealii.org/trunk@21452 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-20/doc/intro.dox b/deal.II/examples/step-20/doc/intro.dox index 0f992b7362..91c7b3ca77 100644 --- a/deal.II/examples/step-20/doc/intro.dox +++ b/deal.II/examples/step-20/doc/intro.dox @@ -397,46 +397,7 @@ matrix-vector products with it. We can do so in steps: to compute $Sv$, we positive definite and symmetric mass matrix $M$;
  • form $z=By$ to obtain $Sv=z$. -We will implement a class that does that in the program. Before showing its -code, let us first note that we need to multiply with $M^{-1}$ in several -places here: in multiplying with the Schur complement $S$, forming the right -hand side of the first equation, and solving in the second equation. From a -coding viewpoint, it is therefore appropriate to relegate such a recurring -operation to a class of its own. We call it InverseMatrix. As far as -linear solvers are concerned, this class will have all operations that solvers -need, which in fact includes only the ability to perform matrix-vector -products; we form them by using a CG solve (this of course requires that the -matrix passed to this class satisfies the requirements of the CG -solvers). Here are the relevant parts of the code that implements this: - -@code -class InverseMatrix -{ - public: - InverseMatrix (const SparseMatrix &m); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > matrix; - // ... -}; - - -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const -{ - SolverControl solver_control (src.size(), 1e-8*src.l2_norm()); - SolverCG<> cg (solver_control, vector_memory); - - cg.solve (*matrix, dst, src, PreconditionIdentity()); -} -@endcode - -Once created, objects of this class can act as matrices: they perform -matrix-vector multiplications. How this is actually done is irrelevant to the -outside world. +This is accomplished by using the class IterativeInverse. Using this class, we can then write a class that implements the Schur complement in much the same way: to act as a matrix, it only needs to offer a diff --git a/deal.II/examples/step-20/doc/results.dox b/deal.II/examples/step-20/doc/results.dox index bf1d2e814e..4ea4ae7634 100644 --- a/deal.II/examples/step-20/doc/results.dox +++ b/deal.II/examples/step-20/doc/results.dox @@ -13,7 +13,7 @@ examples/step-20> make run Number of active cells: 64 Total number of cells: 85 Number of degrees of freedom: 208 (144+64) -10 CG Schur complement iterations to obtain convergence. +15 CG Schur complement iterations to obtain convergence. Errors: ||e_p||_L2 = 0.178055, ||e_u||_L2 = 0.0433435 @endcode diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 36f38faf29..febeb0602f 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2005, 2006, 2007, 2008 by the deal.II authors */ +/* Copyright (C) 2005, 2006, 2007, 2008, 2010 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -28,6 +28,13 @@ #include #include #include + // For our Schur complement solver, + // we need two new objects. One is a + // matrix object which acts as the + // inverse of a matrix by calling an + // iterative solver. +#include + #include #include #include @@ -719,121 +726,6 @@ void MixedLaplaceProblem::assemble_system () // 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 and we - // don't even have to do anything - // special about it (so this comment - // is purely educational). 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 - // GrowingVectorMemory - // that, every time a vector is - // requested, allocates one from a - // pool that is shared by all - // GrowingVectorMemory - // objects. -template -class InverseMatrix : public Subscriptor -{ - public: - InverseMatrix (const Matrix &m); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer matrix; -}; - - -template -InverseMatrix::InverseMatrix (const Matrix &m) - : - matrix (&m) -{} - - - // Here now is the function that - // implements multiplication with the - // inverse matrix by calling a CG - // solver. Note 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 -{ - SolverControl solver_control (src.size(), 1e-8*src.l2_norm()); - SolverCG<> cg (solver_control); - - dst = 0; - - cg.solve (*matrix, dst, src, PreconditionIdentity()); -} - // @sect4{The SchurComplement class template} @@ -882,21 +774,21 @@ class SchurComplement : public Subscriptor { public: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix > &Minv); + const IterativeInverse > &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 IterativeInverse > &Minv) : system_matrix (&A), m_inverse (&Minv), @@ -920,7 +812,7 @@ void SchurComplement::vmult (Vector &dst, // and preconditioner system is the // class that approximates the Schur // complement so we can form a - // InverseMatrix@ + // an InverseIterate // object that approximates the // inverse of the Schur // complement. It follows the same @@ -932,6 +824,12 @@ void SchurComplement::vmult (Vector &dst, // step. Consequently, the class also // does not have to store a pointer // to an inverse mass matrix object. + // + // Since InverseIterate follows the + // standard convention for matrices, + // we need to provide a + // Tvmult function here as + // well. class ApproximateSchurComplement : public Subscriptor { public: @@ -939,6 +837,8 @@ class ApproximateSchurComplement : public Subscriptor void vmult (Vector &dst, const Vector &src) const; + void Tvmult (Vector &dst, + const Vector &src) const; private: const SmartPointer > system_matrix; @@ -964,6 +864,15 @@ void ApproximateSchurComplement::vmult (Vector &dst, } +void ApproximateSchurComplement::Tvmult (Vector &dst, + const Vector &src) const +{ + system_matrix->block(1,0).Tvmult (dst, tmp2); + system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1); + system_matrix->block(0,1).Tvmult (tmp1, src); +} + + // @sect4{MixedLaplace::solve} @@ -985,8 +894,14 @@ void ApproximateSchurComplement::vmult (Vector &dst, template void MixedLaplaceProblem::solve () { - const InverseMatrix > - m_inverse (system_matrix.block(0,0)); + PreconditionIdentity identity; + IterativeInverse > + m_inverse; + m_inverse.initialize(system_matrix.block(0,0), identity); + m_inverse.solver.select("cg"); + ReductionControl inner_control(1000, 0., 1.e-13); + m_inverse.solver.control = inner_control; + Vector tmp (solution.block(0).size()); // Now on to the first @@ -1016,8 +931,11 @@ void MixedLaplaceProblem::solve () ApproximateSchurComplement approximate_schur_complement (system_matrix); - InverseMatrix - preconditioner (approximate_schur_complement); + IterativeInverse > + preconditioner; + preconditioner.initialize(approximate_schur_complement, identity); + preconditioner.solver.select("cg"); + preconditioner.solver.control = inner_control; SolverControl solver_control (solution.block(1).size(),