debug_iterator.h

Go to the documentation of this file.

/*
 * Copyright (c) 2014-2015:  G-CSC, Goethe University Frankfurt
 * Author: Arne Nägel
 * 
 * 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__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__DEBUG_ITERATOR__
#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__DEBUG_ITERATOR__
 
#include "lib_algebra/operator/interface/preconditioner.h"
#include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
#include "lib_algebra/algebra_common/sparsematrix_util.h"
#include "lib_algebra/operator/debug_writer.h"
 
 
#ifdef UG_PARALLEL
  #include "lib_algebra/parallelization/parallelization.h"
  #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
#endif
 
namespace ug{
 
 
template <typename TAlgebra>
class DebugIterator :
public virtual ILinearIterator<typename TAlgebra::vector_type>,
public VectorDebugWritingObject <typename TAlgebra::vector_type>
{
  public:
  //  Algebra type
    typedef TAlgebra algebra_type;
 
  //  Vector type
    typedef typename TAlgebra::vector_type vector_type;
 
  //  Matrix type
    typedef typename TAlgebra::matrix_type matrix_type;
 
  //  Base type
    typedef ILinearIterator<typename TAlgebra::vector_type> base_type;
 
    typedef IPreconditionedLinearOperatorInverse<vector_type> solver_type;
  private:
    typedef VectorDebugWritingObject <typename TAlgebra::vector_type> vdwo_type;
 
  public:
    DebugIterator() : from(0.0), to (1.0)
    {
      MatrixOperator<matrix_type, vector_type> *ptr =new MatrixOperator<matrix_type, vector_type>();
      m_pOperator  = SmartPtr<MatrixOperator<matrix_type, vector_type> >(ptr);
    };
 
    DebugIterator(SmartPtr<base_type> pprecond, SmartPtr<solver_type> psolver) : from(0.0), to (1.0)
    {
      set_preconditioner(pprecond);
      set_solver(psolver);
      MatrixOperator<matrix_type, vector_type> *ptr =new MatrixOperator<matrix_type, vector_type>();
      m_pOperator  = SmartPtr<MatrixOperator<matrix_type, vector_type> >(ptr);
    };
 
    virtual SmartPtr<ILinearIterator<vector_type> > clone()
    {
      SmartPtr<DebugIterator<algebra_type> > newInst(new DebugIterator<algebra_type> ());
      newInst->set_damp(this->damping());
      newInst->set_preconditioner(this->get_preconditioner()->clone());  // clone preconditioner
      newInst->set_solver(this->get_solver());                 // keep solver (initialized below)
      newInst->set_debug(this->vdwo_type::vector_debug_writer());
      newInst->from = this->from;
      newInst->to = this->to;
      return newInst;
    }
 
    virtual bool supports_parallel() const
    {
      if(m_pprecond.valid())
        return m_pprecond->supports_parallel();
      return true;
    }
 
    void set_preconditioner(SmartPtr<base_type> pprecond) {m_pprecond=pprecond;}
    void set_solver(SmartPtr<solver_type> psolver) {m_solver=psolver;}
    void set_random_bounds(double a, double b){from =a; to=b;}
 
    void set_solution(SmartPtr<vector_type> sol)
    {m_spSolVector = sol;}
 
  protected:
    SmartPtr<base_type> get_preconditioner() {return m_pprecond;}
    SmartPtr<solver_type> get_solver() {return m_solver;}
 
 
    virtual const char* name() const
    {return "DebugIterator";}
 
 
    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);
    }
 
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > L)
    {
      //  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> > pOp)
    {
      m_pOperator->get_matrix() = pOp->get_matrix();
      //m_pOperator->get_matrix().set_as_copy_of(pOp->get_matrix());
#ifdef UG_PARALLEL
      //m_pOperator->get_matrix().set_storage_Type(pOp->get_matrix()->get_storage_type());
#endif
      m_pprecond->init(pOp);
 
      //write_debug(pOp->get_matrix(), "DebugMatrix");
      find_smooth_error();
 
      return true;
    }
 
 
    bool find_smooth_error()
    {
      if (m_solver.invalid())
      {
        UG_LOG("WARNING: cannot find smooth error; no solver supplied!");
        return false;
      }
 
      if (m_spSolVector.invalid())
      {
        UG_LOG("WARNING: cannot find smooth error; no valid vector !");
        return false;
      }
 
      //vector_type myRhs(m_pOperator->get_matrix().num_rows());
      //vector_type myError(m_pOperator->get_matrix().num_rows());
      
      SmartPtr<vector_type> myError = m_spSolVector;
      SmartPtr<vector_type> myRhs = m_spSolVector->clone_without_values();
 
      myError->set(0.0);
      myError->set_random(from, to);
 
      myRhs->set(0.0);
 
      this->write_debug(*myError, "DebugIterError0");
      m_solver->init(m_pOperator);
      m_solver->apply(*myError, *myRhs);
 
      this->write_debug(*myError, "DebugIterErrorS");
 
      return true;
    }
 
    virtual bool apply(vector_type& c, const vector_type& d)
    {
      return m_pprecond->apply(c, d);
    }
 
    virtual bool apply_update_defect(vector_type& c, vector_type& d)
    {
      return (m_pprecond->apply_update_defect(c, d));
    }
 
  protected:
 
    SmartPtr<base_type> m_pprecond;
 
    SmartPtr<solver_type> m_solver;
    SmartPtr<MatrixOperator<matrix_type, vector_type> > m_pOperator;
    SmartPtr<vector_type> m_spSolVector;
 
    double from, to;
 
};
 
 
 
} // end namespace ug
 
#endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__TRANSFORMING__