From 7b77102db4c697e7d1b3e3325453f5c7300e580e Mon Sep 17 00:00:00 2001 From: kanschat Date: Fri, 17 Jun 2011 23:23:15 +0000 Subject: [PATCH] undo erroneous commit git-svn-id: https://svn.dealii.org/trunk@23843 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-20/step-20.cc | 92 ++++++++++++++++++++++++++--- 1 file changed, 83 insertions(+), 9 deletions(-) diff --git a/deal.II/examples/step-20/step-20.cc b/deal.II/examples/step-20/step-20.cc index c24995dbfc..73646d2e33 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,14 +901,7 @@ void MixedLaplaceProblem::solve () m_inverse.solver.select("cg"); static ReductionControl inner_control(1000, 0., 1.e-13); m_inverse.solver.set_control(inner_control); - - SchurComplement >, - SparseMatrix, - SparseMatrix, - SparseMatrix > - 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 @@ -853,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); -- 2.39.5