preconditioner.h

Go to the documentation of this file.

/*
 * Copyright (c) 2013-2015:  G-CSC, Goethe University Frankfurt
 * Author: Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__
#define __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__
 
#include "linear_iterator.h"
#include "matrix_operator.h"
#include "lib_algebra/operator/damping.h"
#include "lib_algebra/operator/debug_writer.h"
 
#ifdef UG_PARALLEL
  #include "lib_algebra/parallelization/parallelization.h"
#endif
 
namespace ug{
// Matrix Iterator Operator (Preconditioner)
 
 
template <typename TAlgebra>
class IPreconditioner :
  public ILinearIterator<typename TAlgebra::vector_type>,
  public DebugWritingObject<TAlgebra>
{
  public:
    typedef TAlgebra algebra_type;
 
    typedef typename TAlgebra::vector_type vector_type;
 
    typedef typename TAlgebra::matrix_type matrix_type;
 
    typedef MatrixOperator<matrix_type, vector_type> matrix_operator_type;
    using DebugWritingObject<TAlgebra>::set_debug;
  protected:
    using ILinearIterator<vector_type>::damping;
 
    using DebugWritingObject<TAlgebra>::debug_writer;
    using DebugWritingObject<TAlgebra>::write_debug;
 
  public:
    IPreconditioner() :
      m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
    {};
 
    IPreconditioner(SmartPtr<IDebugWriter<algebra_type> > spDebugWriter) :
      DebugWritingObject<TAlgebra>(spDebugWriter),
      m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
    {};
 
    IPreconditioner( const IPreconditioner<TAlgebra> &parent ) :
      ILinearIterator<vector_type>(parent),
      DebugWritingObject<TAlgebra>(parent),
      m_spDefectOperator(NULL), m_spApproxOperator(NULL), m_bInit(false), m_bOtherApproxOperator(false)
    {
    }
  protected:
 
    virtual const char* name() const = 0;
 
 
    virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp) = 0;
 
 
    virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)  = 0;
 
    virtual bool postprocess() = 0;
 
  public:
 
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > J,
                      const vector_type& u)
    {
    //  cast to matrix based operator
      SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
          J.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
 
    //  Check that matrix if of correct type
      if(pOp.invalid())
        UG_THROW(name() << "::init': Passed Operator is "
            "not based on matrix. This Preconditioner can only "
            "handle matrix-based operators.");
 
    //  forward request to matrix based implementation
      return init(pOp);
    }
 
 
    bool init(SmartPtr<ILinearOperator<vector_type> > L)
    {
    //  Remember operator
      m_spDefectOperator = L;
      if(m_bOtherApproxOperator) return true;
 
    //  cast to matrix based operator
      SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
          L.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
 
    //  Check that matrix if of correct type
      if(pOp.invalid())
        UG_THROW(name() << "::init': Passed Operator is "
            "not based on matrix. This Preconditioner can only "
            "handle matrix-based operators.");
 
    //  forward request to matrix based implementation
      return init(pOp);
 
    }
 
 
 
    bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
    {
    //  Remember operator
      m_spApproxOperator = Op;
      m_spDefectOperator = Op;
 
    //  Check that matrix exists
      if(m_spApproxOperator.invalid())
        UG_THROW(name() << "::init': Passed Operator is invalid.");
 
    //  Preprocess
      if(!preprocess(m_spApproxOperator))
      {
        UG_LOG("ERROR in '"<<name()<<"::init': Preprocess failed.\n");
        return false;
      }
 
    //  Remember, that operator has been initialized
      m_bInit = true;
 
    //  we're done
      return true;
    }
 
 
    virtual bool apply(vector_type& c, const vector_type& d)
    {
    //  Check that operator is initialized
      if(!m_bInit)
      {
        UG_LOG("ERROR in '"<<name()<<"::apply': Iterator not initialized.\n");
        return false;
      }
 
    //  Check parallel status
      #ifdef UG_PARALLEL
      if(!d.has_storage_type(PST_ADDITIVE))
        UG_THROW(name() << "::apply: Wrong parallel "
                       "storage format. Defect must be additive.");
      #endif
 
    //  Check sizes
      THROW_IF_NOT_EQUAL_4(c.size(), d.size(),
          m_spApproxOperator->num_rows(), m_spApproxOperator->num_cols());
 
    //  apply iterator: c = B*d
      if(!step(m_spApproxOperator, c, d))
      {
        UG_LOG("ERROR in '"<<name()<<"::apply': Step Routine failed.\n");
        return false;
      }
 
    //  apply scaling
      const number kappa = damping()->damping(c, d, m_spApproxOperator);
      if(kappa != 1.0){
        c *= kappa;
      }
 
    //  Correction is always consistent
      #ifdef UG_PARALLEL
      if(!c.change_storage_type(PST_CONSISTENT))
        UG_THROW(name() << "::apply': Cannot change "
            "parallel storage type of correction to consistent.");
      #endif
 
    //  we're done
      return true;
    }
 
 
    virtual bool apply_update_defect(vector_type& c, vector_type& d)
    {
    //  compute new correction
      if(!apply(c, d)) return false;
 
    //  update defect d := d - A*c
      m_spDefectOperator->apply_sub(d, c);
 
    //  we're done
      return true;
    }
 
    virtual void set_approximation(SmartPtr<MatrixOperator<matrix_type,vector_type> > approx)
    {
      UG_COND_THROW(!approx.valid(), "");
      m_spApproxOperator = approx;
      if(!preprocess(m_spApproxOperator))
      {
        UG_THROW("ERROR in '"<<name()<<"::init': Preprocess failed.\n");
        return;
      }
      m_bOtherApproxOperator = true;
    }
 
    virtual ~IPreconditioner() {};
 
    SmartPtr<MatrixOperator<matrix_type, vector_type> > defect_operator()
    {
      return m_spDefectOperator;
    }
 
    SmartPtr<MatrixOperator<matrix_type, vector_type> > approx_operator()
    {
      return m_spApproxOperator;
    }
 
  protected:
    SmartPtr<ILinearOperator<vector_type> > m_spDefectOperator;
    SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spApproxOperator;
 
 
    bool m_bInit;
 
    bool m_bOtherApproxOperator;
};
 
 
 
}  // end namespace ug
#endif /* __H__LIB_ALGEBRA__OPERATOR__INTERFACE__PRECONDITIONER__ */