Commit 74f3645f authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

update error system

parent 3ddc6337
......@@ -46,3 +46,4 @@
* TODO error handling system
* TODO report convergence in logger
* TODO rework C++ api for multiple algorithm
* TODO review chameleon sub descriptor solution
......@@ -14,7 +14,6 @@ template<class S> class BCG;
#include "fabulous/data/Base.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
......@@ -136,9 +135,7 @@ public:
* \param[in] B Block containing Right Hand Side
* \param algo instance of Algorithm as returned by ::fabulous::bcg::std()
* \param max_mvp Maximum number of Matrix Vector Product.
* \param max_krylov_space_size maximum size of Krylov space
* \param[in] epsilon Target accuracy
* \param ortho OrthoParam
* (through QR factorization) or vector wise arnoldi. <br/>
* (In distributed, only the vector wise will works)
*
......@@ -157,13 +154,14 @@ public:
print_start_info( dim, nbRHS, max_mvp, epsilon, A, algo );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
normB = compute_norms(B, A);
inv_normB = array_inverse(normB);
X.scale(inv_normB); B.scale(inv_normB);
scale(X, inv_normB);
scale(B, inv_normB);
using BcgRunner = typename Algo::template algo<S>;
BcgRunner bcg{_logger, dim, nbRHS, max_mvp};
bool convergence = bcg.run( A, X, B, epsilon );
bcg.run( A, X, B, epsilon );
_nb_iteration += bcg.get_nb_iteration();
_mvp = bcg.get_nb_mvp();
print_end_info();
......@@ -174,9 +172,9 @@ public:
"You may want to check the algorithm parameters"
);
}
(void) /*FIXME */ convergence;
X.scale(normB); B.scale(normB);
scale(X, normB);
scale(B, normB);
return _mvp;
}
......
......@@ -129,7 +129,7 @@ public:
int solve( const Matrix &A,
int dim, int nbRHS, S *B_, int ldb, S *X_, int ldx,
const std::vector<P> &epsilon, Algo algo,
const int max_mvp, const int max_space, OrthoParam ortho)
const int max_mvp, const int max_space, OrthoParam ortho )
{
Block<S> X{dim, nbRHS, ldx, X_};
Block<S> B{dim, nbRHS, ldb, B_};
......@@ -175,9 +175,10 @@ public:
print_start_info( dim, nbRHS, max_mvp, max_krylov_space_size,
ortho, epsilon, A, algo );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
normB = compute_norms(B, A);
inv_normB = array_inverse(normB);
X.scale(inv_normB); B.scale(inv_normB);
scale(X, inv_normB);
scale(B, inv_normB);
Base<S> V{dim, max_krylov_space_size+nbRHS};
Base<S> Z{dim, max_krylov_space_size+nbRHS};
......@@ -211,7 +212,8 @@ public:
);
}
X.scale(normB); B.scale(normB);
scale(X, normB);
scale(B, normB);
return _mvp;
}
......
......@@ -196,9 +196,10 @@ public:
print_start_info( dim, nbRHS, max_mvp, max_krylov_space_size,
ortho, restart_param, epsilon, A, algo );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
normB = compute_norms(B, A);
inv_normB = array_inverse(normB);
X.scale(inv_normB); B.scale(inv_normB);
scale(X, inv_normB);
scale(B, inv_normB);
const int nb_eigen_pair = restart_param.get_k();
Base<S> base{dim, max_krylov_space_size+nbRHS};
......@@ -233,7 +234,8 @@ public:
);
}
X.scale(normB); B.scale(normB);
scale(X, normB);
scale(B, normB);
return _mvp;
}
......
......@@ -373,6 +373,8 @@ public:
get_everyone(XJ, X, J);
}
print_footer();
_logger.set_convergence(_nb_kept_direction == 0);
return (_nb_kept_direction == 0);
}
};
......
......@@ -124,8 +124,8 @@ private:
P_.get_ptr(), P_.get_leading_dim(),
jpvt, tau);
(void) A; /* FIXME :: this GEQPF must be replaced with distributed-friendly QPF */
if ( err != 0) {
FABULOUS_FATAL_ERROR("geqp3 failed with error code='"<<err<<"'");
if ( err != 0 ) {
FABULOUS_THROW(Kernel, "geqp3 failed with error code='"<<err<<"'");
}
int nb_kept_direction = 0;
......@@ -148,7 +148,7 @@ private:
err = lapacke::orgqr(_dim, _nbRHS, nb_kept_direction, P_.get_ptr(), P_.get_leading_dim(), tau.data());
if ( err != 0 ) {
FABULOUS_FATAL_ERROR("orgqr failed with error code='"<<err<<"'");
FABULOUS_THROW(Kernel, "orgqr failed with error code='"<<err<<"'");
}
int old_nb_kept_direction = _nb_kept_direction;
......@@ -238,14 +238,14 @@ private:
err = lapacke::potrf('L', _nb_kept_direction,
PtQ.get_ptr(), PtQ.get_leading_dim());
if (err != 0) {
FABULOUS_FATAL_ERROR("potrf failed with error code='"<<err<<"'");
FABULOUS_THROW(Kernel, "potrf failed with error code='"<<err<<"'");
}
err = lapacke::potrs('L', _nb_kept_direction, _nbRHS,
PtQ.get_ptr(), PtQ.get_leading_dim(),
alpha.get_ptr(), alpha.get_leading_dim());
if (err != 0) {
FABULOUS_FATAL_ERROR("potrs failed with error code='"<<err<<"'");
FABULOUS_THROW(Kernel, "potrs failed with error code='"<<err<<"'");
}
}
......@@ -264,7 +264,7 @@ private:
PtQ.get_ptr(), PtQ.get_leading_dim(),
beta.get_ptr(), beta.get_leading_dim());
if ( err != 0 ) {
FABULOUS_FATAL_ERROR("potrs failed with error code='"<<err<<"'");
FABULOUS_THROW(Kernel, "potrs failed with error code='"<<err<<"'");
}
}
......@@ -309,7 +309,6 @@ private:
S{-1.0}, S{1.0} );
}
public:
Arnoldi(Logger<P> &logger, int dim, int nbRHS, int max_mvp):
_logger{logger},
......@@ -381,6 +380,7 @@ public:
}
#endif
print_footer();
_logger.set_convergence(_nb_kept_direction == 0);
return (_nb_kept_direction == 0);
}
};
......
......@@ -288,6 +288,7 @@ public:
}
print_footer();
_logger.set_convergence(convergence);
return convergence;
}
......
......@@ -251,7 +251,7 @@ private:
compute_real_residual(A, X, B, Q0);
A.QRFacto(Q0, Lambda0); L0.copy(Lambda0);
_nb_direction_kept = L0.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_kept = decomposition_svd(L0, U, epsilon, inv_epsilon);
_old_nb_direction_kept = _nb_direction_kept;
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
......@@ -338,7 +338,7 @@ private:
_F.compute_proj_residual(PR, Y);
Block<S> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_kept = decomposition_svd(PR, U, epsilon, inv_epsilon);
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
_U = U;
......@@ -349,7 +349,8 @@ private:
_cold_restart = true;
} else {
// Even more suspicious
FABULOUS_FATAL_ERROR(
FABULOUS_THROW(
Numeric,
"Directions kepts during SVD of projected residual is "
" equal to 0, but residual itself has not converged. "
"That should not be possible."
......@@ -462,6 +463,7 @@ public:
#endif
print_footer();
_logger.set_convergence(convergence);
return convergence;
}
......
......@@ -266,7 +266,7 @@ private:
compute_real_residual(A, X, B, Q0);
A.QRFacto(Q0, Lambda0); L0.copy(Lambda0);
_nb_direction_kept = L0.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_kept = decomposition_svd(L0, U, epsilon, inv_epsilon);
_old_nb_direction_kept = _nb_direction_kept;
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
......@@ -356,7 +356,7 @@ private:
_logger.notify_ib_begin();
_F.compute_proj_residual(PR, Y);
Block<S> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_kept = decomposition_svd(PR, U, epsilon, inv_epsilon);
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
_U = U;
......@@ -370,7 +370,8 @@ private:
_cold_restart = true;
} else {
// Even more suspicious
FABULOUS_FATAL_ERROR(
FABULOUS_THROW(
Numeric,
"Number of directions kepts during SVD of projected "
"residual is equal to 0, but residual itself has not "
"converged. That should not be possible."
......@@ -494,6 +495,7 @@ public:
#endif
print_footer();
_logger.set_convergence(convergence);
return convergence;
}
......
......@@ -188,31 +188,6 @@ public:
return std::sqrt(fabulous::real(res));
}
/*! \brief compute norms of vector in the Block */
template< class Matrix >
std::vector<P> compute_norms(Matrix &A) const
{
std::vector<P> norms;
norms.resize(_n);
for (int j = 0; j < _n; ++j) {
norms[j] = get_norm(j, A);
if (norms[j] == 0.0) {
FABULOUS_WARNING("B["<<j<<"] (RHS column "<<j<<") have norm 0.0");
norms[j] = 1.0;
}
}
return norms;
}
/*! \brief scale vectors in the Block with given coefficients */
template<class U>
void scale(const std::vector<U> &coef)
{
for (int j = 0; j < _n; ++j)
if (!std::isnan(std::norm(coef[j])))
lapacke::scal(_m, S{coef[j]}, get_vect(j), 1);
}
/*! \brief compute min and max norms of vector in the Block */
template<class Matrix>
std::pair<P,P> get_min_max_norm(const Matrix &A) const
......@@ -399,83 +374,123 @@ public:
return std::make_tuple(mm.first, mm.second, mm.second < epsilon[0]);
}
/*! \brief Wrapper for the other DecompositionSVD method */
int decomposition_SVD(Block &U1,
const std::vector<P> &epsilon,
const std::vector<P> &inv_epsilon)
/*! \brief return number of flops performed by last operation */
int64_t get_last_flops() const
{
if (epsilon.size() == 1)
return decomposition_SVD(U1, epsilon.front());
scale(inv_epsilon);
return decomposition_SVD(U1, P{1.0});
return _last_flops;
}
/*!
* \brief Compute SVD of current block
*
* \f$ M = \mathbb{U}_1 \Sigma_1 \mathbb{V}_1 + \mathbb{U}_2 \Sigma \mathbb{V}_2 \f$
*
* \param[out] U1 vectors associated with singular values
* superior to threshold.
}; // end class Block
* \param threshold minimum norm of kept singular values
*
* \return count of singular values greater(in module)
* than the epsilon threshold
*/
int decomposition_SVD(Block &U1, P threshold)
{
std::vector<P> sigma, superb;
sigma.resize(_n);
superb.resize(_n-1);
Block U{_m, _n};
Block V{1, 0};
int err = lapacke::gesvd(
'S', /* only first _n vectors computes */
'N', /* V not computed */
_m, _n,
_ptr, _ld,
sigma.data(),
U.get_ptr(), U.get_leading_dim(),
V.get_ptr(), V.get_leading_dim(),
superb.data()
);
if (err != 0) {
FABULOUS_ERROR("Problem ... in SVD, did not converge err="<<err);
if (err > 0) {
for (int i=0; i < err; ++i)
std::cout << superb[i] << "\t";
std::cout<<"\n";
}
FABULOUS_FATAL_ERROR("SVD");
}
/*! \brief scale vectors in the Block with given coefficients */
template<class S, class P>
void scale(Block<S> &B, const std::vector<P> &coef)
{
int M = B.get_nb_row();
int N = B.get_nb_col();
// Count the number of Singular Value superior to threshold
int nb_krylov_direction = 0;
while (nb_krylov_direction < _n
&& sigma[nb_krylov_direction] > threshold) {
++ nb_krylov_direction;
}
for (int j = 0; j < N; ++j)
if (!std::isnan(std::norm(coef[j])))
lapacke::scal(M, S{coef[j]}, B.get_vect(j), 1);
}
// Test if U1 is already allocated
if (U1.get_nb_col() == 0)
U1.realloc(_m, nb_krylov_direction);
// Fill the U with vectors from U
U1.copy(U);
/*! \brief Wrapper for the other DecompositionSVD method */
template<class S, class P>
int decomposition_svd(Block<S> &B, Block<S> &U1,
const std::vector<P> &epsilon,
const std::vector<P> &inv_epsilon)
{
if ( epsilon.size() == 1 ) {
return decomposition_svd(B, U1, epsilon.front());
} else {
scale(B, inv_epsilon);
return decomposition_svd(B, U1, P{1.0});
}
}
return nb_krylov_direction;
/*!
* \brief Compute SVD of Block B
*
* \f$ M = \mathbb{U}_1 \Sigma_1 \mathbb{V}_1 + \mathbb{U}_2 \Sigma \mathbb{V}_2 \f$
*
* \param[in,out] B the block to be decomposed
* \param[out] U1 vectors associated with singular values
* superior to threshold.
*
* \param threshold minimum norm of kept singular values
*
* \return count of singular values greater(in module)
* than the epsilon threshold
*/
template<class S, class P>
int decomposition_svd(Block<S> &B, Block<S> &U1, P threshold)
{
std::vector<P> sigma, superb;
int M = B.get_nb_row();
int N = B.get_nb_col();
sigma.resize(N);
superb.resize(N-1);
Block<S> U{M, N};
Block<S> V{};
int err = lapacke::gesvd(
'S', /* only first _n vectors computes */
'N', /* V not computed */
M, N,
B.get_ptr(), B.get_leading_dim(),
sigma.data(),
U.get_ptr(), U.get_leading_dim(),
V.get_ptr(), V.get_leading_dim(),
superb.data()
);
if (err != 0) {
if (err > 0) {
for (int i=0; i < err; ++i)
std::cout << superb[i] << "\t";
std::cout<<"\n";
}
FABULOUS_THROW(Kernel, "SVD kernel err="<<err);
}
/*! \brief return number of flops performed by last operation */
int64_t get_last_flops() const
{
return _last_flops;
// Count the number of Singular Value superior to threshold
int nb_krylov_direction = 0;
while (nb_krylov_direction < N && sigma[nb_krylov_direction] > threshold) {
++ nb_krylov_direction;
}
}; // end class Block
// Test if U1 is already allocated
if (U1.get_nb_col() == 0)
U1.realloc(M, nb_krylov_direction);
// Fill the U with vectors from U
U1.copy(U);
return nb_krylov_direction;
}
/*! \brief compute norms of vector in the Block */
template<class Matrix, class S, class P = typename Arithmetik<S>::primary_type>
std::vector<P> compute_norms(const Block<S> &B, Matrix &A)
{
int N = B.get_nb_col();
std::vector<P> norms;
norms.resize(N);
for (int j = 0; j < N; ++j) {
norms[j] = B.get_norm(j, A);
if (norms[j] == 0.0) {
FABULOUS_WARNING("B["<<j<<"] (RHS column "<<j<<") have norm 0.0");
norms[j] = 1.0;
}
}
return norms;
}
template<class U>
std::ostream& operator<<(std::ostream &o, const Block<U> &b)
......
......@@ -105,9 +105,9 @@ public:
#ifdef FABULOUS_DEBUG_MODE
Hm = _original.sub_block(0, 0, m, n);
#else
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED. "
" You must activate FABULOUS_DEBUG_MODE to use this function"
FABULOUS_THROW(
Unsupported,
"You must define the FABULOUS_DEBUG_MODE preprocessor constant to use this function"
);
(void) Hm;
(void) m;
......@@ -182,8 +182,9 @@ private:
MORSE_desc_t *T = _tau[i].get();
int err = chameleon::ormqr<S>(trans, A, T, C, _seq.get());
if (err != 0)
FABULOUS_FATAL_ERROR("return value of ormqr step "<<i<<" is "<<err);
if (err != 0) {
FABULOUS_THROW(Kernel, "ormqr 'step' err="<<err);
}
}
/* STEP 2: Call QR on the very last block of last col */
......@@ -191,8 +192,9 @@ private:
MORSE_desc_t *A = get_sub_hess(k, k, 2, 1);
int err = chameleon::geqrf<S>(A, tau, _seq.get());
if (err != 0)
FABULOUS_FATAL_ERROR("return values of geqrf (last block) is"<<err);
if (err != 0) {
FABULOUS_THROW(Kernel, "geqrf 'last block' err="<<err);
}
_tau.push_back(tau);
FABULOUS_DEBUG("current_block="<<k);
......@@ -200,8 +202,9 @@ private:
/* STEP 3: Apply Q^H generated at step 2 to last block of RHS */
MORSE_desc_t *C = get_sub_rhs(k, 0, 2, 1);
err = chameleon::ormqr<S>(trans, A, tau.get(), C, _seq.get());
if (err != 0)
FABULOUS_FATAL_ERROR("return value of ormqrf (RHS) is "<<err);
if (err != 0) {
FABULOUS_THROW(Kernel, "ormqr 'RHS' err="<<err);
}
_last_column_factorized = true;
_logger.notify_facto_end();
......@@ -258,8 +261,7 @@ public:
Block<S> get_H1new()
{
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED: HessChamQR not compatible with DeflatedRestarting");
FABULOUS_THROW(Unsupported, "HessChamQR not compatible with DeflatedRestarting");
return Block<S>{};
}
......@@ -301,8 +303,7 @@ public:
void compute_proj_residual(const Block<S>&, const Block<S>&)
{
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED: HessChamQR not compatible with DeflatedRestarting");
FABULOUS_THROW(Unsupported, "HessChamQR not compatible with DeflatedRestarting");
}
std::tuple<P,P,bool>
......@@ -341,8 +342,7 @@ public:
Block<S> get_hess_block()
{
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED: HessChamQR not compatible with DeflatedRestarting");
FABULOUS_THROW(Unsupported, "HessChamQR not compatible with DeflatedRestarting");
return Block<S>{};
}
......@@ -359,8 +359,7 @@ public:
void reserve_DR(int)
{
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED: HessChamQR not compatible with DeflatedRestarting");
FABULOUS_THROW(Unsupported, "HessChamQR not compatible with DeflatedRestarting");
}
inline int get_nb_vect() const { return (_nb_block_col) * _nbRHS; }
......
......@@ -131,9 +131,10 @@ public:
H.get_ptr(), H.get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim()
);
if (err != 0) {
FABULOUS_THROW(Kernel, "gels (least square) err="<<err);
}
if (err != 0)
FABULOUS_FATAL_ERROR("Error in GELS err="<<err);
_solution_computed = true;
_logger.notify_least_square_end();
}
......
......@@ -128,9 +128,10 @@ public:
H.get_ptr(), H.get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim()
);
if (err != 0) {
FABULOUS_THROW(Kernel, "gels (least square) err="<<err);
}
if (err != 0)
FABULOUS_FATAL_ERROR("Error in GELS err="<<err);
_solution_computed = true;
_logger.notify_least_square_end();
}
......
......@@ -82,11 +82,12 @@ public:
void init_lambda(const Block<S> &lambda)
{
//Check
if (lambda.get_nb_col() != _nbRHS)
std::cout<<"Lambda nb_col is wrong while initiating lambda in Hess\n";
if (lambda.get_nb_row() != _nbRHS)
std::cout<<"Lambda nb_row is wrong while initiating lambda in Hess\n";
if (lambda.get_nb_col() != _nbRHS || lambda.get_nb_row() != _nbRHS) {
FABULOUS_THROW(
Internal,
"Lambda dimensions are wrong while initiating lambda in Hess\n"
);
}
_lambda1 = lambda;
}
......@@ -96,8 +97,9 @@ public:
*/
void init_phi(int p1)
{
if (_phi)
FABULOUS_FATAL_ERROR("THIS FUNCTION CANNOT BE CALLED TWICE!");
if (_phi) {
FABULOUS_THROW(Internal, "This function cannot be called twice!");
}
int N = _lambda1.get_nb_col();
_phi.reset(new Block<S>{N + p1, N});
......@@ -269,8 +271,9 @@ public:
F.get_ptr(), F.get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim()
);
if (err != 0)
FABULOUS_FATAL_ERROR("Problem on Least Square err="<<err);
if (err != 0) {
FABULOUS_THROW(Kernel, "gels (least square) err="<<err);
}
_solution_computed = true;
_logger.notify_least_square_end();
}
......
......@@ -91,11 +91,13 @@ public:
void init_lambda(const Block<S> &lambda)
{
// Check
if (lambda.get_nb_col() != _nbRHS)
FABULOUS_FATAL_ERROR("Lambda nb_col is wrong while initiating lambda in Hess\n");
if (lambda.get_nb_row() != _nbRHS + _nb_eigen_pair)
FABULOUS_FATAL_ERROR("Lambda nb_row is wrong while initiating lambda in Hess\n");
if (lambda.get_nb_col() != _nbRHS
|| lambda.get_nb_row() != _nbRHS + _nb_eigen_pair) {
FABULOUS_THROW(
Internal,
"Lambda dimensions are wrong while initiating lambda in Hess\n"
);
}
_lambda1 = lambda;
}
......@@ -105,8 +107,9 @@ public:
*/
void init_phi(int p1)
{
if (_phi)
FABULOUS_FATAL_ERROR("THIS FUNCTION CANNOT BE CALLED TWICE!");
if (_phi) {
FABULOUS_THROW(Internal, "This function cannot be called twice!");
}
int N = _lambda1.get_nb_col();
_phi.reset(new Block<S>{N + p1, N});
......@@ -121,12 +124,13 @@ public:
*/
void init_phi_restart(int p1)
{
assert(p1 == _nb_eigen_pair);
assert( _lambda1.get_nb_row() == p1+_nbRHS );
FABULOUS_ASSERT( p1 == _nb_eigen_pair );
FABULOUS_ASSERT( _lambda1.get_nb_row() == p1+_nbRHS );
_restarted = true;
if (_phi)
FABULOUS_FATAL_ERROR("THIS FUNCTION CANNOT BE CALLED TWICE!");
if (_phi) {
FABULOUS_THROW(Internal, "This function cannot be called twice!");
}
int N = _lambda1.get_nb_col();
_phi.reset(new Block<S>{N + p1, N});
......@@ -341,8 +345,9 @@ public:
F.get_ptr(), F.get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim()
);
if (err != 0)
FABULOUS_FATAL_ERROR("Problem on Least Square err="<<err);
if (err != 0) {
FABULOUS_THROW(Kernel, "gels (least square) err="<<err);
}
_solution_computed = true;
_logger.notify_least_square_end();
}
......
<
......@@ -79,8 +79,9 @@ private:
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);
if (err != 0) {
FABULOUS_THROW(Kernel, "ormqr 'step' err="<<err);
}
}
/* STEP 2: Call QR on the very last block of last col */
_tau.emplace_back();
......@@ -91,8 +92,9 @@ private: