#ifndef BLOCKDIAGONALPRECONDITIONER #define BLOCKDIAGONALPRECONDITIONER BLOCKDIAGONALPRECONDITIONER /** \file blockdiagonalpreconditioner.hh * \brief Sequential block-diagonal preconditioners * * Contains several block-diagonal preconditioners. In contrast to the * block-recursive preconditioners in ISTL, a depth of 1 is assumed and * the blocks are inverted with LAPACK routines or iteratively. In the * matrix-based preconditioners, the LU- or Cholesky factorisation is * precomputed and stored to avoid a repeated factorisation in every * iteration which is currently done in the ISTL block-preconditioners. */ #include "inversematrix.hh" namespace Dune { /** \class SeqMatrixfreeBlockJac * \brief Matrix-free Jacobi preconditioner with iterative block-inverse * * Given an operator \f\$A\f\$ and it's block-inverse \f\$D^{-1}\f\$, * the preconditioner implements the following block-Jacobi operation * (which can be iterated multiply times): * * \f[ * v \mapsto v + wD^{-1}(d-Av) * \f] * * Here the block-inverse \f\$D^{-1}\f\$ is applied in a matrix-free * way, using the wrapper class * IterativeInverseBlockDiagonalLocalOperatorWrapper. The operator \f\$A\f\$ * is also applied in a matrix-free way. */ template class SeqMatrixfreeBlockJac : public Dune::Preconditioner { public: /** \brief preconditioner category */ enum { category=SolverCategory::sequential }; /** \brief Domain type */ typedef typename GO::Traits::Domain domain_type; /** \brief Range type */ typedef typename GO::Traits::Domain range_type; /** \brief matrixfree operator */ typedef Dune::PDELab::OnTheFlyOperator OP; /** \brief matrixfree operator for inverse diagonal */ typedef Dune::PDELab::OnTheFlyOperator OPInvDiag; /** \brief create a new instance * * \param[in] go_ Grid operator \f\$A\f\$ * \param[in] goinvdiag_ Inverse diagonal gridoperator \f\$D^{-1}\f\$ * \param[in] niter_ Number of block-Jacobi iterations * \param[in] weight_ Scalar weight \f\$w\f\$ */ SeqMatrixfreeBlockJac(const GO& go_, const GOInvDiag& goinvdiag_, const size_t niter_, const double weight_) : niter(niter_), weight(weight_), op(go_), opinvdiag(goinvdiag_) {} /** \brief Prepare preconditioner (do nothing) * * \param[inout] v Output vector * \param[inout] d Input vector */ virtual void pre(domain_type& v, range_type& d) {} /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * * If the solution vector is zero in the first application we * can avoid one operator application. * * \param[inout] v output vector * \param[in] d input vector * \param[in] solution_is_zero Is the input solution vector zero? */ virtual void apply (domain_type& v, const range_type& d, const bool solution_is_zero) { range_type r(d); // r <- d for (int i=0;i0) r = d; // Reassign d if ( (i==0) and (solution_is_zero) ) { r *= weight; opinvdiag.apply(r,v); } else { op.applyscaleadd(-1.0,v,r); // r <- d - A.v opinvdiag.applyscaleadd(weight,r,v); // v <- v + w.D^{-1}.r } } } /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * * \param[inout] v output vector * \param[in] d input vector */ virtual void apply (domain_type& v, const range_type& d) { apply(v,d,false); } /** \brief Clean up after preconditioner application (do nothing) * * \param[inout] v Input vector */ virtual void post(domain_type& v) {} private: /** \brief Number of iterations */ const size_t niter; /** \brief Jacobi weight \f\$w\f\$ */ const double weight; /** \brief Operator \f\$A\f\$ */ const OP op; /** \brief Inverse diagonal part \f\$D^{-1}\f\$ of operator \f\$A\f\$ */ const OPInvDiag opinvdiag; }; /** \class SeqPartiallyMatrixfreeBlockJac * \brief Partially matrix-free Jacobi preconditioner with explicit * block-inverse * * Given an operator \f\$A\f\$ and the matrix-representation of it's * block-diagonal \f\$D\f\$, this preconditioner implements the following * block-Jacobi operation (which can be iterated multiply times): * * \f[ * v \mapsto v + wD^{-1}(d-Av) * \f] * * In contrast to SeqIterativeBlockJac, the block-inverse \f\$D^{-1}\f\$ is * assembled, factorised and stored explicitly. The operator \f\$A\f\$ is * still applied in a matrix-free way. */ template class SeqPartiallyMatrixfreeBlockJac : public Dune::Preconditioner { public: /** \brief preconditioner category */ enum { category=SolverCategory::sequential }; /** \brief Domain type */ typedef typename GO::Traits::Domain domain_type; /** \brief Range type */ typedef typename GO::Traits::Range range_type; /** \brief matrixfree operator */ typedef Dune::PDELab::OnTheFlyOperator OP; /** \brief create a new instance * * \param[in] go_ Grid operator \f\$A\f\$ * \param[in] D_ Diagonal matrix \f\$D\f\$ * \param[in] niter_ Number of iterations * \param[in] weight_ Jacobi weight \f\$w\f\$ * \param[in] symmetric_blockdiag_ Are the block-diagonal entries * symmetric? */ SeqPartiallyMatrixfreeBlockJac(const GO& go_, const M& D_, const size_t niter_, const double weight_, const bool symmetric_blockdiag_) : op(go_), D(D_), niter(niter_), weight(weight_), symmetric_blockdiag(symmetric_blockdiag_) { // check that matrix is square assert (D.N() == D.M()); // Construct inverse of block-diagonal invD = new InverseBlockDiagonalMatrix>(PDELab::Backend::native(D),symmetric_blockdiag); } /*! \brief Destructor */ ~SeqPartiallyMatrixfreeBlockJac() { delete invD; } /** \brief Prepare preconditioner (do nothing) * * \param[inout] v Output vector * \param[inout] d Input vector */ virtual void pre(domain_type& v, range_type& d) {} /** \brief Apply preconditioner * * If the solution vector is zero in the first application we * can avoid one operator application. * * Apply preconditioner to a vector d to obtain another vector v * \param[inout] v output vector * \param[in] d input vector * \param[in] solution_is_zero Is the incoming solution zero? */ virtual void apply (domain_type& v, const range_type& d, const bool solution_is_zero) { range_type r(d); // r <- d for (int i=0;i0) r = d; // Reassign d if ( (i==0) and (solution_is_zero) ) { // r <- w*d r *= weight; // v <- w*D^{-1}.d invD->mv(PDELab::Backend::native(r), PDELab::Backend::native(v)); } else { // r <- d - A.v op.applyscaleadd(-1.0,v,r); // v <- v + w*D^{-1}.r = v + w*D^{-1}.(d-A.v) invD->usmv(weight, PDELab::Backend::native(r), PDELab::Backend::native(v)); } } } /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * * \param[inout] v output vector * \param[in] d input vector */ virtual void apply (domain_type& v, const range_type& d) { apply(v,d,false); } /** \brief Clean up after preconditioner application (do nothing) * * \param[inout] v Input vector */ virtual void post(domain_type& v) {} private: /** \brief Operator \f\$A\f\$ */ const OP op; /** \brief Matrix for diagonal part \f\$D\f\$ of operator \f\$A\f\$ */ const M& D; /** \brief Number of iterations */ const size_t niter; /** \brief Jacobi weight \f\$w\f\$ */ const double weight; //! \brief Are the block-diagonal matrices symmetric? const bool symmetric_blockdiag; // Cache for inverse block-diagonal InverseBlockDiagonalMatrix> *invD; }; /** \class SeqBlockJac * \brief Sequential matrix-based Jacobi preconditioner * * Compared to the DUNE-ISTL default implementation, the inverse of the * block-diagonal is not recomputed in every iteration. * * \tparam M The matrix type to operate on * \tparam X Type of the update * \tparam Y Type of the defect * \tparam level The block level to invert. Default is 1 */ template class SeqBlockJac : public Preconditioner { public: //! \brief The matrix type the preconditioner is for. typedef M matrix_type; //! \brief The domain type of the preconditioner. typedef X domain_type; //! \brief The range type of the preconditioner. typedef Y range_type; //! \brief The field type of the preconditioner. typedef typename X::field_type field_type; typedef typename X::field_type domain_value_type; typedef typename Y::field_type range_value_type; // define the category enum { //! \brief The category the preconditioner is part of category=SolverCategory::sequential }; /*! \brief Constructor. Constructor gets all parameters to operate the prec. \param A The matrix to operate on. \param n The number of iterations to perform. \param w The relaxation factor. \param symmetric_blockdiag Are the block-diagonal entries symmetric? */ SeqBlockJac (const M& A, int niter, field_type weight, const bool symmetric_blockdiag) : _A_(A), _niter(niter), _weight(weight), _symmetric_blockdiag(symmetric_blockdiag) { CheckIfDiagonalPresent::check(_A_); // check that matrix is square assert (_A_.N() == _A_.M()); // Construct inverse of block-diagonal _A_blockdiaginv = new InverseBlockDiagonalMatrix(_A_,symmetric_blockdiag); } /*! \brief Destructor */ ~SeqBlockJac() { delete _A_blockdiaginv; } /*! \brief Prepare the preconditioner. */ virtual void pre (X& x, Y& b) { DUNE_UNUSED_PARAMETER(x); DUNE_UNUSED_PARAMETER(b); } /*! \brief Apply the preconditioner. */ virtual void apply (X& v, const Y& d, const bool solution_is_zero) { X r(d); for (int i=0; i<_niter; i++) { if (i>0) r = d; // reassign RHS if ( (i==0) and (solution_is_zero) ) { // r = w*d r *= _weight; // v <- w*d _A_blockdiaginv->mv(r,v); // v <- w*D^{-1}.v } else { // r <- r - A.v = d - A.v _A_.mmv(v,r); // v <- v + w*D^{-1}.tmp = v + w*D^{-1}.(d-A.v) _A_blockdiaginv->usmv(_weight,r,v); } } } /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * * \param[inout] v output vector * \param[in] d input vector */ virtual void apply (domain_type& v, const range_type& d) { apply(v,d,false); } /*! \brief Clean up. */ virtual void post (X& x) { DUNE_UNUSED_PARAMETER(x); } private: //! \brief The matrix we operate on. const M& _A_; //! \brief The number of steps to perform during apply. int _niter; //! \brief The relaxation parameter to use. field_type _weight; //! \brief Are the block-diagonal matrices symmetric? const bool _symmetric_blockdiag; // Cache for inverse block-diagonal InverseBlockDiagonalMatrix *_A_blockdiaginv; }; /** \brief (partially) Matrix-free SSOR preconditioner * * Given an operator \f\$A\f\$ and it's block-inverse \f\$D^{-1}\f\$, * the preconditioner implements the following block-SSOR operation * (which can be iterated multiply times): * * \f[ * v_i \mapsto v_i += wD_{ii}^{-1}(d_i-A_{i,j}v_j) * \f] * The operator \f\$A\f\$ is always applied in a matrix-free way, but the * block-diagonal matrix might either be stored explicitly or * applied iteratively, depending on the local operator that implements * the SOR operation. */ template class SeqMatrixfreeBlockSSOR : public Dune::Preconditioner { public: /** \brief preconditioner category */ enum { category=SolverCategory::sequential }; /** \brief Domain type */ typedef typename GO::Traits::Domain domain_type; /** \brief Range type */ typedef typename GO::Traits::Range range_type; /** \brief create a new instance * * \param[in] go_ Grid operator for SOR operation \f\$A\f\$ * \param[in] niter_ Number of block-SSOR iterations */ SeqMatrixfreeBlockSSOR(const GO& go_, const size_t niter_) : niter(niter_), go(go_) {} /** \brief Prepare preconditioner (do nothing) * * \param[inout] v Output vector * \param[inout] d Input vector */ virtual void pre(domain_type& v, range_type& d) {} /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * \param[inout] v output vector * \param[in] d input vector */ virtual void apply (domain_type& v, const range_type& d, const bool solution_is_zero) { if (solution_is_zero) v=0.0; for (int i=0;i class SeqBlockSSOR : public Preconditioner { public: //! \brief The matrix type the preconditioner is for. typedef M matrix_type; //! \brief The domain type of the preconditioner. typedef X domain_type; //! \brief The range type of the preconditioner. typedef Y range_type; //! \brief The field type of the preconditioner. typedef typename X::field_type field_type; typedef typename X::field_type domain_value_type; typedef typename Y::field_type range_value_type; // define the category enum { //! \brief The category the preconditioner is part of category=SolverCategory::sequential }; /*! \brief Constructor. Constructor gets all parameters to operate the prec. \param A The matrix to operate on. \param n The number of iterations to perform. \param w The relaxation factor. \param symmetric_blockdiag Are the block-diagonal entries symmetric? */ SeqBlockSSOR (const M& A, int niter, field_type weight, const bool symmetric_blockdiag) : _A_(A), _niter(niter), _weight(weight), _symmetric_blockdiag(symmetric_blockdiag) { CheckIfDiagonalPresent::check(_A_); // check that matrix is square assert (_A_.N() == _A_.M()); // Construct inverse of block-diagonal _A_blockdiaginv = new InverseBlockDiagonalMatrix(_A_,symmetric_blockdiag); } /*! \brief Destructor */ ~SeqBlockSSOR() { delete _A_blockdiaginv; } /*! \brief Prepare the preconditioner. */ virtual void pre (X& x, Y& b) { DUNE_UNUSED_PARAMETER(x); DUNE_UNUSED_PARAMETER(b); } /*! \brief Apply the preconditioner. */ virtual void apply (X& v, const Y& d, const bool solution_is_zero) { if (solution_is_zero) v=0.0; // Loop over all iterations for (int k=0; k<_niter; k++) { iterate(v,d,false); // forward iteration iterate(v,d,true); // backward iteration } } /** \brief Apply preconditioner * * Apply preconditioner to a vector d to obtain another vector v * * \param[inout] v output vector * \param[in] d input vector */ virtual void apply (domain_type& v, const range_type& d) { apply(v,d,false); } /*! \brief Clean up. */ virtual void post (X& x) { DUNE_UNUSED_PARAMETER(x); } private: /*! \brief Iterate over all rows in a specific order */ virtual void iterate (X& v, const Y& d, const bool reverse) { typedef typename M::ConstRowIterator rowiterator; typedef typename M::ConstColIterator coliterator; typedef typename Y::block_type bblock; typedef typename X::block_type xblock; bblock rhs; xblock x; rowiterator begini,endi; int increment; if (reverse) { begini=_A_.beforeEnd(); endi=_A_.beforeBegin(); increment=-1; } else { begini=_A_.begin(); endi=_A_.end(); increment=+1; } // Loop over all rows for (rowiterator i=begini; i!=endi;i+=increment) { rhs = d[i.index()]; // rhs = b_i coliterator endj=(*i).end(); // iterate over a_ij with j < i for (coliterator j=(*i).begin(); j!=endj; ++j) { if (j.index()==i.index()) continue; (*j).mmv(v[j.index()],rhs); // rhs -= sum_j a_ij * v_j } // Divide by block-diagonal auto& invdiag = (*_A_blockdiaginv)[i.index()]; invdiag.apply(rhs,x); x.axpy(-1.0,v[i.index()]); v[i.index()].axpy(_weight,x); // v_i += w / a_ii * (b_i - sum_{j=i} a_ij * vold_j) } } //! \brief The matrix we operate on. const M& _A_; //! \brief The number of steps to perform during apply. int _niter; //! \brief The relaxation parameter to use. field_type _weight; //! \brief Are the block-diagonal matrices symmetric? const bool _symmetric_blockdiag; // Cache for inverse block-diagonal InverseBlockDiagonalMatrix *_A_blockdiaginv; }; } #endif // BLOCKDIAGONALPRECONDITIONER