iter.h

Go to the documentation of this file.

#ifndef ITER_H
#define ITER_H
 
#include <iomanip>
 
namespace algebra
{
enum algoStatus
    {
    UNDEFINED = -1,
    CONVERGED = 0,
    ITER_OVERFLOW = 1,
    CANNOT_CONVERGE = 2
    };
 
template <typename T>
class iteration
    {
    protected:
 
    T rhsn;
 
    const int maxiter;
 
    const bool noise;
 
    int nit;
 
    T res;
 
    bool written;
 
    private:
    const std::string solver_name;
 
    std::string str_status(void) const
        {
        std::string s;
        switch(status)
            {
            case UNDEFINED:
                s="UNDEFINED";
            break;
            case CONVERGED:
                s="CONVERGED";
            break;
            case ITER_OVERFLOW:
                s="ITER_OVERFLOW";
            break;
            case CANNOT_CONVERGE:
                s="CANNOT_CONVERGE";
            break;
            }
        return s;
        }
 
    public :
    iteration(const std::string name, T r, bool _noise, int _maxiter): rhsn(1.0), maxiter(_maxiter),
        noise(_noise), nit(0), res(std::numeric_limits<T>::max()), written(false), solver_name(name), resmax(r)
        { status = UNDEFINED; }
 
    algebra::algoStatus status;
 
    const T resmax;
 
    void reset(void)
        {
        rhsn=1.0;
        nit = 0;
        res = std::numeric_limits<T>::max();
        written = false;
        status = UNDEFINED;
        }
 
    std::string infos(void) const
        {
        std::stringstream sstr;
        sstr << solver_name << " status " << str_status() << " after " << nit << " iterations, residu= " << res;
        return sstr.str();
        }
 
    void operator ++(int)
        {
        nit++;
        written = false;
        if (nit >= maxiter)
            { status = ITER_OVERFLOW; }
        }
 
    void operator ++() { (*this)++; }
 
    T get_res() const { return res; }
 
    int get_iteration() const { return nit; }
 
    T get_rhsnorm() const { return rhsn; }
 
    void set_rhsnorm(T r) { rhsn = r; }
 
    bool converged(T nr)
        {
        bool cv(false);
        res = std::fabs(nr);// fabs should be useless
        if (std::isnan(res))
            { status = CANNOT_CONVERGE; }
        else
            {
            cv = (res <= rhsn * resmax);
            if(cv)
                { status = CONVERGED; }
            }
        return cv;
        }
 
    bool finished(T nr)
        {
        if (noise && !written)
            {
            T a = (rhsn == 0) ? 1.0 : rhsn;
            converged(nr);
            std::cout << " iter " << std::setw(3) << nit << " residual " << std::setw(12) << std::fabs(nr) / a << std::endl;
            written = true;
            }
        return converged(nr);
        }
    };// end class
}//end namespace
#endif