Commit a1cce305 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

fix Hm_ transpose bug(DeflatedRestart) and update R plot script

parent 70e62ca4
......@@ -6,15 +6,30 @@ if (length(args) < 1) {
}
library(ggplot2)
library(reshape2)
library(tools)
data <- read.table(args[1], header=T)
data$ID <- basename(file_path_sans_ext(args[1]))
if (length(args) >= 2) {
for (i in 2:length(args)) {
data <- rbind(data, read.table(args[i], header=T))
d <- read.table(args[i], header=T)
d$DI <- basename(file_path_sans_ext(args[i]))
data <- rbind(data, d)
}
}
p <- ggplot(data, aes(nb_mvp), log10='y') +
id = data$DI;
data$ID <- 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=ID))
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()
......@@ -8,7 +8,6 @@ template<class> class Arnoldi;
#include "Base.hpp"
#include "Block.hpp"
#include "HessStandard.hpp"
#include "HessStorage.hpp"
#include "Utils.hpp"
#include "Logger.hpp"
#include "DeflatedRestart.hpp"
......@@ -140,9 +139,9 @@ private:
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
auto MinMaxConv = R.check_precision(epsilon, A);
printf("real residual=(%g;%g)\n",
std::get<0>(MinMaxConv),
std::get<1>(MinMaxConv));
auto min = std::get<0>(MinMaxConv);
auto max = std::get<1>(MinMaxConv);
printf("real residual=(%g;%g)\n", min, max);
if (std::get<2>(MinMaxConv)) {
_solution_computed = true;
......@@ -158,7 +157,7 @@ public:
Arnoldi(Logger<P> &logger, Restarter &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_hess{max_krylov_space_size+1, nbRHS},
_hess{max_krylov_space_size, nbRHS},
_solution_computed(false),
_base(restarter.get_base()),
_dim(dim),
......@@ -185,7 +184,6 @@ public:
* or vector wise arnoldi. (In distributed, only the vector wise will works)
*
* \return information about the iterations
*
*/
template<class Matrix, class Block, class Restarter>
ArnReturn run( const Matrix &A, Block &X, const Block &B,
......@@ -197,46 +195,40 @@ public:
{
print_header(max_krylov_space_size);
Block R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
bool convergence = false;
Block R1{}; // R1 = RHS of projected problem
int iteration_count = 0;
_logger.notify_iteration_begin(iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
int k = restarter.run(R1);
restarter.restore_hess(_hess, k, _nbRHS);
_size_krylov_space = k;
std::cout<<Color::bold<<Color::yellow;
std::cout<<"\t\tDR done\n"<<Color::reset;
restarter.run(R1, _hess);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_logger.notify_iteration_begin(iteration_count, _size_krylov_space);
Block R0 = _base.get_W(_nbRHS);
Block R0 = _base.get_W(_nbRHS); _base.increase(_nbRHS);
compute_real_residual(A, X, B, R0);
/* auto minmaxR0 = R0.get_min_max_norm(A); */
A.QRFacto(R0, R1);
auto MinMaxConv = R1.check_precision(epsilon);
_base.increase(_nbRHS);
_size_krylov_space = _base.get_nb_vect();
convergence = std::get<2>(MinMaxConv);
_logger.notify_iteration_end(
iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
R1.realloc(_nbRHS, _nbRHS);
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
}
auto 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, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X);
_hess.set_rhs(R1);
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);
Block Vj = _base.get_Vj(); // V_j
assert( Vj.get_nb_col() == _nbRHS );
Block W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
Arnoldi_Ortho(_hess, _base, W, ortho_type, ortho_scheme, A);
// W.check_ortho();
//W.check_ortho("Wcandidate");
_base.increase(_nbRHS); // add W to base
_base.increase(_nbRHS);
_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);
......@@ -248,21 +240,23 @@ public:
Block PR{_hess.get_nb_vect()+_nbRHS, _nbRHS};
_hess.compute_proj_residual(Y, PR);
restarter.store_LS(PR);
auto MinMaxConv2 = PR.check_precision(epsilon);
FABULOUS_DEBUG("PR min: "<<std::scientific<<std::get<0>(MinMaxConv2)
<<" max: "<<std::scientific<<std::get<1>(MinMaxConv2));
restarter.save_LS(PR);
convergence = check_real_residual_if(proj_convergence, A, B, X, Y, epsilon);
//check_arnoldi_formula(A, _base, _hess);
_logger.notify_iteration_end(
iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, Y);
// check_arnoldi_formula(A, _base, _hess);
}
// No convergence, we need to write current sol into X0
if (_solution_computed)
X.copy(_solution);
else
compute_solution(A, X, _lastY);
if (iteration_count) {
if (_solution_computed)
X.copy(_solution);
else
compute_solution(A, X, _lastY);
}
if (_cold_restart)
_base.reset();
......
......@@ -6,7 +6,6 @@
#include "HessExtended.hpp"
#include "Utils.hpp"
#include "Logger.hpp"
#include "HessStorage.hpp"
namespace fabulous {
......@@ -380,6 +379,7 @@ public:
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X);
if (_nb_direction_kept == 0)
return ArnReturn{0, 0, true};
......
......@@ -145,12 +145,14 @@ void ArnoldiBlock_IMGS(Base &base, Block &H, Block &W)
assert( H.get_nb_col() == W.get_nb_col() );
assert( H.get_nb_row() == base.get_nb_vect() );
FABULOUS_DEBUG("base.nb_vect="<<base.get_nb_vect());
// Loop over different blocks
int size_block_sum = 0;
for (int k = 0; k < nb_block; ++k) {
int size_block = base.get_block_size(k);
Block tmp{size_block, W_size};
FABULOUS_DEBUG("base.block["<<k<<"].size="<<size_block);
LapackKernI::Tgemm( // tmp = V_i^{h} * W
size_block, W_size, dim,
......
......@@ -76,14 +76,12 @@ public:
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
inv_normB = array_inverse(normB);
X.scale(inv_normB);
B.scale(inv_normB);
X.scale(inv_normB); B.scale(inv_normB);
bool finished = false;
while (!finished && _krylov_space_spanned < maxMVP) {
while (!finished && _mvp < maxMVP) {
//Compute nb of mat vect product to give to Arnoldi procedure
int size_to_span = std::min(max_krylov_space_size, maxMVP-_krylov_space_spanned);
int size_to_span = std::min(max_krylov_space_size, maxMVP-_mvp);
size_to_span = std::max(size_to_span, nbRHS);
print_iteration_start_info(size_to_span);
......@@ -95,15 +93,13 @@ public:
finished = r.has_converged;
_nb_iteration += r.nb_iteration;
_krylov_space_spanned += r.size_krylov_space;
_mvp += r.size_krylov_space;
++ _nb_restart;
print_iteration_end_info(r.size_krylov_space);
}
X.scale(normB);
B.scale(normB);
return _krylov_space_spanned;
X.scale(normB); B.scale(normB);
return _mvp;
}
};
}; // namespace fabulous
......
......@@ -7,7 +7,6 @@
#include "Block.hpp"
#include "Utils.hpp"
#include "Arnoldi.hpp"
#include "HessStorage.hpp"
#include "DeflatedRestart.hpp"
#include "Arithmetic.hpp"
#include "Color.hpp"
......@@ -84,8 +83,7 @@ public:
* (cumulated size of the spanned Krylov Spaces)
*
*/
template<template <class> class ARNOLDI, class Matrix, class Block,
class P = typename Block::primary_type >
template<template <class> class ARNOLDI, class Matrix, class Block >
int solve_block(Matrix &A, Block &B, Block &X,
const int maxMVP, const int max_krylov_space_size,
const std::vector<P> &epsilon,
......@@ -95,7 +93,6 @@ public:
{
const int nbRHS = B.get_nb_col();
const int dim = B.get_nb_row();
reset();
print_start_info( dim, nbRHS, maxMVP, max_krylov_space_size,
ortho_scheme, ortho_type, epsilon, A );
......@@ -104,18 +101,14 @@ public:
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
inv_normB = array_inverse(normB);
X.scale(inv_normB); B.scale(inv_normB);
X.scale(inv_normB);
B.scale(inv_normB);
Base<S> base{dim, max_krylov_space_size+nbRHS+1};
HessStorage<S> hessSave;
DeflatedRestarting<S> restarter{nb_eigen_pair, target, base, hessSave};
Base<S> base{dim, max_krylov_space_size+nbRHS};
DeflatedRestarting<S> restarter{nb_eigen_pair, target, base};
bool finished = false;
while (!finished && _krylov_space_spanned < maxMVP) { // Main loop
while (!finished && _mvp < maxMVP) { // Main loop
// Compute nb of mat vect product to give to Arnoldi procedure
int size_to_span = std::min(max_krylov_space_size, maxMVP-_krylov_space_spanned);
int size_to_span = std::min(max_krylov_space_size, maxMVP-_mvp);
// Test if SizeTospan is big enough to hold 1 block + Deflation directions
if (size_to_span < (nbRHS+nb_eigen_pair)) {
......@@ -123,23 +116,19 @@ public:
break;
}
print_iteration_start_info(size_to_span);
ARNOLDI<S> arnoldi{_logger, restarter, dim, nbRHS, size_to_span};
auto r = arnoldi.run( A, X, B,
size_to_span, epsilon, restarter,
auto r = arnoldi.run( A, X, B, size_to_span, epsilon, restarter,
ortho_scheme, ortho_type );
finished = r.has_converged;
_nb_iteration += r.nb_iteration;
_krylov_space_spanned += r.size_krylov_space;
_mvp += r.size_krylov_space;
++ _nb_restart;
print_iteration_end_info(r.size_krylov_space);
}
X.scale(normB);
B.scale(normB);
return _krylov_space_spanned;
X.scale(normB); B.scale(normB);
return _mvp;
}
};
}; // namespace fabulous
......
......@@ -17,7 +17,7 @@ namespace fabulous {
template<class S> class Base : public Block<S>
{
private:
int _current_block; // id of current used block
int _nb_block; // id of current used block
int _nb_vect;
int _max_krylov_space_size;
std::vector<int> _block_sizes; // (block size may change)
......@@ -34,7 +34,7 @@ public:
*/
Base(int dim, int max_krylov_space_size):
Block<S>{dim, max_krylov_space_size},
_current_block{0},
_nb_block{0},
_nb_vect{0},
_max_krylov_space_size{max_krylov_space_size}
{
......@@ -46,14 +46,14 @@ public:
*/
void reset()
{
_current_block = 0;
_nb_block = 0;
_nb_vect = 0;
_block_sizes.clear();
_block_sizes_sum.clear();
_block_sizes_sum.emplace_back(0);
}
int get_nb_block() const { return _current_block; }
int get_nb_block() const { return _nb_block; }
int get_nb_vect() const { return _nb_vect; }
int get_block_size(int k) const { return _block_sizes[k]; }
......@@ -71,9 +71,14 @@ public:
return sub_block(0, _nb_vect-block_size, get_nb_row(), block_size);
}
void check_ortho(std::string name="Base")
{
sub_block(0, 0, get_nb_row(), _nb_vect).check_ortho(name);
}
void increase(int block_size)
{
++ _current_block;
++ _nb_block;
_nb_vect += block_size;
_block_sizes.push_back(block_size);
_block_sizes_sum.push_back(_nb_vect);
......
......@@ -145,7 +145,7 @@ public:
/* \brief return norm of k-th vector inside the Block
*/
template < class Matrix >
P get_norm(int k, Matrix &A) const
P get_norm(int k, const Matrix &A) const
{
const S *ptr = get_vect(k);
S res = S{0.0};
......@@ -188,7 +188,7 @@ public:
/* @ brief compute min and max norms of vector in the Block
*/
template< class Matrix >
std::pair<P,P> get_min_max_norm(Matrix &A)
std::pair<P,P> get_min_max_norm(const Matrix &A) const
{
std::pair<P, P> mm;
mm.first = std::numeric_limits<P>::max();
......@@ -204,7 +204,7 @@ public:
/* \brief compute min and max (euclidian) norms of vector in the Block
*/
std::pair<P,P> get_min_max_norm()
std::pair<P,P> get_min_max_norm() const
{
std::pair<P,P> mm;
mm.first = std::numeric_limits<P>::max();
......@@ -307,7 +307,7 @@ public:
* otherwise the method shall return false;
*/
std::tuple<P,P,bool>
check_precision(const std::vector<P> &epsilon)
check_precision(const std::vector<P> &epsilon) const
{
if (epsilon.size() > 1) {
bool convergence = true;
......
......@@ -26,7 +26,6 @@ set(FABULOUS_XX_LIB_SRC
HessExtended.hpp
HessQR.hpp
HessStandard.hpp
HessStorage.hpp
krylov_algo.hpp
LapackInterface.hpp
Logger.hpp
......
......@@ -14,7 +14,6 @@ template <class> class ClassicRestarting;
#include "EigenUtils.hpp"
#include "Base.hpp"
#include "Timer.hpp"
#include "HessStorage.hpp"
namespace fabulous {
......@@ -36,7 +35,7 @@ public:
Base<S> &get_base() { return _base; }
template< class HESS > void save_hess(HESS &) const {}
template< class HESS > void restore_hess(HESS &,int,int) const {}
void store_LS(Block<S> &) const {}
void save_LS(Block<S> &) const {}
template<class ... Params> int run(Params ... ) const { return -0x0DEAD0; }
};
......@@ -53,57 +52,57 @@ private:
int _k; /*!<Number of eigen pair to keep*/
S _target; /*!< targeted eigen value*/
Base<S> &_base; /*!< Last Base generated */
HessStorage<S> &_hess; /*!< Last Hessenberg generated*/
HessStandard<S> _hess; /*!< Last Hessenberg generated*/
Block<S> _LS; /*!< Least square residual*/
public:
DeflatedRestarting(int k, Base<S> &base, HessStorage<S> &hess):
DeflatedRestarting(k, S{0.0}, hess, base)
DeflatedRestarting(int k, Base<S> &base):
DeflatedRestarting(k, S{0.0}, base)
{
}
DeflatedRestarting(int k, S target, Base<S> &base, HessStorage<S> &hess):
DeflatedRestarting(int k, S target, Base<S> &base):
_k{k}, _target{target},
_base{base}, _hess{hess}, _LS{}
_base{base}, _LS{}
{
}
operator bool () { return true; }
Base<S> &get_base() { return _base; }
void store_LS(Block<S> &LS) { _LS = LS; }
template< class HESS > void save_hess(HESS &hess) { _hess.store_hess(hess); }
template< class HESS > void restore_hess(HESS &hess, int k, int nbRHS)
{
for (int i=0; i<k; ++i)
memcpy( hess.get_vect(i), _hess.get_vect(i), (nbRHS+k)*sizeof(S) );
hess.reserve_DR(k);
}
void save_LS(Block<S> &LS) { _LS = LS; }
template< class HESS > void save_hess(HESS &hess) { _hess = hess; }
/**
* \brief This function handle the deflation of eigenvalues at restart.
*
* \param[out] R1 the rhs block for least square problems
* \param[out] hess the new hessenberg to be filled with H1new
*/
int run(Block<S> &R1)
template<class HESS >
int run(Block<S> &R1, HESS &hess)
{
double tic = Timer::get_time();
int p = _LS.get_nb_col();
int nm = _base.get_nb_vect() - p;
int dim = _base.get_nb_row();
std::cout<<"About to DR ################## 2 !!!\n"
assert( hess.get_ptr() != _hess.get_ptr() ) ;
std::cout<<"About to DR ################## !!!\n"
<<"Cste ::\n"
<<"\tk :: "<<_k<<"\n"
<<"\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"
<<"\t_Base nbCols :: "<<_base.get_nb_vect()<<"\n"
<<"\tldh :: "<<_hess.get_leading_dim()<<"\n"
<<"\tdim :: "<<dim<<"\n"
<<"\tLS : leading dim : "<<_LS.get_nb_row()<<"\n"
<<"\tLS : leading dim : "<<_LS.get_leading_dim()<<"\n"
<<"\tLS : nbrow : "<<_LS.get_nb_row()<<"\n"
<<std::endl;
Block<S> A{nm, nm};
......@@ -119,13 +118,13 @@ public:
//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(i, j));
B.at(i,j) = std::conj(_hess.at(j, i));
Block<S> vleft{1, 0}; // hold left eigen vectors that won't be computed
Block<S> vright{nm, nm}; // hold right eigen vectors before sorting them
/* complex alpha in order to store the eigen value that
* might be complex even in real case */
::fabulous::Block<std::complex<P>> alpha{nm, 1};
Block<std::complex<P>> alpha{nm, 1};
Block<S> beta{nm, 1};
int err = LapackKernI::ggev(
......@@ -153,8 +152,11 @@ public:
Block<S> Gleft = G.sub_block(0, 0, nm, newK);
Block<S> Gright = G.sub_block(0, newK, nm+p, p);
Gleft.copy(vright); // copy eigen vector in leftpart
Gright.copy(_LS); // Copy PR_m (last projected residual) inside right part of G
for (int j = 0; j < newK; ++j) {
//NOTGOOD: Gleft.copy(vright); // copy eigen vector in leftpart
memcpy(Gleft.get_vect(j), vright.get_vect(ind[j]), sizeof(S)*nm); // BETTER
}
R1.realloc(p+newK, p); // Init R1
......@@ -163,25 +165,23 @@ public:
Q.InPlaceQRFacto(R);
// Copy last columns of R inside R1
Block<S> Rp = R.sub_block(newK, newK, p, p);
Block<S> Rp = R.sub_block(0, newK, newK+p, p);
R1.copy(Rp);
Block<S> V1new{dim, newK};
Block<S> V2new{dim, p};
Block<S> Qa = Q.sub_block(0, 0, nm, newK);
Block<S> Qb = Q.sub_block(0, newK, nm+p, p);
LapackKernI::gemm( // V1 new = Vm * Qg{1:nm,1:k}
dim, newK, nm,
_base.get_ptr(), _base.get_leading_dim(),
Qa.get_ptr(), Qa.get_leading_dim(),
Q.get_ptr(), Q.get_leading_dim(),
V1new.get_ptr(), V1new.get_leading_dim(),
S{1.0}, S{0.0}
);
LapackKernI::gemm(// V2 new = Vm+1 * Qg{1:nm+p,k+1:k+p}
dim, p, nm+p,
_base.get_ptr(), _base.get_leading_dim(),
Qb.get_ptr(), Qb.get_leading_dim(),
Q.get_vect(newK), Q.get_leading_dim(),
V2new.get_ptr(), V2new.get_leading_dim(),
S{1.0}, S{0.0}
);
......@@ -193,15 +193,15 @@ public:
LapackKernI::gemm(
nm+p, newK, nm,
_hess.get_ptr(), _hess.get_leading_dim(),
Q.get_ptr(), Q.get_nb_row(),
Q.get_ptr(), Q.get_leading_dim(),
tmp_HmXQg.get_ptr(), tmp_HmXQg.get_leading_dim(),
S{1.0}, S{0.0}
);
LapackKernI::Tgemm(
newK+p, newK, nm+p,
Q.get_ptr(), Q.get_nb_row(),
Q.get_ptr(), Q.get_leading_dim(),
tmp_HmXQg.get_ptr(), tmp_HmXQg.get_leading_dim(),
H1new.get_ptr(), H1new.get_nb_row(),
H1new.get_ptr(), H1new.get_leading_dim(),
S{1.0}, S{0.0}
);
}
......@@ -218,12 +218,15 @@ public:
_base.increase(V2new.get_nb_col());
// Modify Hess ...
_hess.init_back(H1new);
hess.zero();
hess.copy(H1new);
hess.reserve_DR(H1new.get_nb_col());
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;
return newK;
}
};
......
......@@ -4,7 +4,6 @@
#include <complex>
#include <type_traits>
#include "HessStorage.hpp"
#include "LapackInterface.hpp"
namespace fabulous {
......@@ -14,6 +13,8 @@ namespace fabulous {
*
* \param k number of eigen pair wanted
* \param target the targeted eigen value
* \param[in] alpha as returned by the GGEV Lapack Kernel
* \param[in] beta as returned by the GGEV Lapack Kernel
* \param[out] ind the indices of the found vectors
*
* \return the indexes of the k nearest ones to target
......@@ -174,7 +175,7 @@ void CheckEigenResults(Block<S>& A, Block<S>& B, Block<S>& vr,
S{1.0}, S{0.0}
);
LapackKernI::gemm( // res_i = -res_i + (alpha[i]/beta[i])B * vr_i
LapackKernI::gemm( // res_i = (alpha[i]/beta[i])B * vr_i - res_i
dim, 1, dim,
B.get_vect(), B.get_leading_dim(),
vr.get_vect(i), vr.get_leading_dim(),
......
......@@ -29,15 +29,14 @@ public:
FABULOUS_INHERITS_BLOCK(S);
private:
int _nb_eigen_vector;
int _nbRHS; /**< number of right hand sides */
int _nb_vect; /**< number of vector */
Block<S> _R1; /**< right hand sides of projected problem */
public:
HessStandard() = default;
HessStandard(int max_krylov_space_size, int nbRHS):
super{ max_krylov_space_size+nbRHS, max_krylov_space_size},
_nb_eigen_vector{0},
_nbRHS{nbRHS},
_nb_vect{0},
_R1{}
......@@ -73,7 +72,6 @@ public:
void reserve_DR(int nb_eigen_vector)
{
_nb_eigen_vector = nb_eigen_vector;
_nb_vect = nb_eigen_vector;
}
......@@ -97,10 +95,13 @@ public:
int M = _nb_vect + _nbRHS;
int N = _nb_vect;
Y.copy(_R1);
FABULOUS_DEBUG("R1.nb_row()="<<_R1.get_nb_row());
Block<S> H{M, N};
H.copy(*this);
assert( Y.get_nb_row() == M );
FABULOUS_DEBUG("M="<<M<<" N="<<N);
int err = LapackKernI::gels(
M, N, _nbRHS,
H.get_ptr(), H.get_leading_dim(),
......@@ -115,7 +116,13 @@ public:
{
int M = _nb_vect + _nbRHS;
int N = _nb_vect;
assert( Y.get_nb_row() == N );
assert( PR.get_nb_row() == M );
PR.copy(_R1);
assert(PR.get_nb_col() == _R1.get_nb_col());
assert(PR.get_nb_row() >= _R1.get_nb_row());
FABULOUS_DEBUG("R1.nb_row()="<<_R1.get_nb_row());
LapackKernI::gemm(
M, _nbRHS, N,
......
#ifndef FABULOUS_HESS_STORAGE_HPP
#define FABULOUS_HESS_STORAGE_HPP
namespace fabulous {
/**
* \brief Save the hessenberg for restarting purpose
*
* This class will receive ownership of Hessenberg Matrix at the end of
* Arnoldi procedure, in order to use it for restarting purposes.
*/
template< class S > class HessStorage : public Block<S>
{
public:
FABULOUS_INHERITS_BLOCK(S);
HessStorage() = default;
template< class HESS >
void store_hess(HESS &hess)
{
int N = hess.get_nb_vect();
int M = hess.get_nb_hess_line();
Block<S> sub = hess.sub_block(0, 0, M, N);
*(Block<S>*)this = sub;
std::cout<<"Init Storage-Hessenberg with :\n"
<<"\t nbLine :: "<<M<<"\n"
<<"\t nbCols ::"<< N<<"\n"
<<"\t Nb tot values :: "<< M*N<<"\n";
}
void init_back(Block<S>& H1new)
{
this->zero();
this->copy(H1new);
}
};
}; // namespace fabulous
#endif // FABULOUS_HESS_STORAGE_HPP
......@@ -15,7 +15,7 @@ namespace fabulous {
typedef KrylovAlgo<TYPE__> super;
#define FABULOUS_INHERITS_KRYLOV_ALGO_PROTECTED \
using super::_krylov_space_spanned; \