Commit 6bd09317 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

fix a bug in DR due to sign problem when computing Projected residual

parent a1cce305
......@@ -88,8 +88,8 @@ else()
endif()
if(FABULOUS_BUILD_WARNINGS)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fdiagnostics-color=always -Wall -W -Wextra")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fdiagnostics-color=always -Wall -W -Wextra")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fdiagnostics-color=always -Wall -Wextra")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fdiagnostics-color=always -Wall -Wextra")
endif()
if(BUILD_SHARED_LIBS)
......@@ -107,14 +107,11 @@ if(BUILD_SHARED_LIBS)
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
endif()
#Ajout d'un repertoire src
add_subdirectory(src)
#Ajout du meme repertoire pour les includes des tests
include_directories(src)
#On vire les ; qui viennent des findBlas lapack...
list(REMOVE_DUPLICATES CMAKE_EXE_LINKER_FLAGS)
string(REPLACE ";" " " CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}")
......
......@@ -84,7 +84,6 @@ public:
};
}; // namespace fabulous
#endif // FABULOUS_ALGORITHM_H
......@@ -24,16 +24,6 @@ namespace fabulous {
*/
template< class S > class Arnoldi
{
public:
/**
* \brief information about a finished Arnoldi procedure
*/
struct ArnReturn {
int size_krylov_space; /*!< Size Krylov Space reached */
int nb_iteration; /*!< Number of iterations done */
bool has_converged; /*!< Flag set to true if convergence has been reached. */
};
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
......@@ -48,9 +38,9 @@ private:
o << "######## Mat Vect product scheduled "<< nb_mvp <<"###########\n";
}
void print_footer(int iteration_count, std::ostream &o = std::cout)
void print_footer(std::ostream &o = std::cout)
{
o << "################# Iterations done: "<<iteration_count<<"(+1)\n";
o << "################# Iterations done: "<<_iteration_count<<"(+1)\n";
o << "#######################################################\n";
}
......@@ -61,11 +51,21 @@ private:
Block<S> _solution;
Block<S> _lastY;
Base<S> &_base;
int _dim;
int _nbRHS;
int _size_krylov_space;
int _nb_mvp;
int _iteration_count;
// DR SPECIFIC
bool _cold_restart;
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:
template< class Matrix, class Block, class Base >
void X_plus_MVY(Matrix &A, Block &X, Base &V, const Block &Y)
......@@ -76,11 +76,10 @@ private:
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
A.PrecondBlockVect(VY, MVY); // MVY := M*V*Y
} else {
V.compute_product(Y, MVY); // MVY = V*Y
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);
......@@ -90,7 +89,8 @@ private:
void compute_real_residual(Matrix &A, const 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
A.MatBlockVect(X, R, S{-1.0}, S{1.0}); // R <- B-A*X
_nb_mvp += X.get_nb_col();
}
template< class Matrix, class Block >
......@@ -99,21 +99,11 @@ private:
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
A.MatBlockVect(MV, W); // W = A*MV
} else {
A.MatBlockVect(V, W); // W = A * V
}
}
template<class Matrix, class Block>
void compute_solution(Matrix &A, Block &X)
{
int nb_vect = _hess.get_nb_vect();
Block YY{nb_vect+_nbRHS, _nbRHS};
_hess.compute_least_square_sol(YY);
Block Y = YY.sub_block(0, 0, nb_vect, _nbRHS);
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
_nb_mvp += V.get_nb_col();
}
template<class Matrix, class Block>
......@@ -163,6 +153,8 @@ public:
_dim(dim),
_nbRHS(nbRHS),
_size_krylov_space{0},
_nb_mvp{0},
_iteration_count{0},
_cold_restart{false}
{
}
......@@ -186,23 +178,22 @@ public:
* \return information about the iterations
*/
template<class Matrix, class Block, class Restarter>
ArnReturn 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::ICGS,
OrthoType ortho_type = OrthoType::RUHE )
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::ICGS,
OrthoType ortho_type = OrthoType::RUHE )
{
print_header(max_krylov_space_size);
Block R1{}; // R1 = RHS of projected problem
int iteration_count = 0;
_logger.notify_iteration_begin(iteration_count, _size_krylov_space);
_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); _base.increase(_nbRHS);
compute_real_residual(A, X, B, R0);
R1.realloc(_nbRHS, _nbRHS);
......@@ -212,23 +203,23 @@ public:
bool convergence = std::get<2>(MinMaxConv);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_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);
++ _iteration_count;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
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("Wcandidate");
_base.increase(_nbRHS); // add W to base
_base.increase(_nbRHS); // add W to base
_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);
......@@ -240,29 +231,24 @@ public:
Block PR{_hess.get_nb_vect()+_nbRHS, _nbRHS};
_hess.compute_proj_residual(Y, 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);
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, Y);
}
if (iteration_count) {
if (_solution_computed)
X.copy(_solution);
else
compute_solution(A, X, _lastY);
}
if (_solution_computed)
X.copy(_solution);
else if (_iteration_count)
compute_solution(A, X, _lastY);
if (_cold_restart)
_base.reset();
print_footer(iteration_count);
print_footer();
restarter.save_hess(_hess);
return ArnReturn{_size_krylov_space, iteration_count, convergence};
return convergence;
}
};
......
......@@ -19,15 +19,6 @@ namespace fabulous {
*/
template <class S> class Arnoldi_IB
{
public:
/**
* \brief information about a finished Arnoldi_IB procedure
*/
struct ArnReturn {
int size_krylov_space; /*!< Size Krylov Space reached */
int nb_iteration; /*!< Number of iterations done */
bool has_converged; /*!< Flag set to true if convergence has been reached. */
};
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
......@@ -56,8 +47,10 @@ private:
Block<S> _last_Y;
Base<S> &_base;
// GENERIC
int _dim;
int _nbRHS;
int _nb_mvp;
int _size_krylov_space;
int _iteration_count;
......@@ -71,6 +64,12 @@ private:
// DR specific
bool _cold_restart;
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:
template< class Matrix, class Block, class Base >
void X_plus_MVY(Matrix &A, Block &X, Base &V, Block &Y)
......@@ -109,6 +108,7 @@ private:
} else {
A.MatBlockVect(V, W); // W = A * V
}
_nb_mvp += V.get_nb_col();
}
template<class Matrix, class Block>
......@@ -116,6 +116,7 @@ private:
{
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();
}
template<class Matrix, class Block>
......@@ -321,6 +322,12 @@ public:
_base{restarter.get_base()},
_dim(dim),
_nbRHS(nbRHS),
_nb_mvp{0},
_size_krylov_space{0},
_iteration_count{0},
_nb_direction_kept{nbRHS},
_old_nb_direction_kept{nbRHS},
_nb_direction_discarded{0},
_IB_happened{false},
_cold_restart{false}
{
......@@ -346,25 +353,20 @@ public:
*
*/
template< class Matrix, class Block, class DRParam >
ArnReturn run(
Matrix &A, Block &X, Block &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
DRParam &DR_struct,
OrthoScheme ortho_scheme = OrthoScheme::ICGS,
OrthoType ortho_type = OrthoType::RUHE)
bool run( Matrix &A, Block &X, Block &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
DRParam &DR_struct,
OrthoScheme ortho_scheme = OrthoScheme::ICGS,
OrthoType ortho_type = OrthoType::RUHE)
{
(void) DR_struct;
bool convergence = false;
print_header(max_krylov_space_size);
_base.reset(); // DR not implemented ; we force a cold restart
_nbRHS = B.get_nb_col(); // Initial Size of block
_dim = B.get_nb_row();
_size_krylov_space = 0;
_iteration_count = 0;
_nb_direction_kept = _nbRHS;
_old_nb_direction_kept = _nbRHS;
_nb_direction_discarded = 0;
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 );
auto inv_epsilon = array_inverse(epsilon);
......@@ -382,7 +384,7 @@ public:
_logger.log_real_residual(A, B, X);
if (_nb_direction_kept == 0)
return ArnReturn{0, 0, true};
return true;
// Should not be true (because we already tested it):
convergence = (_nb_direction_discarded == _nbRHS);
......@@ -442,7 +444,7 @@ public:
#endif
print_footer();
return ArnReturn{_size_krylov_space, _iteration_count, convergence};
return convergence;
}
};
......
......@@ -108,7 +108,7 @@ void ArnoldiBlock_MGS(Base &base, Block &H, Block &W)
Block H_k = H.sub_block(size_block_sum, 0, size_block, W_size);
size_block_sum += size_block;
LapackKernI::Tgemm( // H_{ij} += V_i^{h} * W
LapackKernI::Tgemm( // H_{ij} = V_i^{h} * W
size_block, W_size, dim,
base.get_block_ptr(k), base.get_leading_dim(),
W.get_ptr(), W.get_leading_dim(),
......@@ -145,14 +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());
//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);
//FABULOUS_DEBUG("base.block["<<k<<"].size="<<size_block);
LapackKernI::Tgemm( // tmp = V_i^{h} * W
size_block, W_size, dim,
......
......@@ -25,16 +25,6 @@ namespace fabulous {
*/
template<class S> class Arnoldi_QRInc {
public:
/**
* \brief information about a finished Arnoldi_QRInc procedure
*/
struct ArnReturn {
int size_krylov_space; /*!< Size Krylov Space reached */
int nb_iteration; /*!< Number of iterations done */
bool has_converged; /*!< Flag set to true if convergence has been reached. */
};
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
......@@ -45,8 +35,17 @@ private:
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:
......@@ -57,13 +56,12 @@ private:
o << "######## Mat Vect product scheduled "<< nb_mvp <<"###########\n";
}
void print_footer(int iteration_count, std::ostream &o = std::cout)
void print_footer(std::ostream &o = std::cout)
{
o << "################# Iterations done: "<<iteration_count<<"(+1)\n";
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
{
......@@ -74,6 +72,7 @@ private:
} else {
A.MatBlockVect(V, W); // W = A * V
}
_nb_mvp += V.get_nb_col();
}
// X_n = X_0 + M * V * Y
......@@ -133,6 +132,7 @@ private:
{
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:
......@@ -144,7 +144,10 @@ public:
_solution_computed(false),
_base{restarter.get_base()},
_dim(dim),
_nbRHS(nbRHS)
_nbRHS(nbRHS),
_size_krylov_space{0},
_nb_mvp{0},
_iteration_count{0}
{
}
......@@ -167,43 +170,41 @@ public:
* \return Struct containing information about the procedure
*/
template<class Matrix, class Block, class Restarter>
ArnReturn run(
const Matrix &A, Block &X, const Block &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter&,
OrthoScheme ortho_scheme = OrthoScheme::MGS,
OrthoType ortho_type = OrthoType::BLOCK)
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 = OrthoScheme::MGS,
OrthoType ortho_type = OrthoType::BLOCK)
{
const int nbRHS = B.get_nb_col();
int iteration_count = 0, size_krylov_space = 0;
assert( _nbRHS == B.get_nb_col() );
_base.reset(); // DR not implemented ; force a cold restart
_logger.notify_iteration_begin(iteration_count, size_krylov_space);
Block R0 = _base.get_W(nbRHS);
_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};
Block Lambda1{_nbRHS, _nbRHS};
A.QRFacto(R0, Lambda1); // R0.InPlaceQRFacto_UserDot(R1, A);
_base.increase(nbRHS);
_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();
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
iteration_count, size_krylov_space, nbRHS, MinMaxConv);
_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);
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);
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();
_base.increase(_nbRHS);
_size_krylov_space = _base.get_nb_vect();
_hess.factorize_last_column();
auto MinMaxConv = _hess.check_convergence(epsilon);
......@@ -212,14 +213,14 @@ public:
proj_convergence, A, X, B, epsilon);
_logger.notify_iteration_end(
iteration_count, size_krylov_space, nbRHS, MinMaxConv);
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
//check_arnoldi_formula(A, _base, hess);
}
FABULOUS_DEBUG("loop end: "<<iteration_count);
FABULOUS_DEBUG("loop end: "<<_iteration_count);
if (_solution_computed)
X.copy(_solution);
else if (iteration_count)
else if (_iteration_count)
compute_solution(A, X);
#ifdef FABULOUS_DEBUG_MODE
......@@ -229,9 +230,8 @@ public:
}
#endif
print_footer(iteration_count);
ArnReturn returned_struct{size_krylov_space, iteration_count, convergence};
return returned_struct;
print_footer();
return convergence;
}
};
......
......@@ -88,14 +88,13 @@ public:
Base<S> base{dim, size_to_span+nbRHS};
ClassicRestarting<S> restarter{base};
ARNOLDI<S> arnoldi{_logger, restarter, dim, nbRHS, size_to_span};
auto r = arnoldi.run( A, X, B, size_to_span, epsilon, restarter,
ortho_scheme, ortho_type );
finished = arnoldi.run( A, X, B, size_to_span, epsilon, restarter,
ortho_scheme, ortho_type );
finished = r.has_converged;
_nb_iteration += r.nb_iteration;
_mvp += r.size_krylov_space;
_nb_iteration += arnoldi.get_nb_iteration();
_mvp += arnoldi.get_nb_mvp();
++ _nb_restart;
print_iteration_end_info(r.size_krylov_space);
print_iteration_end_info(arnoldi.get_krylov_space_size());
}
X.scale(normB); B.scale(normB);
return _mvp;
......
......@@ -117,14 +117,13 @@ public:
}
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,
ortho_scheme, ortho_type );
finished = arnoldi.run( A, X, B, size_to_span, epsilon, restarter,
ortho_scheme, ortho_type );
finished = r.has_converged;
_nb_iteration += r.nb_iteration;
_mvp += r.size_krylov_space;
_nb_iteration += arnoldi.get_nb_iteration();
_mvp += arnoldi.get_nb_mvp();
++ _nb_restart;
print_iteration_end_info(r.size_krylov_space);
print_iteration_end_info(arnoldi.get_krylov_space_size());
}
X.scale(normB); B.scale(normB);
return _mvp;
......
......@@ -36,7 +36,9 @@ public:
Block<S>{dim, max_krylov_space_size},
_nb_block{0},
_nb_vect{0},
_max_krylov_space_size{max_krylov_space_size}
_max_krylov_space_size{max_krylov_space_size},
_block_sizes{},
_block_sizes_sum{}
{
_block_sizes_sum.emplace_back(0);
}
......
......@@ -120,7 +120,7 @@ public:
for (int i = 0; i < nm; ++i)
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> vleft{}; // 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 */
......@@ -149,28 +149,24 @@ public:
// Create block to hold ~G
Block<S> G{nm+p, newK+p};
Block<S> Gleft = G.sub_block(0, 0, nm, newK);
Block<S> Gright = G.sub_block(0, newK, nm+p, p);
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
Gright.copy(_LS); // Copy last residual inside right part of G
// Copy back the corresponding eigen vectors inside left part G.
Block<S> Gleft = G.sub_block(0, 0, nm, newK);
for (int j = 0; j < newK; ++j)
memcpy(Gleft.get_vect(j), vright.get_vect(ind[j]), nm*sizeof(S));
Block<S> Q = G;
Block<S> R{newK+p, newK+p};
Q.InPlaceQRFacto(R);
// Copy last columns of R inside R1
Block<S> Rp = R.sub_block(0, newK, newK+p, p);
R1.copy(Rp);
Block<S> Rleft = R.sub_block(0, newK, newK+p, p);
R1.realloc(p+newK, p); // Init R1
R1.copy(Rleft);
Block<S> V1new{dim, newK};
Block<S> V2new{dim, p};
LapackKernI::gemm( // V1 new = Vm * Qg{1:nm,1:k}
dim, newK, nm,
_base.get_ptr(), _base.get_leading_dim(),
......
......@@ -31,6 +31,7 @@ public:
private:
int _nbRHS; /**< number of right hand sides */
int _nb_vect; /**< number of vector */
int _nb_eigen_pair;
Block<S> _R1; /**< right hand sides of projected problem */
public:
......@@ -39,6 +40,7 @@ public:
super{ max_krylov_space_size+nbRHS, max_krylov_space_size},
_nbRHS{nbRHS},
_nb_vect{0},
_nb_eigen_pair{0},
_R1{}
{
}
......@@ -73,17 +75,26 @@ public:
void reserve_DR(int nb_eigen_vector)
{
_nb_vect = nb_eigen_vector;
_nb_eigen_pair = nb_eigen_vector;
}
void display(std::string name = "")
{
Block<S> sub = this->sub_block(0, 0, _nb_vect+_nbRHS, _nb_vect);
sub.display(name);
}
void display_hess_bitmap(std::string name = "")
{
int M = _nb_vect + _nbRHS;
int N = _nb_vect;
std::cout << name<<"\n";
std::cout << name<<" (M="<<M<<"; N="<<N<<")\n";
for (int i = 0; i < M; ++i) {
for (int j = 0; j < N; ++j) {
bool b = at(i, j) != S{0.0};
std::cout << (b ? ". " : " ");
if (i < _nb_eigen_pair || j < _nb_eigen_pair)
std::cout << Color::red;
std::cout << (b ? ". " : " ") << Color::reset;
}
std::cout<<"\n";
}
......@@ -95,13 +106,13 @@ public:
int M = _nb_vect + _nbRHS;
int N = _nb_vect;
Y.copy(_R1);
FABULOUS_DEBUG("R1.nb_row()="<<_R1.get_nb_row());
//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);
//FABULOUS_DEBUG("M="<<M<<" N="<<N);
int err = LapackKernI::gels(
M, N, _nbRHS,
H.get_ptr(), H.get_leading_dim(),
......@@ -121,15 +132,14 @@ public:
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());
//FABULOUS_DEBUG("R1.nb_row()="<<_R1.get_nb_row());
LapackKernI::gemm(
M, _nbRHS, N,
get_ptr(), get_leading_dim(),
Y.get_ptr(), Y.get_leading_dim(),
PR.get_ptr(), PR.get_leading_dim(),
S{1.0}, S{-1.0}
S{-1.0}, S{1.0}
);
}
};
......
......@@ -15,7 +15,7 @@ namespace fabulous {