Commit 08b773c9 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Make HessQR, HessQRDR and Hess_cham_QR follow the same concept

ArnodiStandard now have Hessenberg as a template parameter and
so that HessStandard, HessQR, HessQRDR, and Hess_cham_QR can be
used with ArnoldiStandard
parent 22d5c0ee
......@@ -91,6 +91,7 @@ endif()
include_directories(include)
include_directories(src/api/include)
add_subdirectory(include)
add_subdirectory(src)
if(FABULOUS_BUILD_DOC)
......
#+TITLE: fabulous (Fast Accurate Block Linear krylOv Solver, a.k.a IB-BGMRES-DR)
#+AUTHOR: Thomas Mijieux
#+AUTHOR: Cyrille Piacibello
#+AUTHOR: Luc Giraud
* Purpose
Block Krylov iterative solver API.
This Library is implementing Block Krylov General Minimum Residual iterative solver
(denoted as BGMRES) which means that this library can solve linear equations
with multiple right hand sides.
This library also implements several algorithm variants that may improve the
performance of the library (Incremental QR factorization, IB-BGMRES, BGMRES-DR,
and IB-BGMRES-DR (not implemented yet)
* References
IB-BGMRES: M. Robbé and M. Sadkane.
Exact and inexact breakdowns in the block GMRES method.
Linear Algebra and its Applications, 419:265–285, 2006
BGMRES-DR: R. B. Morgan.
Restarted block GMRES with deflation of eigenvalues.
Applied Numerical Mathematics, 54(2):222–236, 2005.
IB-BGMRES-DR: E. Agullo, L. Giraud, Y.-F. Jing
Block GMRES method with inexact breakdowns and deflated restarting
SIAM Journal on Matrix Analysis and Applications
35, 4, November 2014, p. 1625–1651
* Details
The Matrix x Vector product is not handled by the library and must be
implemented as a callback by the user.
This allow the user to implement a algorithm adapted to his/her problem.
(whether the matrix is sparse/dense or have a specific structure)
The user must also implements a DotProduct callback.
This allow the user to use the library in a distributed fashion.
(The input solution and right hand sides may be distributed with block
lines style.) If the user intend to distribute its vector over several nodes,
the callbacks must be adapted to perform the necessary communications
and/or reductions.
* Installing
** Dependencies
*** LAPACKE (LAPACK, CBLAS and BLAS)
Any implementation of C interface for LAPACK and BLAS (intel mkl, netlib, openblas, eigen, ...)
*** chameleon (optional)
As fabulous is a C++ library, it requires that chameleon supports redefining
complex type to C++ complex types in its 'client' headers. for C++ application.
This is currently supported in a fork of the main chameleon development branch,
(hopefully merge any time soon)
#+BEGIN_SRC sh
cd ..
git clone https://gitlab.inria.fr/tmijieux/chameleon.git
#+END_SRC
** fabulous
A basic installation can be performed with the following commands:
#+BEGIN_SRC sh
prefix="/usr/"
mkdir -p build/
rm -rf build/*
cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix}
make
make install
#+END_SRC
Some options can be passed to the compilation stage:
#+BEGIN_SRC sh
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DFABULOUS_BUILD_SHARED_LIBS=ON \
-DFABULOUS_BUILD_SHARED_LIBS=ON -DFABULOUS_BUILD_DOC=ON \
-DFABULOUS_DEBUG_MODE=OFF -DFABULOUS_LAPACKE_NANCHECK=OFF
#+END_SRC
......@@ -32,7 +32,7 @@ DOXYFILE_ENCODING = UTF-8
# title of most generated pages and in a few other places.
# The default value is: My Project.
PROJECT_NAME = "Ib-BGMRes-Dr"
PROJECT_NAME = "Fabulous"
# The PROJECT_NUMBER tag can be used to enter a project or revision number. This
# could be handy for archiving the generated documentation or if some version
......@@ -771,7 +771,7 @@ WARN_LOGFILE =
# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
# Note: If this tag is empty the current directory is searched.
INPUT = @CMAKE_CURRENT_SOURCE_DIR@/../src/ @CMAKE_CURRENT_SOURCE_DIR@/../Api/src @CMAKE_CURRENT_SOURCE_DIR@/src_dox
INPUT = @CMAKE_CURRENT_SOURCE_DIR@/../src/api/ @CMAKE_CURRENT_SOURCE_DIR@/../src/test_api/ @CMAKE_CURRENT_SOURCE_DIR@/../include/fabulous/ @CMAKE_CURRENT_SOURCE_DIR@/src_dox
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
......@@ -802,7 +802,7 @@ FILE_PATTERNS =
# be searched for input files as well.
# The default value is: NO.
RECURSIVE = NO
RECURSIVE = YES
# The EXCLUDE tag can be used to specify files and/or directories that should be
# excluded from the INPUT source files. This way you can easily exclude a
......@@ -838,7 +838,7 @@ EXCLUDE_PATTERNS =
# Note that the wildcards are matched against the file with absolute path, so to
# exclude all test directories use the pattern */test/*
EXCLUDE_SYMBOLS =
EXCLUDE_SYMBOLS = fabulous::experimental
# The EXAMPLE_PATH tag can be used to specify one or more files or directories
# that contain example code fragments that are included (see the \include
......
/*!
\mainpage Ib-BGMRes-Dr: Inexact Breakdown - BlockGMRes - Deflated Restarting.
\mainpage Fabulous (Fast Accurate Block Linear krylOv Solver)
(IB-BGMRES-DR): Inexact Breakdown - BlockGMRes - Deflated Restarting.
This library implements a block GMRes linear solver for the solution of
multi right-hand sides given simulateneously.
......
set(FABULOUS_XX_LIB_SRC
Algorithm.hpp
Arithmetic.hpp
Arnoldi.hpp
Arnoldi_IB.hpp
Arnoldi_Ortho_BLOCK.hpp
Arnoldi_Ortho_Enum.hpp
Arnoldi_Ortho.hpp
Arnoldi_Ortho_RUHE.hpp
Arnoldi_Ortho_RUHE_IB.hpp
Arnoldi_QRDR.hpp
Arnoldi_QRInc.hpp
Base.hpp
BGMResDR.hpp
BGMRes.hpp
Block.hpp
BlockWP.hpp
Chameleon.hpp
ChameleonInterface.hpp
Color.hpp
DeflatedRestart.hpp
EigenUtils.hpp
Error.hpp
HessExtended.hpp
HessQRDR.hpp
HessQR.hpp
HessStandard.hpp
KrylovAlgo.hpp
LapackInterface.hpp
Logger.hpp
Timer.hpp
Utils.hpp
)
install(
FILES ${FABULOUS_XX_LIB_SRC}
DESTINATION include/fabulous
DIRECTORY fabulous
DESTINATION include
)
......@@ -91,8 +91,10 @@ private:
}
template<class Matrix, class Block>
void compute_solution(Matrix &A, Block &X, Block &Y)
void compute_solution(Matrix &A, Block &X)
{
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
}
......@@ -121,7 +123,7 @@ private:
template<class Matrix, class Block>
bool check_real_residual_if(bool projected_residual_converged,
Matrix &A, Block &X, Block &Y, const Block &B,
Matrix &A, Block &X, const Block &B,
const std::vector<primary_type> &epsilon)
{
if (!projected_residual_converged)
......@@ -129,7 +131,7 @@ private:
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution, Y);
compute_solution(A, _solution);
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
......@@ -268,10 +270,15 @@ private:
template<class Matrix>
void handle_no_convergence(
int proj_convergence, BlockWP<S> &WP, Block<S> &PR, Block<S> &Y,
int proj_convergence, BlockWP<S> &WP,
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> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
......@@ -288,7 +295,7 @@ private:
// relation would be false from the begining
_solution = Block<S>{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution, Y);
compute_solution(A, _solution);
_solution_computed = true;
_cold_restart = true;
}
......@@ -399,39 +406,29 @@ public:
// Ortho and filling of L (the (IB) Hessenberg)
Arnoldi_Ortho(L, _base, WP, ortho_type, ortho_scheme, A);
// WP.check_ortho_WP();
// WP.get_W().check_ortho();
// WP.check_ortho_WP();// WP.get_W().check_ortho();
//L.display_hess_extended_bitmap();
Block PR{L.get_nb_hess_line(), _nbRHS};
Block YY{L.get_nb_hess_line(), _nbRHS};
Block Y = YY.sub_block(0, 0, L.get_nb_vect(), _nbRHS);
Block res = YY.sub_block(L.get_nb_vect(), 0, _nbRHS, _nbRHS);
_last_Y = Y;
// Y is the solution of LeastSquare problem: Y = argmin(||F*Y-\Lambda||)
// PR is the projected residual PR <- F*Y-\Lambda
L.compute_proj_residual(YY, PR);
//auto MinMaxConv = PR.check_precision(epsilon);
auto MinMaxConv = res.check_precision(epsilon);
auto MinMaxConv = L.check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, X, Y, B, epsilon);
convergence = check_real_residual_if(proj_convergence, A, X, B, epsilon);
if (!convergence) {
handle_no_convergence(
proj_convergence, WP, PR, Y, epsilon, inv_epsilon, A, X);
handle_no_convergence(proj_convergence, WP,
epsilon, inv_epsilon, A, X);
}
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, WP.get_size_W(), MinMaxConv);
_logger.log_real_residual(A, B, X, _base, Y);
_logger.log_real_residual(A, B, X, _base, L);
}
if (_solution_computed)
X.copy(_solution);
else
compute_solution(A, X, _last_Y);
compute_solution(A, X);
#ifdef FABULOUS_DEBUG_MODE
if (convergence) {
......
......@@ -2,7 +2,7 @@
#define FABULOUS_ARNOLDI_HPP
namespace fabulous {
template<class> class Arnoldi;
template<template<class> class HESSENBERG, class S > class Arnoldi;
};
#include "fabulous/data/Base.hpp"
......@@ -23,7 +23,7 @@ namespace fabulous {
*
* This class support DeflatedRestarting
*/
template< class S > class Arnoldi
template<template<class> class HESSENBERG, class S > class Arnoldi
{
public:
using value_type = typename Arithmetik<S>::value_type;
......@@ -47,7 +47,7 @@ private:
private:
Logger<P> &_logger;
HessStandard<S> _hess;
HESSENBERG<S> _hess;
bool _solution_computed;
Block<S> _solution;
Block<S> _lastY;
......@@ -107,16 +107,17 @@ private:
_nb_mvp += V.get_nb_col();
}
template<class Matrix, class Block>
void compute_solution(Matrix &A, Block &X, const Block &Y)
template<class Matrix>
void compute_solution(Matrix &A, Block<S> &X)
{
Block<S> Y = _hess.alloc_least_square_sol();
_hess.solve_least_square(Y);
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
}
template<class Matrix, class Block>
bool check_real_residual_if(bool projected_residual_converged,
Matrix &A, const Block &B,
Block &X, const Block &Y,
Matrix &A, const Block &B, Block &X,
const std::vector<P> &epsilon)
{
if (!projected_residual_converged)
......@@ -125,14 +126,17 @@ private:
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution, Y);
compute_solution(A, _solution);
Block R{_dim, _nbRHS};
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);
printf("real residual=(%g;%g)\n", min, max);
#endif
if (std::get<2>(MinMaxConv)) {
_solution_computed = true;
......@@ -200,6 +204,7 @@ public:
_base.increase(_nbRHS);
R1.realloc(_nbRHS, _nbRHS);
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
}
auto MinMaxConv = R1.check_precision(epsilon);
bool convergence = std::get<2>(MinMaxConv);
......@@ -207,7 +212,6 @@ public:
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X);
_hess.set_rhs(R1);
while (!convergence && _size_krylov_space < max_krylov_space_size) {
++ _iteration_count;
......@@ -222,33 +226,31 @@ public:
_size_krylov_space = _base.get_nb_vect();
//_hess.display_hess_bitmap("Hess Before Least Square");
int nb_vect = _hess.get_nb_vect();
Block YY{nb_vect+_nbRHS, _nbRHS};
Block Y = YY.sub_block(0, 0, nb_vect, _nbRHS);
Block res = YY.sub_block(nb_vect, 0, _nbRHS, _nbRHS);
_hess.compute_least_square_sol(YY);
_lastY = Y;
auto MinMaxConv = res.check_precision(epsilon);
auto MinMaxConv = _hess.check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
Block PR{_hess.get_nb_vect()+_nbRHS, _nbRHS};
_hess.compute_proj_residual(Y, PR);
restarter.save_LS(PR);
convergence = check_real_residual_if(proj_convergence, A, B, X, Y, epsilon);
convergence = check_real_residual_if(proj_convergence, A, B, X, epsilon);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, Y);
_logger.log_real_residual(A, B, X, _base, _hess);
}
if (_solution_computed)
X.copy(_solution);
else if (_iteration_count)
compute_solution(A, X, _lastY);
compute_solution(A, X);
if (_cold_restart)
if (_cold_restart) {
_base.reset();
} else if (restarter && !convergence) {
int M = _hess.get_nb_vect() + _nbRHS;
Block LS{M, _nbRHS};
_hess.compute_proj_residual(LS);
restarter.save_LS(LS);
Block H = _hess.get_hess_block();
restarter.save_hess(H);
}
print_footer();
restarter.save_hess(_hess);
return convergence;
}
};
......
#ifndef FABULOUS_ARNOLDI_QRINC_HPP
#define FABULOUS_ARNOLDI_QRINC_HPP
namespace fabulous {
template<class> class Arnoldi_QRInc;
};
#include "fabulous/utils/Timer.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/data/Base.hpp"
#include "fabulous/data/Block.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/hessenberg/Hess_Std_QR.hpp" // TODO: templatize HESSENBERGS
namespace fabulous {
/**
* \brief %Arnoldi iterations With incremental QR
*
* Incremental QR: GEQRF, ORMQR and TRSM kernels are used to solve
* the LeastSquare problem on the projected matrix
*
* \warning This class does NOT support DeflatedRestarting (not implemented)
*/
template<class S> class Arnoldi_QRInc {
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
private:
Logger<P> &_logger;
HessQR<S> _hess;
bool _solution_computed;
Block<S> _solution;
Base<S> &_base;
int _dim;
int _nbRHS;
int _size_krylov_space;
int _nb_mvp;
int _iteration_count;
public:
int get_krylov_space_size() const { return _size_krylov_space; }
int get_nb_mvp() const { return _nb_mvp; }
int get_nb_iteration() const { return _iteration_count; }
private:
void print_header(int max_mvp, std::ostream &o = std::cout)
{
o << "#######################################################\n";
o << "################# Arnoldi with QR Inc ################\n";
o << "######## Mat Vect product scheduled "<< max_mvp <<" ###########\n";
}
void print_footer(std::ostream &o = std::cout)
{
o << "################# Iterations done: "<<_iteration_count<<"(+1)\n";
o << "#######################################################\n";
}
template< class Matrix, class Block >
void W_AMV(Matrix &A, Block &V, Block &W) // W = A*M*V
{
if (A.useRightPrecond()) {
Block MV{W.get_nb_row(), V.get_nb_col()};
A.PrecondBlockVect(V, MV); // MV = M * V
A.MatBlockVect(MV, W); // B = A*MV
} else {
A.MatBlockVect(V, W); // W = A * V
}
_nb_mvp += V.get_nb_col();
}
// X_n = X_0 + M * V * Y
template< class Matrix, class Block, class Base >
void X_plus_MVY(Matrix &A, Block &X, Base &V, Block &Y)
{
int dim = V.get_nb_row();
Block MVY{dim, Y.get_nb_col()};
if (A.useRightPrecond()) {
Block VY{dim, Y.get_nb_col()};
V.compute_product(Y, VY); // VY = V*Y
A.PrecondBlockVect(VY, MVY); // MVY = M*V*Y
} else {
V.compute_product(Y, MVY); // MVY = V*Y
}
// X_n = X_0 + MVY
for (int j = 0; j < X.get_nb_col(); ++j)
LapackKernI::axpy(dim, S{1.0}, MVY.get_vect(j), 1, X.get_vect(j), 1);
}
template<class Matrix, class Block>
void compute_solution(Matrix &A, Block &X)
{
Block Y{_hess.get_nb_vect(), _nbRHS};
// in Chameleon QR version: this function can only be called once:
_hess.solve_least_square(Y);
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
}
template<class Matrix, class Block>
bool check_real_residual_if(bool projected_residual_converged,
Matrix &A, Block &X, const Block &B,
const std::vector<primary_type> &epsilon)
{
if (!projected_residual_converged)
return false;
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution);
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
if (std::get<2>(R.check_precision(epsilon, A))) {
_solution_computed = true;
return true;
}
return false;
}
template<class Matrix, class Block>
void compute_real_residual(Matrix &A, Block &X, const Block &B, Block &R)
{
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();
}
public:
template<class Restarter>
Arnoldi_QRInc(Logger<primary_type> &logger, Restarter &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_hess{max_krylov_space_size, nbRHS},
_solution_computed(false),
_base{restarter.get_base()},
_dim(dim),
_nbRHS(nbRHS),
_size_krylov_space{0},
_nb_mvp{0},
_iteration_count{0}
{
}
/**
* \brief Arnoldi iterations with incremental QR for least square resolution
*
* Solution will be stored inside X0 at the end of computation.
* Ref IB-BGMRES-DR (Giraud Jing Agullo): p.21
*
* \param[in] A User Matrix
* \param[in,out] X initial guess (input). Best approx. solution (output)
* \param[in] B the right hand sides
* \param max_krylov_space_size Maximum size of 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) or vector wise arnoldi (RUHE). <br/>
* (In distributed, only the vector wise will works)
*
* \return whether the convergence was reached
*/
template<class Matrix, class Block, class Restarter>
bool run( const Matrix &A, Block &X, const Block &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter&,
OrthoScheme ortho_scheme,
OrthoType ortho_type )
{
assert( _nbRHS == B.get_nb_col() );
_base.reset(); // DR not implemented ; force a cold restart
print_header(max_krylov_space_size);
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
Block R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
Block Lambda1{_nbRHS, _nbRHS};
A.QRFacto(R0, Lambda1); // R0.InPlaceQRFacto_UserDot(R1, A);
_base.increase(_nbRHS);
_hess.set_rhs(Lambda1); // Give R part to Hessenberg; to be used in LeastSquare.
auto MinMaxConv = Lambda1.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, _nbRHS, MinMaxConv);
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();
Block W = _base.get_W(_nbRHS);
W_AMV(A, Vj, W); // W = A*M*Vj
Arnoldi_Ortho(_hess, _base, W, ortho_type, ortho_scheme, A); //W.check_ortho();
_base.increase(_nbRHS);
_size_krylov_space = _base.get_nb_vect();
_hess.factorize_last_column();
auto MinMaxConv = _hess.check_convergence(epsilon);
auto proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(
proj_convergence, A, X, B, epsilon);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
//check_arnoldi_formula(A, _base, hess);
}
FABULOUS_DEBUG("loop end: "<<_iteration_count);
if (_solution_computed)
X.copy(_solution);
else if (_iteration_count)
compute_solution(A, X);
#ifdef FABULOUS_DEBUG_MODE
if (convergence) { // FIXME; this is a debug test
auto mm = eval_current_solution(A, B, X);
printf("real residual=(%g;%g)\n", mm.first, mm.second);
}
#endif
print_footer();
return convergence;
}
};
}; // namespace fabulous
#endif // FABULOUS_ARNOLDI_QRINC_HPP
#ifndef FABULOUS_ARNOLDI_QRDR_H
#define FABULOUS_ARNOLDI_QRDR_H
namespace fabulous {
template<class> class Arnoldi_QRDR;
};
#include "fabulous/data/Base.hpp"
#include "fabulous/data/Block.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/Timer.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/hessenberg/Hess_Std_QRDR.hpp" // TODO templatize HESSENBERGS
namespace fabulous {
/**
* \brief %Arnoldi iterations With incremental QR and Deflated Restarting
*
* Incremental QR: GEQRF, ORMQR and TRSM kernels are used to solve
* the LeastSquare problem on the projected matrix
*/
template<class S> class Arnoldi_QRDR {
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
private:
Logger<P> &_logger;
HessQRDR<S> _hess;
bool _solution_computed;
Block<S> _solution;
Block<S> _lastY;
Base<S> &_base;
int _dim;
int _nbRHS;
int _size_krylov_space;
int _nb_mvp;
int _iteration_count;
public:
int get_krylov_space_size() const { return _size_krylov_space; }
int get_nb_mvp() const { return _nb_mvp; }
int get_nb_iteration() const { return _iteration_count; }
private:
void print_header(int max_mvp, std::ostream &o = std::cout)
{
o << "#######################################################\n";
o << "################# Arnoldi with QR DR #################\n";
o << "######## Mat Vect product scheduled "<< max_mvp <<" ###########\n";
}
void print_footer(std::ostream &o = std::cout)
{
o << "################# Iterations done: "<<_iteration_count<<"(+1)\n";
o << "#######################################################\n";
}