Commit b98831fb authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Add timer to measure independently ortho, MVP and LeastSquare

parent 4b9adb93
......@@ -7,6 +7,7 @@
* DONE Implement DR+QR
* HOLD Reproducible installation
** DONE Instruction in INSTALL.org
** DONE Update spack package
** HOLD merge request to make Chameleon headers compatible with C++ complex types
* TODO reproducible results and visualization
** TODO Improve Logging
......@@ -17,7 +18,7 @@
*** TODO Add timer for Least square
*** TODO Add timer for Matrix Vector Product
** TODO Improve RESULTS.org
Eventually, anyone must be able to gather into the file all interesting results
just by evaluating code blocks and/or tangling file from RESULTS.org
Eventually, anyone must be able to gather all interesting results into RESULTS.org
just by evaluating code blocks from RESULTS.org and/or tangling RESULTS.org
* TODO Implement IB+DR
* TODO Implement IB+DR+QR
......@@ -95,13 +95,16 @@ private:
void compute_real_residual(Matrix &A, const Block &X, const Block &B, Block &R)
{
R.copy(B); // R <- B
_logger.notify_mvp_begin();
A.MatBlockVect(X, R, S{-1.0}, S{1.0}); // R <- B-A*X
_logger.notify_mvp_end();
_nb_mvp += X.get_nb_col();
}
template< class Matrix, class Block >
void W_AMV(Matrix &A, Block &V, Block &W) // W = A*M*V
{
_logger.notify_mvp_begin();
if (A.useRightPrecond()) {
Block MV{W.get_nb_row(), V.get_nb_col()};
A.PrecondBlockVect(V, MV); // MV = M * V
......@@ -109,6 +112,7 @@ private:
} else {
A.MatBlockVect(V, W); // W = A * V
}
_logger.notify_mvp_end();
_nb_mvp += V.get_nb_col();
}
......@@ -152,6 +156,15 @@ private:
return false;
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
_logger.notify_least_square_check_begin();
auto minmaxconv = _hess.check_least_square_residual(epsilon);
_logger.notify_least_square_check_end();
return minmaxconv;
}
public:
template<class RestartParam>
ArnoldiDR(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
......@@ -190,8 +203,8 @@ public:
*
* \return whether the convergence was reached
*/
template<class Matrix, class Block, class Restarter>
bool run( const Matrix &A, Block &X, const Block &B,
template<class Matrix, class Restarter>
bool run( const Matrix &A, Block<S> &X, const Block<S> &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter &restarter,
......@@ -199,14 +212,14 @@ public:
{
print_header(max_krylov_space_size);
Block R1{}; // R1 = RHS of projected problem
Block<S> R1{}; // R1 = RHS of projected problem
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
restarter.run(R1, _hess);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block R0 = _base.get_W(_nbRHS);
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
R1.realloc(_nbRHS, _nbRHS);
......@@ -223,18 +236,20 @@ public:
while (!convergence && _size_krylov_space < max_krylov_space_size) {
++ _iteration_count;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
Block Vj = _base.get_Vj(); // V_j
Block<S> Vj = _base.get_Vj(); // V_j
assert( Vj.get_nb_col() == _nbRHS );
Block W = _base.get_W(_nbRHS); // W (candidate)
Block<S> W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A);
_logger.notify_ortho_end();
_base.increase(_nbRHS); // add W to base
_size_krylov_space = _base.get_nb_vect();
//_hess.display_hess_bitmap("Hess Before Least Square");
auto MinMaxConv = _hess.check_least_square_residual(epsilon);
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);
_logger.notify_iteration_end(
......@@ -251,10 +266,10 @@ public:
_base.reset();
} else if (restarter && !convergence) {
int M = _hess.get_nb_vect() + _nbRHS;
Block LS{M, _nbRHS};
Block<S> LS{M, _nbRHS};
_hess.compute_proj_residual(LS);
restarter.save_LS(LS);
Block H = _hess.get_hess_block();
Block<S> H = _hess.get_hess_block();
restarter.save_hess(H);
}
......
......@@ -49,7 +49,7 @@ private:
private:
Logger<P> &_logger;
HESSENBERG<S> L;
HESSENBERG<S> _F;
bool _solution_computed;
Block<S> _solution;
Block<S> _last_Y;
......@@ -77,7 +77,6 @@ public:
int get_nb_mvp() const { return _nb_mvp; }
int get_nb_iteration() const { return _iteration_count; }
private:
template< class Matrix, class Block, class Base >
void X_plus_MVY(Matrix &A, Block &X, Base &V, Block &Y)
......@@ -99,11 +98,13 @@ private:
}
template<class Matrix, class Block>
void compute_solution(Matrix &A, Block &X)
void compute_real_residual(Matrix &A, Block &X, const Block &B, Block &R)
{
Block Y = L.alloc_least_square_sol();
L.solve_least_square(Y);
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
R.copy(B); // R <- B
_logger.notify_mvp_begin();
A.MatBlockVect(X, R, S{-1.0}, S{1.0}); // R <- B-A*X
_logger.notify_mvp_end();
_nb_mvp += X.get_nb_col();
}
template< class Matrix, class Block, class BlockWP >
......@@ -111,6 +112,7 @@ private:
{
Block W = WP.get_W();
_logger.notify_mvp_begin();
if (A.useRightPrecond()) {
Block MV{W.get_nb_row(), V.get_nb_col()};
A.PrecondBlockVect(V, MV); // MV = M * V
......@@ -118,15 +120,16 @@ private:
} else {
A.MatBlockVect(V, W); // W = A * V
}
_logger.notify_mvp_end();
_nb_mvp += V.get_nb_col();
}
template<class Matrix, class Block>
void compute_real_residual(Matrix &A, Block &X, const Block &B, Block &R)
void compute_solution(Matrix &A, Block &X)
{
R.copy(B); // R <- B
A.MatBlockVect(X, R, S{-1.0}, S{1.0}); // R <- A*X-B
_nb_mvp += X.get_nb_col();
Block Y = _F.alloc_least_square_sol();
_F.solve_least_square(Y);
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
}
template<class Matrix, class Block>
......@@ -182,7 +185,7 @@ private:
W2 = W1_W2.sub_block(0, _nb_direction_kept, _nbRHS, _nb_direction_discarded);
U1_2.OutOfPlaceQRFacto(W1_W2, R); // Qr facto of bottom block of U
L.update_phi(W1_W2, _nb_direction_kept);
_F.update_phi(W1_W2, _nb_direction_kept);
// Update Phi with \W_1 and \W_2, doc Annexe Eq 2.7
Block V = _base.get_W(_nb_direction_kept);
......@@ -192,7 +195,7 @@ private:
// L_{j+1,} = | \W_1^{H} | * H_j
// Gj = | \W_2^{H} | * H_j
L.update_bottom_line(W1_W2);
_F.update_bottom_line(W1_W2);
}
void handle_IB_on_R0( Block<S> &U, BlockWP<S> &WP,
......@@ -235,17 +238,16 @@ private:
Lambda1.get_ptr(), Lambda1.get_leading_dim(),
S{1.0}, S{0.0}
);
L.init_lambda(Lambda1);
L.init_phi(_nb_direction_kept); // Create and initialize Phi (\Phi_1)
_F.init_lambda(Lambda1);
_F.init_phi(_nb_direction_kept); // Create and initialize Phi (\Phi_1)
}
template< class Matrix, class Block, class BlockWP,
class Epsilon, class Hess >
template< class Matrix, class Block, class BlockWP, class Epsilon >
std::tuple<P,P,bool>
inexact_breakdown_on_R0(
const Matrix &A, Block &R0,
const Epsilon &epsilon, const Epsilon &inv_epsilon,
Hess &L, BlockWP &WP)
BlockWP &WP)
{
Block Q0 = R0;
Block Lambda0{_nbRHS, _nbRHS};
......@@ -259,7 +261,7 @@ private:
_old_nb_direction_kept = _nb_direction_kept;
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
// Test if needed to launch algorithm again
// Test if algorithm need to be launched again
if (_nb_direction_kept == 0) {
_IB_happened = true;
std::cout << "\t\tTOTAL INEXACT BREAKDOWN on R0\n";
......@@ -277,7 +279,7 @@ private:
Block W = _base.get_W(_nbRHS);
W.copy(Q0);
_base.increase(_nbRHS);
L.init_lambda(Lambda0);
_F.init_lambda(Lambda0);
}
_block_size_sum.clear();
_block_size_sum.push_back(_nb_direction_kept);
......@@ -291,10 +293,10 @@ private:
const std::vector<P> &epsilon, const std::vector<P> &inv_epsilon,
const Matrix &A, const Block<S> &X)
{
Block<S> Y = L.alloc_least_square_sol();
Block<S> PR{L.get_nb_hess_line(), _nbRHS};
L.solve_least_square(Y);
L.compute_proj_residual(PR);
Block<S> Y = _F.alloc_least_square_sol();
Block<S> PR{_F.get_nb_hess_line(), _nbRHS};
_F.solve_least_square(Y);
_F.compute_proj_residual(PR);
Block<S> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
......@@ -303,8 +305,8 @@ private:
if (_nb_direction_kept == 0) { // Suspicious ...
if (!proj_convergence) { // Even more suspicious
FABULOUS_FATAL_ERROR(
"Directions kepts during SVD of residuals is equal to 0, "
"but residual itself has not converged. "
"Directions kepts during SVD of projected residual is "
" equal to 0, but residual itself has not converged. "
"That should not be possible."
);
} else {
......@@ -335,13 +337,21 @@ private:
_size_krylov_space = _base.get_nb_vect();
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
_logger.notify_least_square_check_begin();
auto minmaxconv = _F.check_least_square_residual(epsilon);
_logger.notify_least_square_check_end();
return minmaxconv;
}
public:
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},
_F{max_krylov_space_size+1, nbRHS},
_solution_computed{false},
_base{restarter.get_base()},
_dim(dim),
......@@ -377,8 +387,8 @@ public:
*
* \return whether the convergence was reached
*/
template< class Matrix, class Block, class Restarter >
bool run( Matrix &A, Block &X, Block &B,
template< class Matrix, class Restarter >
bool run( Matrix &A, Block<S> &X, Block<S> &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter&,
......@@ -397,10 +407,10 @@ public:
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
// STEP1: Compute first block residual
Block R0{_dim, _nbRHS};
Block<S> R0{_dim, _nbRHS};
compute_real_residual(A, X, B, R0);
// auto minmaxR0 = R0.get_min_max_norm(A);
auto MinMaxConv = inexact_breakdown_on_R0(A, R0, epsilon, inv_epsilon, L, WP);
auto MinMaxConv = inexact_breakdown_on_R0(A, R0, epsilon, inv_epsilon, WP);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
......@@ -416,21 +426,23 @@ public:
_size_krylov_space < max_krylov_space_size) // Main Loop
{
++ _iteration_count;
if (! L.check_room( WP.get_size_W()) )
if (! _F.check_room( WP.get_size_W()) )
break;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
Block Vj = _base.get_Vj();
Block<S> Vj = _base.get_Vj();
W_AMV(A, Vj, WP); // W <- A*M*Vj: perform matrix multiplication
// Ortho and filling of L (the (IB) Hessenberg)
ortho.run(L, _base, WP, A);
_logger.notify_ortho_begin();
ortho.run(_F, _base, WP, A);
_logger.notify_ortho_end();
// WP.check_ortho_WP();// WP.get_W().check_ortho();
//L.display_hess_extended_bitmap();
//_F.display_hess_extended_bitmap();
// Y is the solution of LeastSquare problem: Y = argmin(||F*Y-\Lambda||)
// PR is the projected residual PR <- F*Y-\Lambda
auto MinMaxConv = L.check_least_square_residual(epsilon);
auto MinMaxConv = check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, X, B, epsilon);
......@@ -441,7 +453,7 @@ public:
_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, L);
_logger.log_real_residual(A, B, X, _base, _F);
}
if (_solution_computed)
......@@ -464,7 +476,6 @@ public:
}
};
}; // namespace fabulous
#endif // FABULOUS_ARNOLDI_IB_HPP
......@@ -179,7 +179,9 @@ public:
template<class HESS >
int run(Block<S> &R1, HESS &hess)
{
double tic = Timer::get_time();
Timer time;
time.start();
int p = _LS.get_nb_col();
int nm = _base.get_nb_vect() - p;
int dim = _base.get_nb_row();
......@@ -315,11 +317,10 @@ public:
hess.set_rhs(R1);// important to do this before orthogonalization end
hess.notify_orthogonalization_end();
double tac = Timer::get_time();
std::cout<<"\tElapsed Time for Computing the Deflation at restart :: "
<<tac-tic<<" seconds\n";
std::cout<<Color::bold<<Color::yellow;
std::cout<<"\t\tDR done\n"<<Color::reset;
time.stop();
double length = time.get_length();
std::cout<<"\tElapsed Time for Computing the Deflation at restart: "<<length<<" seconds\n";
std::cout<<Color::bold<<Color::yellow<<"\t\tDR done\n"<<Color::reset;
return newK;
}
};
......
......@@ -16,7 +16,8 @@ namespace fabulous {
* \brief Struct containing information about one iterations
*/
template<class P>
struct LogEntry {
struct LogEntry
{
int glob_iteration; /*!< global iteration taking restart into account */
int loc_iteration; /*!< current iteration in arnoldi procedure */
int krylov_space_size; /*!< Index of Current iteration*/
......@@ -26,13 +27,16 @@ struct LogEntry {
P max; /*!< Arnoldi Maximum norm reached at this iteration*/
P minReal; /*!< Real Minimal norm reached at this iteration*/
P maxReal; /*!< Real Maximum norm reached at this iteration*/
double timestamp; /*!< Time spent for this iteration */
double time_spent; /*!< Time spent for this iteration */
double least_square_spent; /*!< Time spent for this iteration in least_square */
double mvp_spent; /*!< Time spent for this iteration in mvp */
double ortho_spent; /*!< Time spent for this iteration in orthogonalization */
static void print_header(std::ostream &o = std::cout)
{
o <<"global_iteration\tlocal_iteration\tkrylov_space_size\tnb_mvp\t"
"current_block_size\tminRes\tmaxRes\tminRealRes\tmaxRealRes\ttime\tID\n";
"current_block_size\tminRes\tmaxRes\tminRealRes\tmaxRealRes\ttime\t"
"least_square_time\tmvp_spent\tortho_spent\tID\n";
}
void print(std::string ID, std::ostream &o = std::cout)
......@@ -46,7 +50,10 @@ struct LogEntry {
<< max<<"\t"
<< minReal<<"\t"
<< maxReal<<"\t"
<< timestamp<<"\t"
<< time_spent<<"\t"
<< least_square_spent<<"\t"
<< mvp_spent<<"\t"
<< ortho_spent<<"\t"
<< ID << "\n";
}
};
......@@ -55,24 +62,36 @@ struct LogEntry {
* \brief Class that will store informations during the solve at each
* iteration.
*/
template<class P, class Entry=LogEntry<P> >
class Logger {
Timer _timer;
template<class P>
class Logger
{
using Entry = LogEntry<P>;
bool _started;
double _total_time;
Timer _global_timer;
Timer _iterations_timer;
Timer _least_square_timer;
Timer _mvp_timer;
Timer _ortho_timer;
bool _log_real_residual;
int _global_iteration;
int _total_mvp;
int _local_mvp;
int _last_mvp;
std::vector<Entry> _iterations;
std::vector<double> _array;
double _last_timestamp;
public:
Logger(bool log_real_residual = false):
_started{false},
_total_time{0.0},
_log_real_residual{log_real_residual},
_global_iteration{0},
_total_mvp{0},
_local_mvp{0}
_last_mvp{0}
{
}
......@@ -81,51 +100,109 @@ public:
void reset()
{
_started = false;
_global_timer.reset();
_iterations_timer.reset();
_least_square_timer.reset();
_mvp_timer.reset();
_ortho_timer.reset();
_iterations.clear();
_array.clear();
_global_iteration = 0;
_total_mvp = 0;
_local_mvp = 0;
_last_mvp = 0;
}
const Entry &get_entry(int i) const { return _iterations.at(i); }
const std::vector<Entry> &get_iterations() const { return _iterations; }
double get_total_time() const { return _total_time; }
/* -------------------- notify BEGIN --------------------------------------*/
void notify_iteration_begin(int iteration_count, int size_krylov_space)
{
_last_timestamp = _timer.get_time();
_iterations_timer.reset();
_least_square_timer.reset();
_mvp_timer.reset();
_ortho_timer.reset();
std::cout<<"\nStart of Iteration ["
<<Color::note<<iteration_count<<" :: "
<<size_krylov_space<<Color::reset<<"]\n";
if (iteration_count == 0)
_local_mvp = 0;
if (iteration_count == 0) {
_last_mvp = 0;
if (!_started) {
_global_timer.start();
_started = true;
}
}
_iterations_timer.start();
}
void notify_mvp_begin()
{
_mvp_timer.start();
}
void notify_ortho_begin()
{
_ortho_timer.start();
}
void notify_least_square_check_begin()
{
_least_square_timer.start();
}
/* -------------------- notify END --------------------------------------*/
void notify_mvp_end()
{
_mvp_timer.stop();
}
void notify_ortho_end()
{
_ortho_timer.stop();
}
void notify_least_square_check_end()
{
_least_square_timer.stop();
}
/**
* \brief Add an iteration to the list.
*/
void notify_iteration_end(int id, int size_krylov_space, int local_mvp,
void notify_iteration_end(int id, int size_krylov_space, int mvp,
std::tuple<P, P, bool> min_max_conv)
{
_iterations_timer.stop();
double elapsed = _iterations_timer.get_length();
double least_square_elapsed = _least_square_timer.get_length();
double mvp_elapsed = _mvp_timer.get_length();
double ortho_elapsed = _ortho_timer.get_length();
_total_time = _global_timer.get_length();
auto min = std::get<0>(min_max_conv);
auto max = std::get<1>(min_max_conv);
int mvp = local_mvp - _local_mvp;
_local_mvp = local_mvp;
_total_mvp += mvp;
Entry entry = Entry{
_global_iteration+1, id, size_krylov_space,
_total_mvp, mvp, min, max, 0, 0,
_timer.get_time() - _last_timestamp
};
_iterations.push_back(entry);
int mvp_diff = mvp - _last_mvp; _last_mvp = mvp;
_total_mvp += mvp_diff;
++ _global_iteration;
_iterations.push_back(
Entry{ _global_iteration, id, size_krylov_space,
_total_mvp, mvp_diff, min, max, 0, 0,
elapsed, least_square_elapsed, mvp_elapsed, ortho_elapsed }
);
FABULOUS_NOTE("\t\t "<<Color::red<<"MIN="<<Color::reset<<std::scientific<<min
<<"\t"<<Color::red<<"MAX="<<Color::reset<<std::scientific<<max);
if (std::isnan(min) || std::isnan(max))
FABULOUS_FATAL_ERROR("nan encoutered");
std::cout<<"End of Iteration ["<<Color::note<<id<<" :: "
std::cout<< "End of Iteration ["<<Color::note<<id<<" :: "
<< size_krylov_space<<Color::reset<<"]\n\n";
}
......@@ -168,7 +245,6 @@ public:
_iterations.back().maxReal = pair.second;
}
/**
* \brief iterator over the iterations with everything.
*/
......@@ -180,7 +256,7 @@ public:
it.current_block_size,
it.min, it.max,
it.minReal, it.maxReal,
it.timestamp);
it.time_spent);
}
}
......@@ -195,7 +271,7 @@ public:
it.current_block_size,
it.min, it.max,
it.minReal, it.maxReal,
it.timestamp,
it.time_spent,
it.loc_iteration,
it.krylov_space_size);
}
......@@ -212,7 +288,7 @@ public:
it.current_block_size,
it.min,
it.max,
it.timestamp);
it.time_spent);
}
}
......@@ -227,7 +303,7 @@ public:
for (int i = 0; i < size; ++i) {
_array[3*i + 0] = _iterations[i].min;
_array[3*i + 1] = _iterations[i].max;
_array[3*i + 2] = _iterations[i].timestamp;
_array[3*i + 2] = _iterations[i].time_spent;
}
return _array.data();
}
......
......@@ -3,11 +3,12 @@
#include <time.h>
#include <sys/time.h>
#include "fabulous/utils/Error.hpp"
namespace fabulous {
/**
* \brief time measurements facilities
* \brief time measurements facilities (chronometer)
*/
struct Timer {
public:
......@@ -15,12 +16,55 @@ public:
/**
* \brief get a time stamp.
*/
static double get_time()
static double get_timestamp()
{
timeval t;
gettimeofday(&t, NULL);
return double(t.tv_sec) + (double(t.tv_usec)/1000000.0);
}
private:
bool _running;
double _last_start;
double _length;
public:
Timer() { reset(); }
void reset()
{
_running = false;
_last_start = -1.0;
_length = 0.0;
}
void start()
{
if (_running) {
FABULOUS_FATAL_ERROR("this timer is already started");
}
_last_start = get_timestamp();
_running = true;
}
void stop()
{
if (!_running) {
FABULOUS_FATAL_ERROR("this timer is not running");
}
double tsp = get_timestamp();
_length += (tsp - _last_start);
_running = false;
}
double get_length()
{
if (_running) {
double tsp = get_timestamp();
return _length + (tsp - _last_start);
}
return _length;
}
};
}; // namespace fabulous
......
......@@ -121,12 +121,14 @@ TestRet<P> run_test_BGMRES_filelog(
const int dim = RHS.get_nb_row();
Block X{dim, nbRHS};
double tic = Timer::get_time();
Timer time;
time.start();
auto eq = equation(mat, X, RHS, epsilon);
auto param = parameters(maxMVP, maxKrylovSpaceSize);
auto logger = solve(eq, algo, param, ortho, restart);
double elapsed = Timer::get_time() - tic;
time.stop();
double elapsed = time.get_length();
int nb = logger.get_nb_mvp();
std::cout<<"Return value / max_mvp : " << nb << "/" << maxMVP << std::endl;
......
......@@ -63,9 +63,9 @@ int main(int argc, char *argv[])