Commit 8be01a3c authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Global improvement of parameters

Add iteration parameter to Orthogonalizer;
parent 24ec0b74
......@@ -10,7 +10,7 @@
auto ortho = fabulous::orthogonalization(/*iteration_count=*/3)
+ fabulous::OrthoType::RUHE
+ fabulous::OrthoScheme::IMGS;
auto solution = fabulous::solve(eq) + dr + ortho + arn;
auto solution = fabulous::solve(eq, arn, ortho, dr);
fabulous::print_info(solution);
return 0;
}
......
#ifndef FABULOUS_ARNOLDI_HESSENBERG_PAIR_H
#define FABULOUS_ARNOLDI_HESSENBERG_PAIR_H
#ifndef FABULOUS_ALGO_TYPE_HPP
#define FABULOUS_ALGO_TYPE_HPP
#include "fabulous/utils/Meta.hpp"
#include <fabulous/arnoldi/ArnoldiIB.hpp>
#include <fabulous/arnoldi/ArnoldiDR.hpp>
......@@ -11,13 +13,11 @@
#include <fabulous/hessenberg/HessChamQR.hpp>
#include <fabulous/hessenberg/HessIB.hpp>
namespace fabulous {
template<class S> using ARNOLDI_STD = ArnoldiDR<HessDR, S>;
template<class S> using ARNOLDI_QR = ArnoldiDR<HessQR, S>;
template<class S> using ARNOLDI_QRDR = ArnoldiDR<HessQRDR, S>;
template<class S> using ARNOLDI_IB = ArnoldiIB<HessIB, S>;
#ifdef FABULOUS_USE_CHAMELEON
......@@ -25,6 +25,53 @@ template<class S> using ARNOLDI_CHAM_QR = ArnoldiDR<HessChamQR, S>;
template<class S> using ARNOLDI_CHAM_TL = ArnoldiDR<HessChamTLDR, S>;
#endif
inline meta::Template<ARNOLDI_STD> arnoldi_std() { return meta::Template<ARNOLDI_STD>{}; }
inline meta::Template<ARNOLDI_QR> arnoldi_qr() { return meta::Template<ARNOLDI_QR>{}; }
inline meta::Template<ARNOLDI_QRDR> arnoldi_qrdr() { return meta::Template<ARNOLDI_QRDR>{}; }
inline meta::Template<ARNOLDI_IB> arnoldi_ib() { return meta::Template<ARNOLDI_IB>{}; }
#ifdef FABULOUS_USE_CHAMELEON
inline meta::Template<ARNOLDI_CHAM_QR> arnoldi_cham_qr()
{ return meta::Template<ARNOLDI_CHAM_QR>{}; }
inline meta::Template<ARNOLDI_CHAM_TL> arnoldi_cham_tl()
{ return meta::Template<ARNOLDI_CHAM_TL>{}; }
#endif
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_STD> &)
{
return (o << "STD");
}
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_QR> &)
{
return (o << "QR");
}
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_QRDR> &)
{
return (o << "QR");
}
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_IB> &)
{
return (o << "IB");
}
#ifdef FABULOUS_USE_CHAMELEON
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_CHAM_QR> &)
{
return (o << "CHAM_QR");
}
std::ostream& operator<<(std::ostream &o, const meta::Template<ARNOLDI_CHAM_TL> &)
{
return (o << "CHAM_TL");
}
#endif
};
#endif // FABULOUS_ARNOLDI_HESSENBERG_PAIR_H
#endif // FABULOUS_ALGO_TYPE_HPP
......@@ -2,48 +2,136 @@
#define FABULOUS_BGMRES_HPP
#include <algorithm>
#include <iostream>
namespace fabulous {
template< class S > class BGMRes;
};
#include "fabulous/data/Block.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/algo/KrylovAlgo.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
namespace fabulous {
/**
* \brief %Block General Minimum Residual (%BGMRes) algorithm
*/
template< class RestartParam, class S > class BGMRes : public KrylovAlgo<S>
template< class S > class BGMRes
{
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
private:
FABULOUS_INHERITS_KRYLOV_ALGO(S);
FABULOUS_INHERITS_KRYLOV_ALGO_PROTECTED;
using P = primary_type;
Logger<P> _logger;
int _nb_restart;
int _nb_iteration;
int _mvp;
private:
template<class Matrix>
void print_precond(const Matrix &A, std::ostream &o = std::cout)
{
using fabulous::Color;
o << "Right Pre Conditionner : [";
if (A.useRightPrecond())
o<<Color::green<< "USED";
else
o<<Color::red<< "NOT USED";
o<<Color::reset<<"]\n";
}
void print_iteration_start_info(int size_to_span, std::ostream &o = std::cout)
{
o << "################## Restart number: " <<_nb_restart <<" #################\n";
o << "################## SizeToSpan : " <<size_to_span <<" ##############\n";
}
void print_iteration_end_info(int size_krylov_space, std::ostream &o = std::cout)
{
o <<"######################################################\n";
o <<" Restart Done : "<< _nb_restart <<"\n";
o <<" Total iterations : "
<< _nb_iteration<<"(+"<<_nb_restart<<")\n";
o <<" Size of Krylov Space : "<< size_krylov_space <<"\n";
o <<" Total Number of Matrix vector Product : "<< _mvp<<"\n";
o <<"######################################################\n";
}
template< class Matrix, class RestartParam >
void print_start_info(
int dim, int nbRHS, int maxMVP, int max_krylov_space_size,
const Orthogonalizer &ortho, const RestartParam &restart_param,
const std::vector<P> &epsilon, Matrix &A, std::ostream &o = std::cout)
{
Arithmetik<S> arith;
o << "#######################################################\n";
o << "#################### BGMRES #########################\n";
o << "#######################################################\n";
o << "Dimension of Problem : "<<dim<<"\n";
o << "Number of Right Hand Side : "<<nbRHS<<"\n";
o << "Maximum nb of Mat Vect product: "<<maxMVP<<"\n";
o << "Max Size Krylov space before restart : "<<max_krylov_space_size<<"\n";
o << "Orthogonalization Scheme : "<<ortho<<"\n";
o << "Arithmetic used : "<<arith<<"\n";
if (epsilon.size() == 1)
o << "Targeted Tolerance : "<<epsilon.front()<<"\n";
else
o << "Multi precision is being used\n";
print_precond(A, o);
o << "#######################################################\n";
o << restart_param;
}
public:
FABULOUS_INHERITS_KRYLOV_ALGO_PUBLIC;
std::string get_name() const override { return "Block GMRES"; }
template< template<class> class ARNOLDI_HESS, class RestartParam, class Matrix >
int solve(const Matrix &A,
int m, int n,
S *B_, int ldb,
S *X_, int ldx,
const int maxMVP,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
const RestartParam &restart_param,
OrthoScheme ortho_scheme = OrthoScheme::ICGS,
OrthoType ortho_type = OrthoType::RUHE)
BGMRes(bool log_real_residual = false):
_logger{log_real_residual},
_nb_restart{0},
_nb_iteration{0},
_mvp{0}
{
}
void reset()
{
_nb_restart = 0;
_nb_iteration = 0;
_mvp = 0;
_logger.reset();
}
template< class Algo, class RestartParam, class Matrix >
int solve( const Matrix &A,
int dim, int nbRHS, S *B_, int ldb, S *X_, int ldx,
const std::vector<P> &epsilon,
Algo algo, const int max_mvp, const int max_space,
Orthogonalizer ortho, RestartParam restart_param)
{
Block<S> B{m, n, ldb, B_};
Block<S> X{m, n, ldx, X_};
Block<S> X{dim, nbRHS, ldx, X_};
Block<S> B{dim, nbRHS, ldb, B_};
return solve_block<ARNOLDI_HESS>(
A, B, X, maxMVP, max_krylov_space_size,
epsilon, restart_param, ortho_scheme, ortho_type
return solve(
A, X, B, algo, max_mvp, max_space,
epsilon, ortho, restart_param
);
}
template<class Algo, class Equation, class RestartParam>
int solve(Equation &eq, Algo &algo, Parameters param,
Orthogonalizer ortho, RestartParam restart_param)
{
return solve(eq._A, eq._X, eq._B, algo,
param._max_mvp, param._max_space,
eq._epsilon, ortho, restart_param);
}
/**
* \brief BGMRes with Inexact Breakdowns, restart and preconditionner
*
......@@ -52,33 +140,30 @@ public:
* \param[in] A Matrix, need to provide <br/>
* MatBlockVect(), size(), useRightPreCond(), PrecondBlockVect(), DotProduct()
*
* \param[in] B Block containing Right Hand Side
* \param[in,out] X Initial guess solution and solution in ouput
* \param[in] B Block containing Right Hand Side
* \param maxMVP Maximum number of Matrix Vector Product.
* \param max_krylov_space_size maximum size of Krylov space
* \param[in] epsilon Target accuracy
* \param ortho_scheme Orthogonalization schema
* \param ortho_type Choose between BLOCK-wise arnoldi
* \param[in] restart_param DeflatedRestart or ClassicRestart instance
* \param ortho Orthogonalizer
* (through QR factorization) or vector wise arnoldi. <br/>
* (In distributed, only the vector wise will works)
*
* \return Total Number of Matrix vector product done
* (cumulated size of the spanned Krylov Spaces)
*/
template< template<class> class ARNOLDI_HESS,
class Matrix, class Block >
int solve_block( Matrix &A, Block &B, Block &X,
const int maxMVP, const int max_krylov_space_size,
const std::vector<P> &epsilon,
const RestartParam &restart_param,
OrthoScheme ortho_scheme = OrthoScheme::ICGS,
OrthoType ortho_type = OrthoType::RUHE)
template< class Algo, class RestartParam, class Matrix >
int solve( Matrix &A, Block<S> &X, Block<S> &B,
Algo &, const int maxMVP, const int max_krylov_space_size,
const std::vector<P> &epsilon,
Orthogonalizer ortho, RestartParam restart_param )
{
const int nbRHS = B.get_nb_col();
const int dim = B.get_nb_row();
reset();
print_start_info( dim, nbRHS, maxMVP, max_krylov_space_size,
ortho_scheme, ortho_type, epsilon, A );
ortho, restart_param, epsilon, A );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
inv_normB = array_inverse(normB);
......@@ -92,10 +177,10 @@ public:
int size_to_span = std::min(max_krylov_space_size, maxMVP-_mvp);
size_to_span = std::max(size_to_span, nbRHS);
print_iteration_start_info(size_to_span);
ARNOLDI_HESS<S> arnoldi{_logger, restarter, dim, nbRHS, size_to_span};
convergence = arnoldi.run( A, X, B, size_to_span, epsilon, restarter,
ortho_scheme, ortho_type );
using Arnoldi = typename Algo::template t3mpl4te<S>;
Arnoldi arnoldi{_logger, restarter, dim, nbRHS, size_to_span};
convergence = arnoldi.run( A, X, B, size_to_span, epsilon,
restarter, ortho );
_nb_iteration += arnoldi.get_nb_iteration();
_mvp += arnoldi.get_nb_mvp();
++ _nb_restart;
......@@ -104,6 +189,11 @@ public:
X.scale(normB); B.scale(normB);
return _mvp;
}
Logger<P> &get_logger()
{
return _logger;
}
};
}; // namespace fabulous
......
#ifndef FABULOUS_EQUATION_HPP
#define FABULOUS_EQUATION_HPP
#include "fabulous/utils/Arithmetic.hpp"
#include "fabulous/data/Block.hpp"
namespace fabulous {
template< class S > class BGMRes;
/**
* \brief Tuple container for A, X, B and epsilon
*/
template<class Matrix, class S>
class Equation {
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
private:
friend class BGMRes<S>;
using P = primary_type;
Matrix &_A;
Block<S> _X;
Block<S> _B;
std::vector<P> _epsilon;
public:
Equation(Matrix &A, Block<S> &X, Block<S> &B, const std::vector<P> &epsilon):
_A(A), _X(X), _B(B), _epsilon(epsilon)
{
}
Equation(Matrix &A, int dim, int nbRHS,
S *X, int ldx, S *B, int ldb,
P *epsilon, int Nepsilon):
_A(A), _X(dim, nbRHS, X, ldx), _B(dim, nbRHS, B, ldb),
_epsilon{epsilon, epsilon+Nepsilon}
{
}
};
/**
* \brief Helper function to create Equation object
*/
template<class Matrix, class S,
class P = typename Arithmetik<S>::primary_type >
inline Equation<Matrix, S> equation(
Matrix &A, Block<S> &X, Block<S> &B, std::vector<P> epsilon)
{
return Equation<Matrix,S>{A, X, B, epsilon};
}
}; // namespace fabulous
#endif // FABULOUS_EQUATION_HPP
#ifndef FABULOUS_KRYLOV_ALGO_H
#define FABULOUS_KRYLOV_ALGO_H
namespace fabulous {
template<class> class KrylovAlgo;
};
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/Arithmetic.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho_Enum.hpp"
namespace fabulous {
#define FABULOUS_INHERITS_KRYLOV_ALGO(TYPE__) \
typedef KrylovAlgo<TYPE__> super;
#define FABULOUS_INHERITS_KRYLOV_ALGO_PROTECTED \
using super::_mvp; \
using super::_nb_restart; \
using super::_nb_iteration; \
using super::_logger;
#define FABULOUS_INHERITS_KRYLOV_ALGO_PUBLIC \
using super::get_logger; \
using super::reset; \
using super::print_start_info; \
using super::print_dr_start_info; \
using super::print_iteration_start_info; \
using super::print_iteration_end_info; \
using typename super::value_type; \
using typename super::primary_type; \
using typename super::P; \
/**
* \brief %base class for BGMRes and BGMResDR
*
* This class is not intended to be instancied directly
* and shall keep its constructor as protected.
*/
template <class S> class KrylovAlgo
{
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
protected:
Logger<P> _logger;
int _nb_restart;
int _nb_iteration;
int _mvp;
private:
template<class Matrix>
void print_precond(const Matrix &A, std::ostream &o = std::cout)
{
using fabulous::Color;
o << "Right Pre Conditionner : [";
if (A.useRightPrecond())
o<<Color::green<< "USED";
else
o<<Color::red<< "NOT USED";
o<<Color::reset<<"]\n";
}
protected:
KrylovAlgo():
_logger{false},
_nb_restart(0),
_nb_iteration(0),
_mvp(0)
{
}
virtual std::string get_name() const { return "generic krylov algo"; }
void reset()
{
_nb_restart = 0;
_nb_iteration = 0;
_mvp = 0;
}
void print_iteration_end_info(int size_krylov_space)
{
std::ostream& o = std::cout;
o <<"######################################################\n";
o <<" Restart Done : "<< _nb_restart <<"\n";
o <<" Total iterations : "<< _nb_iteration<<"(+"
<<_nb_restart<<")\n";
o <<" Size of Krylov Space : "<< size_krylov_space <<"\n";
o <<" Total Number of Matrix vector Product : "<< _mvp<<"\n";
o <<"######################################################\n";
}
void print_iteration_start_info(int size_to_span)
{
std::ostream& o = std::cout;
o << "\n";
o << "################## Restart number: " <<_nb_restart <<" #################\n";
o << "################## SizeToSpan : " <<size_to_span <<" #################\n";
}
template< class Matrix >
void print_start_info(int dim,
int SizeBlock,
int MaxMVP,
int MaxKrylovSpaceSize,
OrthoScheme scheme, OrthoType type,
const std::vector<P>& epsilon,
Matrix &A)
{
std::ostream &o = std::cout;
o << "#######################################################\n";
o << "################## "<< get_name() <<" ###################\n";
o << "#######################################################\n";
o << "Dimension of Problem : "<<dim<<"\n";
o << "Number of Right Hand Side : "<<SizeBlock<<"\n";
o << "Maximum nb of Mat Vect product: "<<MaxMVP<<"\n";
o << "Max Size Krylov space before restart : "<<MaxKrylovSpaceSize<<"\n";
o << "Orthogonalization Scheme : "<<scheme<<" " << type <<"\n";
using ARI = fabulous::Arithmetik<S>;
o << "Arithmetic used : "
<< ARI::name <<" ("<<ARI::type_name<<")\n";
if (epsilon.size() == 1)
o << "Targeted Tolerance : "<<epsilon.front()<<"\n";
else
o << "Multi precision is being used\n";
print_precond(A, o);
o << "#######################################################\n";
}
template<class Scalar>
void print_dr_start_info(int NbEigenVector, Scalar target)
{
std::ostream &o = std::cout;
o << "#################### DR PARAM ####################\n";
o << "Size of recycled space at restart : "<<NbEigenVector<<"\n";
o << "Targeted eigen value for restart : "<<target<<"\n";
o << "#######################################################\n";
}
public:
Logger<P> &get_logger()
{
return _logger;
}
public:
virtual ~KrylovAlgo(){}
};
}; // namespace fabulous
#endif // FABULOUS_KRYLOV_ALGO_H
#ifndef FABULOUS_PARAMETERS_HPP
#define FABULOUS_PARAMETERS_HPP
namespace fabulous {
template< class S > class BGMRes;
/**
* \brief Hold parameters for BGMRES algorithm
*/
class Parameters
{
private:
template<class S> friend class BGMRes;
int _max_mvp; /*!< maximum number of matrix vector product */
int _max_space; /*!< maximum size of krylov search space */
public:
Parameters(int max_mvp, int max_space):
_max_mvp(max_mvp), _max_space(max_space)
{
}
};
/**
* \brief Helper function to build Parameters object
* \param max_mvp maximum number of matrix vector product
* \param max_space maximum size of krylov search space
*/
Parameters parameters(int max_mvp, int max_space)
{
return Parameters{max_mvp, max_space};
}
}; // namespace fabulous
#endif // FABULOUS_PARAMETERS_HPP
#ifndef FABULOUS_SOLVE_HPP
#define FABULOUS_SOLVE_HPP
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
#include "fabulous/algo/BGMRes.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/utils/Logger.hpp"
namespace fabulous {
/**
* \brief main fabulous solver entry point
*/
template<class Algo, class Equation, class Restart,
class S = typename Equation::value_type,
class P = typename Equation::primary_type >
inline Logger<P> solve(Equation &eq, const Algo &algo, Parameters param,
Orthogonalizer ortho, Restart restart)
{
BGMRes<S> solver;
solver.solve(eq, algo, param, ortho, restart);
return solver.get_logger();
}
/**
* \brief main fabulous solver entry point
*/
template<class Algo, class Equation,
class S = typename Equation::value_type,
class P = typename Equation::primary_type >
inline Logger<P> solve(Equation &eq, const Algo &algo, Parameters param,
Orthogonalizer ortho = orthogonalization())
{
BGMRes<S> solver;
solver.solve(eq, algo, param, ortho, classic_restart());
return solver.get_logger();
}
};
#endif // FABULOUS_SOLVE_HPP
......@@ -152,8 +152,8 @@ private:
}
public:
template<template <class> class Restarter>
ArnoldiDR(Logger<P> &logger, Restarter<S> &restarter,
template<class RestartParam>
ArnoldiDR(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_hess{max_krylov_space_size, nbRHS},
......@@ -167,7 +167,7 @@ public:
_cold_restart{false}
{
static_assert(
hessenbergXrestarter<HESSENBERG, Restarter>::value,
hessenbergXrestarter<HESSENBERG, RestartParam>::value,
"The Hessenberg used is not compatible with the kind of restarting selected"
);
}
......@@ -184,8 +184,7 @@ public:
* \param max_krylov_space_size Maximum size of Krylov Search space.
* \param[in] epsilon tolerance (backward error) for residual
* \param[in,out] restarter hold data to deflate at restart
* \param ortho_scheme orthogonalization schema
* \param ortho_type Choose between block wise arnoldi (through QR factorization)
* \param ortho Orthogonalizer
* or vector wise arnoldi. (In distributed, only the vector wise will works)
*
* \return whether the convergence was reached
......@@ -195,8 +194,7 @@ public:
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter &restarter,
OrthoScheme ortho_scheme,
OrthoType ortho_type )
Orthogonalizer &ortho )
{
print_header(max_krylov_space_size);
......@@ -229,7 +227,8 @@ public:
Block W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
Arnoldi_Ortho(_hess, _base, W, ortho_type, ortho_scheme, A);
ortho.run(_hess, _base, W, A);
_base.increase(_nbRHS); // add W to base
_size_krylov_space = _base.get_nb_vect();
//_hess.display_hess_bitmap("Hess Before Least Square");
......
......@@ -327,8 +327,8 @@ private:
public:
template<template <class> class Restarter>
ArnoldiIB(Logger<P> &logger, Restarter<S> &restarter,
template<class RestartParam>
ArnoldiIB(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
L{max_krylov_space_size+1, nbRHS},
......@@ -346,7 +346,7 @@ public:
_cold_restart{false}
{
static_assert(
hessenbergXrestarter<HESSENBERG, Restarter>::value,
hessenbergXrestarter<HESSENBERG, RestartParam>::value,
"The Hessenberg used is not compatible with the kind of restarting selected"
);
}
......@@ -362,8 +362,7 @@ public:
* \param[in] B the right hand sides
* \param max_krylov_space_size maximum size for the Krylov search space
* \param[in] epsilon tolerance for Residual
* \param ortho_scheme orthogonalization schema
* \param ortho_type Choose between BLOCK-wise arnoldi (through QR factorization)
* \param ortho Orthogonalizer
* or vector wise arnoldi(RUHE). (In distributed, only the vector wise will works)
*
* \return whether the convergence was reached
......@@ -373,8 +372,7 @@ public:
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter&,
OrthoScheme ortho_scheme,
OrthoType ortho_type )
Orthogonalizer &ortho)
{
bool convergence = false;
print_header(max_krylov_space_size);
......