Commit 4108a3aa authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Merge branch 'ibdr_improvements' into develop

parents fc2a1b92 1ccdf08a
......@@ -124,9 +124,9 @@ ggplot(df, aes(x=nb_mvp)) +
cd ${WORKDIR}/build/src/test_core/
mkdir -p ../data/res
#nbRHS=6; maxSpace=90
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -u -o "BGMRES"
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -A QR -u -o "BGMRES"
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -A IB -u -o "IB-BGMRES"
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -r DEFLATED -p 5 -u -o "BGMRES-DR"
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -A QRDR -r DEFLATED -p 5 -u -o "BGMRES-DR"
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -A IBDR -r DEFLATED -p 5 -u -o "IB-BGMRES-DR"
#+end_src
*** plot the graphic
......
......@@ -43,27 +43,27 @@ inline Template<ARNOLDI_CHAM_TL> arnoldi_cham_tl()
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_STD> &)
{
return (o << "STD");
return (o << "");
}
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_QR> &)
{
return (o << "QR");
return (o << "QR-");
}
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_QRDR> &)
{
return (o << "QR");
return (o << "QR-");
}
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_IB> &)
{
return (o << "IB");
return (o << "IB-");
}
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_IBDR> &)
{
return (o << "IBDR");
return (o << "IB-");
}
#ifdef FABULOUS_USE_CHAMELEON
......
......@@ -65,16 +65,19 @@ private:
o <<"######################################################\n";
}
template< class Matrix, class RestartParam >
template< class Matrix, class RestartParam, class Algo >
void print_start_info(
int dim, int nbRHS, int maxMVP, int max_krylov_space_size,
const OrthoParam &ortho, const RestartParam &restart_param,
const std::vector<P> &epsilon, Matrix &A, std::ostream &o = std::cout)
const std::vector<P> &epsilon, Matrix &A, Algo &algo, std::ostream &o = std::cout)
{
Arithmetik<S> arith;
o << "#######################################################\n";
o << "#################### BGMRES #########################\n";
std::stringstream ss;
ss << algo << "BGMRES" << get_type(restart_param);
int l = (55-(ss.str().length()+2))/2;
o << std::string(l, '#') << " " << ss.str() << " " << std::string(l, '#') << "\n";
o << "#######################################################\n";
o << "Dimension of Problem : "<<dim<<"\n";
o << "Number of Right Hand Side : "<<nbRHS<<"\n";
......@@ -164,7 +167,7 @@ public:
*/
template< class Algo, class RestartParam, class Matrix >
int solve( Matrix &A, Block<S> &X, Block<S> &B,
Algo&, const int max_mvp, const int max_krylov_space_size,
Algo &algo, const int max_mvp, const int max_krylov_space_size,
const std::vector<P> &epsilon,
OrthoParam ortho, RestartParam restart_param )
{
......@@ -173,7 +176,7 @@ public:
const int dim = B.get_nb_row();
reset();
print_start_info( dim, nbRHS, max_mvp, max_krylov_space_size,
ortho, restart_param, epsilon, A );
ortho, restart_param, epsilon, A, algo );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
inv_normB = array_inverse(normB);
......
......@@ -36,11 +36,10 @@ public:
using P = primary_type;
private:
void print_header(int nb_mvp, std::ostream &o = std::cout)
{
o << "#######################################################\n";
o << "################# Arnoldi Standard ###################\n";
o << "#################### Arnoldi ##########################\n";
o << "######## Mat Vect product scheduled "<< nb_mvp <<"###########\n";
}
......@@ -124,6 +123,15 @@ private:
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
}
template<class Matrix>
void do_ortho(Matrix &A, OrthoParam &orthoparm, Block<S> &W)
{
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A); // Ortho and filling of L (the (IB) Hessenberg)
_logger.notify_ortho_end();
}
template<class Matrix, class Block>
bool check_real_residual_if(bool projected_residual_converged,
Matrix &A, const Block &B, Block &X,
......@@ -140,7 +148,6 @@ private:
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
auto MinMaxConv = R.check_precision(epsilon, A);
#ifdef FABULOUS_DEBUG_MODE
auto min = std::get<0>(MinMaxConv);
auto max = std::get<1>(MinMaxConv);
......@@ -156,6 +163,32 @@ private:
return false;
}
template<class Matrix, class Epsilon, class Restarter>
bool initial_iteration(Matrix &A, Block<S> &X, const Block<S> &B,
const Epsilon &epsilon, Restarter &restarter)
{
std::tuple<P,P,bool> MinMaxConv;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
MinMaxConv = restarter.run(_hess, epsilon);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
Block<S> R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
MinMaxConv = R1.check_precision(epsilon);
}
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
return std::get<2>(MinMaxConv);
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
......@@ -210,28 +243,11 @@ public:
{
print_header(max_krylov_space_size);
std::tuple<P,P,bool> MinMaxConv;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
MinMaxConv = restarter.run(_hess, epsilon);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
Block<S> R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
MinMaxConv = R1.check_precision(epsilon);
}
bool convergence = std::get<2>(MinMaxConv);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
bool convergence = initial_iteration(A, X, B, epsilon, restarter);
while (!convergence && _size_krylov_space < max_krylov_space_size) {
while (!convergence && // Main Loop
_size_krylov_space + _nbRHS < max_krylov_space_size)
{
++ _iteration_count;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
Block<S> Vj = _base.get_Vj(); // V_j
......@@ -239,22 +255,16 @@ public:
Block<S> W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A);
_logger.notify_ortho_end();
do_ortho(A, orthoparm, W);
_base.increase(_nbRHS); // add W to base
_size_krylov_space = _base.get_nb_vect();
//_hess.display_hess_bitmap("Hess Before Least Square");
//_hess.display_hess_bitmap("Hess Before Least Square");
auto MinMaxConv = check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, B, X, epsilon);
//check_arnoldi_formula(A, _base, _hess, _nbRHS);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, _hess);
......@@ -269,10 +279,9 @@ public:
FABULOUS_WARNING("cold restart");
_base.reset();
} else if (restarter && !convergence) {
int M = _hess.get_nb_vect() + _nbRHS;
Block<S> LS{M, _nbRHS};
_hess.compute_proj_residual(LS);
restarter.save_LS(LS);
Block<S> PR{_hess.get_nb_vect() + _nbRHS, _nbRHS};
_hess.compute_proj_residual(PR);
restarter.save_LS(PR);
Block<S> H = _hess.get_hess_block();
restarter.save_hess(H);
}
......
This diff is collapsed.
This diff is collapsed.
......@@ -75,7 +75,7 @@ inline DeflatedRestart<S> deflated_restart(int k, S target)
*/
std::ostream& operator<<(std::ostream &o, const Type<ClassicRestart>&)
{
o << "Classic_Restarting";
o << "";
return o;
}
......@@ -85,7 +85,7 @@ std::ostream& operator<<(std::ostream &o, const Type<ClassicRestart>&)
template<class S>
std::ostream& operator<<(std::ostream &o, const Type<DeflatedRestart<S>>&)
{
o << "Deflated_Restarting";
o << "-DR";
return o;
}
......@@ -265,9 +265,9 @@ public:
operator bool () { return true; }
Base<S> &get_base() { return _base; }
template< class BLOCK > void save_LS(BLOCK &LS) { _LS = LS.copy(); }
template< class BLOCK > void save_hess(BLOCK &hess) { _saved_Hess = hess.copy(); }
template< class BLOCK > void save_WP(BLOCK &WP) { _saved_WP = WP.copy(); }
template< class BLOCK > void save_LS(BLOCK &LS) { _LS = LS; }
template< class BLOCK > void save_hess(BLOCK &hess) { _saved_Hess = hess; }
template< class BLOCK > void save_WP(BLOCK &WP) { _saved_WP = WP; }
/**
* \brief This function handle the deflation of eigenvalues at restart. (STD variant)
......
......@@ -3,6 +3,7 @@
#include <algorithm>
#include <vector>
#include <iostream>
#include <string>
#include <sstream>
#include <complex>
......
......@@ -49,7 +49,7 @@ public:
void QRFacto(Block &Q, Block &R) const
{
// default implementation: ( distributed QR; uses DotProduct )
// ::fabulous::Algorithm::InPlaceQRFactoMGS_User(Q, R, *this);
//::fabulous::Algorithm::InPlaceQRFactoMGS_User(Q, R, *this);
// default implementation: (not distributed)
::fabulous::Algorithm::InPlaceQRFacto(Q, R);
......
......@@ -136,9 +136,9 @@ TestRet<P> run_test_BGMRES_filelog(
name << filename;
if (append_suffix) {
name <<"_" <<ortho.get_scheme()<<"_"<<ortho.get_type()<<"_"
<< algo <<"_"<< get_type(restart);
<< algo <<"BGMRES"<< get_type(restart);
}
filepath << "../data/res/"<< name.str() << ".res";
filepath << "../data/res/" << name.str() << ".res";
std::ofstream file{filepath.str(), std::ios::out};
if (file.is_open()) {
......@@ -158,10 +158,10 @@ static void print_elapsed_time(const TestRet<P> &r, const std::string &test_name
{
std::cout << test_name<<":\n"
<< "\tElapsed time:\t"<<r.elapsed << " seconds for "
<<r.nb_mvp<<"/"<<r.maxMVP<<" iterations\n"
<< r.nb_mvp<<"/"<<r.maxMVP<<" iterations\n"
<< "\tForward error:\n"
<<"\t\tMin:\t"<<r.minmax.first<<"\n"
<<"\t\tMax:\t"<<r.minmax.second<<"\n";
<< "\t\tMin:\t"<<r.minmax.first<<"\n"
<< "\t\tMax:\t"<<r.minmax.second<<"\n";
}
}; // end namespace fabulous
......
Markdown is supported
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