From: Guido Kanschat Date: Mon, 5 Jul 2010 18:16:56 +0000 (+0000) Subject: Oops! X-Git-Tag: v8.0.0~5851 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c1efa1974511e220f7a0c8af0c74d4ceb5bbe999;p=dealii.git Oops! git-svn-id: https://svn.dealii.org/trunk@21455 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index 6cc260a81b..febeb0602f 100644 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@ -34,7 +34,6 @@ // inverse of a matrix by calling an // iterative solver. #include -#include #include #include @@ -727,6 +726,86 @@ void MixedLaplaceProblem::assemble_system () // rather only comment on // implementational aspects. + + // @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. + // + // The vmult function requires + // two temporary vectors that we do + // not want to re-allocate and free + // every time we call this + // function. Since here, we have full + // control over the use of these + // vectors (unlike above, where a + // class called by the vmult + // function required these vectors, + // not the vmult function + // itself), we allocate them + // directly, rather than going + // through the VectorMemory + // mechanism. However, again, these + // member variables do not carry any + // state between successive calls to + // the member functions of this class + // (i.e., we never care what values + // they were set to the last time a + // member function was called), we + // mark these vectors as mutable. + // + // The rest of the (short) + // implementation of this class is + // straightforward if you know the + // order of matrix-vector + // multiplications performed by the + // vmult function: +class SchurComplement : public Subscriptor +{ + public: + SchurComplement (const BlockSparseMatrix &A, + const IterativeInverse > &Minv); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer > system_matrix; + const SmartPointer > > m_inverse; + + mutable Vector tmp1, tmp2; +}; + + +SchurComplement::SchurComplement (const BlockSparseMatrix &A, + const IterativeInverse > &Minv) + : + system_matrix (&A), + m_inverse (&Minv), + tmp1 (A.block(0,0).m()), + tmp2 (A.block(0,0).m()) +{} + + +void SchurComplement::vmult (Vector &dst, + const Vector &src) const +{ + system_matrix->block(0,1).vmult (tmp1, src); + m_inverse->vmult (tmp2, tmp1); + system_matrix->block(1,0).vmult (dst, tmp2); +} + + // @sect4{The ApproximateSchurComplement class template} // The third component of our solver @@ -822,11 +901,7 @@ void MixedLaplaceProblem::solve () m_inverse.solver.select("cg"); ReductionControl inner_control(1000, 0., 1.e-13); m_inverse.solver.control = inner_control; - - SchurComplement schur_complement (m_inverse, - system_matrix.block(1,0), - system_matrix.block(1,0)); - + Vector tmp (solution.block(0).size()); // Now on to the first @@ -850,6 +925,8 @@ void MixedLaplaceProblem::solve () schur_rhs -= system_rhs.block(1); + SchurComplement + schur_complement (system_matrix, m_inverse); ApproximateSchurComplement approximate_schur_complement (system_matrix);