Commit 1ccdf08a authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Make IB look like IB+DR

parent 7898fe9f
......@@ -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";
}
......
This diff is collapsed.
......@@ -38,7 +38,7 @@ private:
void print_header(int nb_mvp, std::ostream &o = std::cout)
{
o << "#######################################################\n";
o << "#################### Arnoldi IB DR ####################\n";
o << "#################### Arnoldi ##########################\n";
o << "######### Mat Vect product scheduled: "<< nb_mvp <<" ###########\n";
}
......@@ -188,7 +188,7 @@ private:
}
template<class Matrix, class Block, class Epsilon, class Restarter>
void initial_iteration(Matrix &A, Block &X, Block &B, Restarter &restarter,
bool initial_iteration(Matrix &A, Block &X, Block &B, Restarter &restarter,
const Epsilon &epsilon, const Epsilon &inv_epsilon)
{
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
......@@ -207,6 +207,7 @@ private:
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
return std::get<2>(MinMaxConv);
}
/************************* IB on R0 ********************************/
......@@ -253,15 +254,15 @@ private:
_F.init_phi(_nb_direction_kept); // Create and initialize Phi (\Phi_1)
}
template< class Matrix, class Block, class Epsilon >
template< class Matrix, class Epsilon >
std::tuple<P,P,bool>
inexact_breakdown_on_R0( Matrix &A, Block &X, Block &B,
inexact_breakdown_on_R0( Matrix &A, Block<S> &X, Block<S> &B,
const Epsilon &epsilon, const Epsilon &inv_epsilon)
{
Block Q0{_dim, _nbRHS};
Block Lambda0{_nbRHS, _nbRHS};
Block L0{_nbRHS, _nbRHS};
Block U{_nbRHS, _nbRHS};
Block<S> Q0{_dim, _nbRHS};
Block<S> Lambda0{_nbRHS, _nbRHS};
Block<S> L0{_nbRHS, _nbRHS};
Block<S> U{_nbRHS, _nbRHS};
compute_real_residual(A, X, B, Q0);
A.QRFacto(Q0, Lambda0); L0.copy(Lambda0);
......@@ -281,7 +282,7 @@ private:
handle_IB_on_R0(U, Q0, Lambda0);
} else {
std::cout<<"\t\tNo INEXACT BREAKDOWN on R0\n";
Block W = _base.get_W(_nbRHS);
Block<S> W = _base.get_W(_nbRHS);
W.copy(Q0);
_base.increase(_nbRHS);
_F.init_lambda(Lambda0);
......@@ -299,15 +300,14 @@ private:
/**
* \brief perform update on WP and Hessenberg due to Inexact Breakdown
*
* \param max_krylov_space_size maximum space (nb_vect) of krylov space base
* \return whether max_krylov_space_size was not reached when adding new vector into base
*/
void do_IB_update()
{
if (!_ib_update_scheduled) {
return;
}
_ib_update_scheduled = false;
Block<S> &U = _U;
assert( U.get_nb_col() == _nb_direction_kept );
......@@ -339,7 +339,6 @@ private:
// L_{j+1,} = | \W_1^{H} | * H_j
// Gj = | \W_2^{H} | * H_j
_F.update_bottom_line(W1_W2);
_ib_update_scheduled = false;
}
template<class P>
......@@ -372,14 +371,7 @@ private:
);
}
}
//if (_nb_direction_discarded > 0) {
//check_arnoldi_formula_IB(A, _base, _F, _WP);
_ib_update_scheduled = true;
/*} else {
Block<S> W = _base.get_W(_nbRHS);
W.copy(_WP);
_base.increase(_nbRHS);
}*/
}
/*********************************************************/
......@@ -387,7 +379,7 @@ private:
public:
template<class RestartParam>
ArnoldiIBDR(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_F{max_krylov_space_size+1, nbRHS},
_solution_computed{false},
......@@ -433,19 +425,15 @@ public:
const std::vector<P> &epsilon,
Restarter &restarter, OrthoParam &orthoparm)
{
bool convergence = false;
print_header(max_krylov_space_size);
assert( _nbRHS == B.get_nb_col() ); // Initial Size of block
assert( _dim == B.get_nb_row() );
// IS_INVARIANT( _nb_direction_discarded + _nb_direction_kept == nbRHS );
const auto inv_epsilon = array_inverse(epsilon);
initial_iteration(A, X, B, restarter, epsilon, inv_epsilon);
// Should not be true (because we already tested it):
convergence = (_nb_direction_discarded == _nbRHS);
bool convergence = initial_iteration(A, X, B, restarter, epsilon, inv_epsilon);
while (!convergence && // Main Loop
while (!convergence && // Main Loop
_size_krylov_space+_nb_direction_kept < max_krylov_space_size)
{
++ _iteration_count;
......
......@@ -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;
}
......
......@@ -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