Commit a828dba6 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

beginning of QRDR implementation

parent 6bd09317
......@@ -15,21 +15,21 @@ data$ID <- basename(file_path_sans_ext(args[1]))
if (length(args) >= 2) {
for (i in 2:length(args)) {
d <- read.table(args[i], header=T)
d$DI <- basename(file_path_sans_ext(args[i]))
d$ID <- basename(file_path_sans_ext(args[i]))
data <- rbind(data, d)
}
}
id = data$DI;
data$ID <- paste(id, "minRes")
id = data$ID;
data$ID1 <- paste(id, "minRes")
data$ID2 <- paste(id, "maxRes")
data$ID3 <- paste(id, "minRealRes")
data$ID4 <- paste(id, "maxRealRes")
p <- ggplot(data, aes(nb_mvp)) +
geom_line(aes(y=minRes, color=ID)) +
geom_line(aes(y=maxRes, color=ID2))+
geom_line(aes(y=minRealRes, color=ID3))+
geom_line(aes(y=minRes, color=ID1)) +
geom_line(aes(y=maxRes, color=ID2)) +
geom_line(aes(y=minRealRes, color=ID3)) +
geom_line(aes(y=maxRealRes, color=ID4))
p + scale_y_log10()
......@@ -194,8 +194,9 @@ public:
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block R0 = _base.get_W(_nbRHS); _base.increase(_nbRHS);
Block R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
R1.realloc(_nbRHS, _nbRHS);
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
}
......@@ -228,7 +229,6 @@ public:
_lastY = Y;
auto MinMaxConv = res.check_precision(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);
......
#ifndef FABULOUS_ARNOLDI_QRDR_H
#define FABULOUS_ARNOLDI_QRDR_H
namespace fabulous {
template<class> class Arnoldi_QRDR;
};
#include "Timer.hpp"
#include "Base.hpp"
#include "Block.hpp"
#include "HessQRDR.hpp"
#include "Utils.hpp"
#include "Logger.hpp"
#include "Arnoldi_Ortho.hpp"
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
*
* \warning This class does NOT support DeflatedRestarting (not implemented)
*/
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";
}
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, const 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, const Block &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,
const std::vector<primary_type> &epsilon)
{
if (!projected_residual_converged)
return false;
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution, Y);
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
auto MinMaxConv = R.check_precision(epsilon, A);
if (std::get<2>(MinMaxConv)) {
_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_QRDR(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),
_solution{},
_lastY{},
_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 Needed for computing real solution if convergence
* \param[in] B Needed for computing real residual if convergence
* \param max_krylov_space_size Maximum size of Krylov Search space.
* \param[in] epsilon tolerance for Residual
* \param ortho_scheme orthogonalization schma
* \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 Struct containing information about the procedure
*/
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 &restarter,
OrthoScheme ortho_scheme = OrthoScheme::MGS,
OrthoType ortho_type = OrthoType::BLOCK)
{
assert( _nbRHS == B.get_nb_col() );
print_header(max_krylov_space_size);
Block Lambda1{};
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
restarter.run(Lambda1, _hess);
_hess.factorize_last_column();
Block ori = _hess.get_original();
check_arnoldi_formula(A, _base, ori);
} else {
Block R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
Lambda1.realloc(_nbRHS, _nbRHS);
A.QRFacto(R0, Lambda1); // R0.InPlaceQRFacto_UserDot(R1, A);
}
_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);
_logger.log_real_residual(A, B, X);
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();
Block Y{_hess.get_nb_vect(), _nbRHS};
Block PR{_hess.get_nb_vect()+_nbRHS, _nbRHS};
_hess.solve_least_square(Y); _lastY = Y;
_hess.compute_proj_residual(Y, PR);
restarter.save_LS(PR);
auto MinMaxConv = _hess.check_convergence(epsilon);
auto proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(
proj_convergence, A, B, X, Y, epsilon);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, Y);
Block ori = _hess.get_original();
check_arnoldi_formula(A, _base, ori);
}
FABULOUS_DEBUG("loop end: "<<_iteration_count);
if (_solution_computed)
X.copy(_solution);
else if (_iteration_count)
compute_solution(A, X, _lastY);
#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
Block ori = _hess.get_original();
restarter.save_hess(ori);
print_footer();
return convergence;
}
};
}; // namespace fabulous
#endif // FABULOUS_ARNOLDI_QRDR_H
......@@ -49,11 +49,11 @@ public:
private:
void print_header(int nb_mvp, std::ostream &o = std::cout)
void print_header(int max_mvp, std::ostream &o = std::cout)
{
o << "#######################################################\n";
o << "################# Arnoldi with QR Inc ################\n";
o << "######## Mat Vect product scheduled "<< nb_mvp <<"###########\n";
o << "######## Mat Vect product scheduled "<< max_mvp <<" ###########\n";
}
void print_footer(std::ostream &o = std::cout)
......@@ -180,6 +180,7 @@ public:
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);
......
......@@ -52,7 +52,7 @@ private:
int _k; /*!<Number of eigen pair to keep*/
S _target; /*!< targeted eigen value*/
Base<S> &_base; /*!< Last Base generated */
HessStandard<S> _hess; /*!< Last Hessenberg generated*/
Block<S> _saved; /*!< Last Hessenberg generated*/
Block<S> _LS; /*!< Least square residual*/
public:
......@@ -71,7 +71,7 @@ public:
operator bool () { return true; }
Base<S> &get_base() { return _base; }
void save_LS(Block<S> &LS) { _LS = LS; }
template< class HESS > void save_hess(HESS &hess) { _hess = hess; }
template< class HESS > void save_hess(HESS &hess) { _saved = hess; }
/**
* \brief This function handle the deflation of eigenvalues at restart.
......@@ -87,7 +87,7 @@ public:
int nm = _base.get_nb_vect() - p;
int dim = _base.get_nb_row();
assert( hess.get_ptr() != _hess.get_ptr() ) ;
assert( hess.get_ptr() != _saved.get_ptr() ) ;
std::cout<<"About to DR ################## !!!\n"
<<"Cste ::\n"
......@@ -95,11 +95,9 @@ public:
<<"\ttarget :: "<<_target<<"\n"
<<"\tp :: "<<p<<"\n"
<<"\tnm :: "<<nm<<"\n"
<<"\tnbCol in Hess :: "<<_hess.get_nb_col()<<"\n"
<<"\tnb_vect in Hess :: "<<_hess.get_nb_vect()<<"\n"
<<"\tnb_hessline in Hess :: "<<_hess.get_nb_hess_line()<<"\n"
<<"\tnbCol in Hess :: "<<_saved.get_nb_col()<<"\n"
<<"\t_Base nbCols :: "<<_base.get_nb_vect()<<"\n"
<<"\tldh :: "<<_hess.get_leading_dim()<<"\n"
<<"\tldh :: "<<_saved.get_leading_dim()<<"\n"
<<"\tdim :: "<<dim<<"\n"
<<"\tLS : leading dim : "<<_LS.get_leading_dim()<<"\n"
<<"\tLS : nbrow : "<<_LS.get_nb_row()<<"\n"
......@@ -110,15 +108,15 @@ public:
LapackKernI::Tgemm( // A <- Hm_^{h} * Hm_
nm, nm, nm+p,
_hess.get_ptr(), _hess.get_leading_dim(),
_hess.get_ptr(), _hess.get_leading_dim(),
_saved.get_ptr(), _saved.get_leading_dim(),
_saved.get_ptr(), _saved.get_leading_dim(),
A.get_ptr(), A.get_leading_dim(),
S{1.0}, S{0.0}
);
//B <- Hm^{h}
for (int j = 0; j < nm; ++j)
for (int i = 0; i < nm; ++i)
B.at(i,j) = std::conj(_hess.at(j, i));
B.at(i,j) = std::conj(_saved.at(j, i));
Block<S> vleft{}; // hold left eigen vectors that won't be computed
Block<S> vright{nm, nm}; // hold right eigen vectors before sorting them
......@@ -159,6 +157,7 @@ public:
Block<S> Q = G;
Block<S> R{newK+p, newK+p};
Q.InPlaceQRFacto(R);
//Q.display("thiz iz Q:");
// Copy last columns of R inside R1
Block<S> Rleft = R.sub_block(0, newK, newK+p, p);
......@@ -188,7 +187,7 @@ public:
Block<S> tmp_HmXQg{nm+p, newK};
LapackKernI::gemm(
nm+p, newK, nm,
_hess.get_ptr(), _hess.get_leading_dim(),
_saved.get_ptr(), _saved.get_leading_dim(),
Q.get_ptr(), Q.get_leading_dim(),
tmp_HmXQg.get_ptr(), tmp_HmXQg.get_leading_dim(),
S{1.0}, S{0.0}
......
......@@ -36,7 +36,6 @@ void GetKIndexesComplex(Block<std::complex<P>> &alpha, Block<S> &beta,
ind.begin(), ind.end(),
[j](const int& a) -> bool { return (a == j); }) == ind.end())
{
Complex diff = (alpha.at(j,0)/beta.at(j,0)) - target;
if (std::norm(diff) < min) // Check if min
{
......
#ifndef FABULOUS_HESS_DR_QR_H
#define FABULOUS_HESS_DR_QR_H
namespace fabulous {
template<class> class HessQRDR;
};
#include "LapackInterface.hpp"
#include "Block.hpp"
#include "Arithmetic.hpp"
namespace fabulous {
/**
* \brief Hessenberg for Incremental QR version
*
* Format: <br/>
\verbatim
| R X X |
| H R X |
| H R |
| H |
\endverbatim
* with: <br/>
* R Triangular Sup <br/>
* H Householders reflectors (trapezoidal)
*
* \warning the whole matrix is allocated, its size is driven by the restart parameter.
*/
template<class S> class HessQRDR : public Block<S>
{
public:
FABULOUS_INHERITS_BLOCK(S);
private:
Block<S> _untouched;
std::vector<std::vector<S>> _tau;
Block<S> _Lambda1;
Block<S> _Lambda;
int _nbRHS;
int _nb_block_col;
int _nb_vect;
std::vector<int> _block_size; /**< size of each block */
std::vector<int> _block_size_sum; /**< sum of size of each block */
FABULOUS_DELETE_COPY_MOVE(HessQRDR);
Block<S> qr_sub_block(int i, int j)
{
int I = _block_size_sum[i];
int J = _block_size_sum[j];
int M = _block_size[i] + _nbRHS;
int N = _block_size[j];
return this->sub_block( I, J, M, N );
}
Block<S> rhs_sub_block(int i)
{
int I = _block_size_sum[i];
int M = _block_size[i] + _nbRHS;
FABULOUS_DEBUG("rhs sub_block I="<<I<<" M="<<M);
return _Lambda.sub_block(I, 0, M, _nbRHS);
}
Block<S> residual_sub_block()
{
int I = _block_size_sum[_nb_block_col];
FABULOUS_DEBUG("I=" <<I);
return _Lambda.sub_block(I, 0, _nbRHS, _nbRHS);
}
public:
HessQRDR(int max_krylov_space_size, int nbRHS):
super{max_krylov_space_size+nbRHS, max_krylov_space_size},
_untouched{max_krylov_space_size+nbRHS, max_krylov_space_size},
_tau{},
_Lambda{max_krylov_space_size+nbRHS, nbRHS},
_nbRHS{nbRHS},
_nb_block_col{0},
_nb_vect{0},
_block_size{},
_block_size_sum{}
{
int max_nb_iter = (max_krylov_space_size / nbRHS) + 1;
_block_size_sum.emplace_back(0);
_tau.reserve( max_nb_iter );
}
void set_rhs(const Block<S> &Lambda1)
{
_Lambda1 = Lambda1;
_Lambda.copy(Lambda1);
}
private:
Block<S> get_H(Block<S> &B)
{
return B.sub_block(0, _nb_vect - _nbRHS, _nb_vect, _nbRHS);
}
Block<S> get_R(Block<S> &B)
{
return B.sub_block(_nb_vect, _nb_vect - _nbRHS, _nbRHS, _nbRHS);
}
public:
Block<S> get_H() { return get_H(*this); }
Block<S> get_R() { return get_R(*this); }
Block<S> get_original() { return _untouched; };
void factorize_last_column()
{
/* STEP 1: Loop over each block in last column to
apply already computed Householders */
char Trans = Arithmetik<S>::ltrans;
{
// save data before qr inc;
int nb_vect = _block_size_sum[_nb_block_col];
assert( nb_vect == _nb_vect );
int block_size = _block_size[_nb_block_col-1];
Block<S> a = this->sub_block(0, nb_vect-block_size, nb_vect+_nbRHS, block_size);
Block<S> b = _untouched.sub_block(
0, nb_vect-block_size, nb_vect+_nbRHS, block_size);
b.copy(a);
FABULOUS_DEBUG("save to ori: J="<<nb_vect-block_size
<<" M="<<nb_vect+_nbRHS<<" N="<<block_size);
}
int k = _nb_block_col-1;
for (int i = 0; i < k; ++i) {
Block<S> A = this->qr_sub_block(i, i);
Block<S> C = this->qr_sub_block(i, k);
FABULOUS_DEBUG("i= "<<i<<" k="<<k);
FABULOUS_DEBUG("A = ("<<A.get_nb_row()<<", "<< A.get_nb_col()<<")");
FABULOUS_DEBUG("C = ("<<C.get_nb_row()<<", "<< C.get_nb_col()<<")");
S *tau = _tau[i].data();
int err = LapackKernI::ormqr(
Trans, C.get_nb_row(), C.get_nb_col(), A.get_nb_col(),
A.get_ptr(), A.get_leading_dim(), tau,
C.get_ptr(), C.get_leading_dim()
);
if (err != 0)
FABULOUS_FATAL_ERROR("return value of ormqr step "<<i<<" is "<<err);
}
/* STEP 2: Call QR on the very last block of last col */
_tau.emplace_back();
std::vector<S> &tau = _tau.back();
Block<S> A = this->qr_sub_block(k, k);
int err = LapackKernI::geqrf( A.get_nb_row(), A.get_nb_col(),
A.get_ptr(), A.get_leading_dim(), tau);
if (err != 0)
FABULOUS_FATAL_ERROR("return values of geqrf (last block) is"<<err);
/* STEP 3: Apply Q^H generated at step 2 to last block of RHS */
Block<S> C = rhs_sub_block(k);
FABULOUS_DEBUG("ormqr k="<<k<<" M="<<C.get_nb_row()<<" N="
<<C.get_nb_col()<<" P="<<A.get_nb_col());
err = LapackKernI::ormqr(
Trans, C.get_nb_row(), C.get_nb_col(), A.get_nb_col(),
A.get_ptr(), A.get_leading_dim(), tau.data(),
C.get_ptr(), C.get_leading_dim()
);
if (err != 0)
FABULOUS_FATAL_ERROR("return value of ormqr (RHS) is "<<err);
}
void solve_least_square(Block<S> &Y)
{
Block<S> R = this->sub_block(0, 0, _nb_vect, _nb_vect);
Block<S> Lambda = _Lambda.sub_block(0, 0, _nb_vect, _nbRHS);
Y.copy(Lambda);
FABULOUS_DEBUG("Y.row = " << Y.get_nb_row());
FABULOUS_DEBUG("Y.col = " << Y.get_nb_col());
FABULOUS_DEBUG("Lambda.row = " << Lambda.get_nb_row());
FABULOUS_DEBUG("R.row = " << R.get_nb_row());
FABULOUS_DEBUG("R.col = " << R.get_nb_col());
LapackKernI::trsm(
R.get_nb_row(), Lambda.get_nb_col(),
R.get_ptr(), R.get_leading_dim(),
Y.get_ptr(), Y.get_leading_dim()
);
Y.display("fresh Y:");
}
void compute_proj_residual(const Block<S> &Y, Block<S> &PR)
{
int M = _nb_vect + _nbRHS;
int N = _nb_vect;
assert( Y.get_nb_row() == N );
assert( PR.get_nb_row() == M );
PR.copy(_Lambda1);
_Lambda1.display("Hey this is Lambda1:");
assert(PR.get_nb_col() == _Lambda1.get_nb_col());
assert(PR.get_nb_row() >= _Lambda1.get_nb_row());
//FABULOUS_DEBUG("R1.nb_row()="<<_R1.get_nb_row());
LapackKernI::gemm(
M, _nbRHS, N,
_untouched.get_ptr(), _untouched.get_leading_dim(),
Y.get_ptr(), Y.get_leading_dim(),
PR.get_ptr(), PR.get_leading_dim(),
S{-1.0}, S{1.0}
);
//PR.display("eat my PR:");
}
std::tuple<P,P,bool>
check_convergence(const std::vector<P> &epsilon)
{
Block<S> RR = residual_sub_block();
return RR.check_precision(epsilon);
}
void increase(int block_size)
{
_nb_vect += block_size;
++ _nb_block_col;
_block_size.emplace_back(block_size);
_block_size_sum.emplace_back(_nb_vect);
}
void reserve_DR(int nb_eigen_vector)
{
// CHECK this is called just after hessenberg initialization:
assert( _nb_block_col == 0 );
assert( _nb_vect == 0 );
assert( _block_size.size() == 0 );
assert( _block_size_sum.size() == 1 );
increase(nb_eigen_vector);
}
int get_nb_vect() const { return _nb_vect; }
};
}; // namespace fabulous
#endif // FABULOUS_HESS_DR_QR_H
......@@ -179,11 +179,11 @@ struct LapackKernI
* \param m Integer number of rows in C
* \param n Integer number of columns in C
* \param k Integer Number of elementary reflectors in Q
* \param A Ptr array returned by geqrf
* \param lda Integer leading dim of a
* \param tau Ptr array returned by geqrf
* \param A array returned by geqrf
* \param lda leading dim of a
* \param tau array returned by geqrf
* \param C the MxN matrix to be multiplied
* \param ldc Integer leading dimension of C
* \param ldc leading dimension of C
*
* Ref: https://software.intel.com/en-us/node/521011
*
......
......@@ -218,10 +218,10 @@ std::pair<P,P> compare_blocks(const Block &block1, const Block &block2)
* \param base Current base (correspond to V in above formula)
* \param hess Current hessenberg, given through parent references.
*/
template< class Hess, class Matrix, class Base,
template< class Matrix, class Base,
class S = typename Base::value_type,
class P = typename Base::primary_type >
void check_arnoldi_formula(Matrix &A, Base &base, Hess &hess)
void check_arnoldi_formula(Matrix &A, Base &base, Block<S> &hess)
{
using Block = ::fabulous::Block<S>;
......
......@@ -12,6 +12,7 @@
#include "../src/Arnoldi_IB.hpp"
#include "../src/Arnoldi.hpp"
#include "../src/Arnoldi_QRInc.hpp"
#include "../src/Arnoldi_QRDR.hpp"
#include "../src/LapackInterface.hpp"
#define MARK_AS_USED(var_) ((void)(var_))
......
......@@ -28,10 +28,10 @@ int main(int argc, char *argv[])
mml.LoadMatrix(mat, filename);
int dim = mat.size();
int nbRHS = 10;
int nbRHS = 4;
int restart = 0;
int nbEigenVector = 20;
int maxMVP = 10000;
int nbEigenVector = 2;
int maxMVP = 100;
S target{0.0};
if (args.size() > 1)
......@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
// OrthoScheme::IMGS, OrthoType::BLOCK
// );
runTest_BGMRES_DR_filelog<Arnoldi>(
runTest_BGMRES_DR_filelog<Arnoldi_QRDR>(
"young1c",