partialDirectPreconditioner.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2012 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef PARTIAL_DIRECT_PRECONDITIONER_HH
#define PARTIAL_DIRECT_PRECONDITIONER_HH
 
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvercategory.hh>
 
#include "fem/istlinterface.hh"
#include "utilities/duneInterface.hh"
 
namespace Kaskade
{
 
  template <class Op>
  class PartialDirectPreconditioner: public Dune::Preconditioner<typename Op::Domain, typename Op::Range>
  {
    typedef typename Op::Domain Domain;
    typedef typename Op::Domain Range;
    typedef typename GetScalar<Domain>::type Scalar;
 
  public:
    static int const category = Dune::SolverCategory::sequential;
 
    PartialDirectPreconditioner(Op const& op, size_t first, size_t last,
                                DirectType directType=DirectType::UMFPACK, MatrixProperties matprop=MatrixProperties::SYMMETRIC)
    {
      MatrixAsTriplet<Scalar> A = op.template get<MatrixAsTriplet<Scalar> >();
 
      // Delete all the off-diagonal entries whose indices are not
      // contained in [first,last)x[first,last).
      size_t out = 0;
      for (size_t in=0; in<A.nnz(); ++in)
        if (A.ridx[in]==A.cidx[in] || (A.ridx[in]>=first && A.ridx[in]<last && A.cidx[in]>=first && A.cidx[in]<last)) {
          A.ridx[out] = A.ridx[in];
          A.cidx[out] = A.cidx[in];
          A.data[out] = A.data[in];
          ++out;
        }
      A.ridx.erase(A.ridx.begin()+out,A.ridx.end());
      A.cidx.erase(A.cidx.begin()+out,A.cidx.end());
      A.data.erase(A.data.begin()+out,A.data.end());
 
      solver.reset(getFactorization(directType,A,typename Factorization<Scalar>::Options(matprop)).release());
    }
 
    virtual void pre(Domain&, Range&) {}
    virtual void post (Domain&) {}
 
    virtual void apply (Domain& x, Range const& y) {
      std::vector<Scalar> b(y.dim());
      y.write(b.begin());
      solver->solve(b);
      x.read(b.begin());
    }
 
  private:
    std::unique_ptr<Factorization<Scalar> > solver;
  };
 
} // namespace Kaskade
#endif