Commit 7c8d93f3 authored by prudhomm's avatar prudhomm
Browse files

added an iterative solver iteration class handler concept

first step before much better framework
parent df0dcb14
/* -*- mode: c++ -*-
This file is part of the LifeV library
Author(s): Christophe Prud'homme <christophe.prudhomme@epfl.ch>
Date: 2005-04-07
Copyright (C) 2005 EPFL
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/**
\file iteration.hpp
\author Christophe Prud'homme <christophe.prudhomme@epfl.ch>
\date 2005-04-07
*/
#ifndef __Iteration_H
#define __Iteration_H 1
#include <boost/shared_ptr.hpp>
namespace LifeV
{
/*!
\class Iteration
\brief brief description
The Iteration object calculates whether the solution has reached
the desired accuracy, or whether the maximum number of iterations
has been reached. The method \c finished() checks both convergence and
number of iterations. The method \c converged() only checks
convergence. The error \c code() method is used to determine the
return value for the iterative solver function. The \c first()
method is used to determine the first iteration of the loop.
The following notation will be used
@li \f$r\f$ the residual
@li \f$\epsilon\f$ the relative precision
@li \f$I\f$ the number of already performed iterations
@li \f$M\f$ the maximum number of iterations allowed
@author Christophe Prud'homme <christophe.prudhomme@epfl.ch>
@see
*/
template<typename Real>
class Iteration
{
public:
/** @name Typedefs
*/
//@{
/**
\brief Numerical Type
*/
typedef Real NumericalType;
typedef Real value_type;
//@}
/** @name Constructors, destructor
*/
//@{
/**
\brief create a new instance
@param l stream for debugging purpose
*/
static Iteration<NumericalType>* New()
{
return new Iteration<NumericalType>;
}
Iteration( Iteration const& iter )
:
__iterations( iter.__iterations ),
__max_iter( iter.__max_iter ),
__residual( iter.__residual ),
__precision( iter.__precision ),
__norm_init( iter.__norm_init )
{
// do nothing here
}
//! destructor
virtual ~Iteration()
{
// do nothing here
}
//@}
/** @name Operators
*/
//@{
//! copy operator
Iteration& operator=( Iteration<NumericalType> const& iter )
{
if ( this == &iter )
{
return *this;
}
__precision = iter.__precision;
__max_iter = iter.__max_iter;
__residual = iter.__residual;
__norm_init = iter.__norm_init;
return *this;
}
//! prefix ++ operator
void operator++() throw()
{
++__iterations;
}
//@}
/** @name Accessors
*/
//@{
/**
\brief get the Residual
*/
int numberOfInterations() const { return __iterations; }
double residual() const throw()
{
return __residual;
}
double relativePrecision() const
{
return __precision;
}
int maximumNumberOfIterations() const
{
return __max_iter;
}
NumericalType initialResidual() const
{
return __norm_init;
}
value_type relaxation () const
{
return _M_relaxation;
}
/** @name Mutators
*/
//@{
//! set the Max number of iterations
/**
@param m max number of iterations to perform
*/
void setMaximumNumberOfIterations( int m ) throw()
{
__max_iter = m;
}
//! set the relative precision to reach
/**
@param p precision
*/
void setRelativePrecision( NumericalType p ) throw()
{
__precision = p;
}
//! initial norm for the residual
/*!
\param ninit initial norm for the residual
*/
void setInitialResidual( NumericalType ninit ) throw()
{
__norm_init = ninit;
}
void setRelaxation ( value_type __w )
{
_M_relaxation = __w;
}
//@}
/** @name Methods
*/
//@{
/**
\brief tells if the iteration finished
Three cases can occur:
@li if \f$ r < \epsilon \f$ then the iterration is over
@li if \f$ I > M \f$ then the iterration is over
@li else the iteration must continue
@param r residual to test the convergence
@return false if not finished and true otherwise
*/
bool isFinished( NumericalType r ) //throw(SExceptionSolverHasNotConverged)
{
bool ret;
if ( this->isConverged( r ) )
{
ret = true;
}
else if ( __iterations < __max_iter)
{
ret = false;
}
else
{
#if 0
handleEvents( true );
std::string why = "Solver has not converged";
SExceptionSolverHasNotConverged __e( why.c_str(), LOCATION );
__e.setNumberOfIterations( __iterations );
__e.setResidual( __residual );
throw __e;
#endif
}
handleEvents( ret );
return ret;
}
template<typename VectorX>
bool isFinished(const VectorX& r) //throw(SExceptionSolverHasNotConverged)
{
bool ret;
if ( this->isConverged( r ) )
{
ret = true;
}
else if ( __iterations < __max_iter)
{
ret = false;
}
else
{
#if 0
handleEvents( true );
std::string why = "Solver has not converged";
SExceptionSolverHasNotConverged __e( why.c_str(), LOCATION );
__e.setNumberOfIterations( __iterations );
__e.setResidual( __residual );
throw __e;
#endif
}
handleEvents( ret );
return ret;
}
bool isConverged( NumericalType r ) throw()
{
__residual = r / __norm_init;
return ( __residual <= __precision );
}
template<typename VectorX> bool isConverged( VectorX const& x ) throw()
{
__residual = ublas::norm_2( x ) / __norm_init;
return ( __residual <= __precision );
}
bool isFirst() const
{
return ( __iterations == 0 );
}
void reset()
{
__iterations = 0;
}
//@}
protected:
/**
Default constructor.
*/
Iteration()
:
__iterations( 0 ),
__max_iter( 0 ),
__residual( 0 ),
__precision( 0 ),
__norm_init( 1.0 ),
_M_relaxation( 1.0 )
{
// do nothing here
}
virtual void handleEvents( bool __is_finished )
{
if ( __iterations == 0 )
{
std::cout << "iteration " << __iterations << " : " << residual() << "\n";
}
if ( __is_finished == true )
{
std::cout << "iteration " << __iterations << " : " << residual() << "\n";
}
else
{
std::cout << "iteration " << __iterations << " : " << residual() << "\n";
}
}
private:
int __iterations;
int __max_iter;
value_type __residual;
value_type __precision;
value_type __norm_init;
value_type _M_relaxation;
};
typedef Iteration<double> iteration_type;
typedef boost::shared_ptr<iteration_type> iteration_ptrtype;
}
#endif /* __Iteration_H */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment