Commit 7898fe9f authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Improve code organization in main loop (beginning)

parent 57a97bdc
......@@ -124,6 +124,15 @@ private:
X_plus_MVY(A, X, _base, Y); // X_n = X_0 + M*V*Y
}
template<class Matrix>
void do_ortho(Matrix &A, OrthoParam &orthoparm, Block<S> &W)
{
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A); // Ortho and filling of L (the (IB) Hessenberg)
_logger.notify_ortho_end();
}
template<class Matrix, class Block>
bool check_real_residual_if(bool projected_residual_converged,
Matrix &A, const Block &B, Block &X,
......@@ -140,7 +149,6 @@ private:
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);
......@@ -156,6 +164,32 @@ private:
return false;
}
template<class Matrix, class Epsilon, class Restarter>
bool initial_iteration(Matrix &A, Block<S> &X, const Block<S> &B,
const Epsilon &epsilon, Restarter &restarter)
{
std::tuple<P,P,bool> MinMaxConv;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
MinMaxConv = restarter.run(_hess, epsilon);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
Block<S> R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
MinMaxConv = R1.check_precision(epsilon);
}
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
return std::get<2>(MinMaxConv);
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
......@@ -210,28 +244,11 @@ public:
{
print_header(max_krylov_space_size);
std::tuple<P,P,bool> MinMaxConv;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
MinMaxConv = restarter.run(_hess, epsilon);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
Block<S> R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
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, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
bool convergence = initial_iteration(A, X, B, epsilon, restarter);
while (!convergence && _size_krylov_space < max_krylov_space_size) {
while (!convergence && // Main Loop
_size_krylov_space + _nbRHS < max_krylov_space_size)
{
++ _iteration_count;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
Block<S> Vj = _base.get_Vj(); // V_j
......@@ -239,22 +256,16 @@ public:
Block<S> W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A);
_logger.notify_ortho_end();
do_ortho(A, orthoparm, W);
_base.increase(_nbRHS); // add W to base
_size_krylov_space = _base.get_nb_vect();
//_hess.display_hess_bitmap("Hess Before Least Square");
//_hess.display_hess_bitmap("Hess Before Least Square");
auto MinMaxConv = check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, B, X, epsilon);
//check_arnoldi_formula(A, _base, _hess, _nbRHS);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, _hess);
......@@ -269,10 +280,9 @@ public:
FABULOUS_WARNING("cold restart");
_base.reset();
} else if (restarter && !convergence) {
int M = _hess.get_nb_vect() + _nbRHS;
Block<S> LS{M, _nbRHS};
_hess.compute_proj_residual(LS);
restarter.save_LS(LS);
Block<S> PR{_hess.get_nb_vect() + _nbRHS, _nbRHS};
_hess.compute_proj_residual(PR);
restarter.save_LS(PR);
Block<S> H = _hess.get_hess_block();
restarter.save_hess(H);
}
......
......@@ -66,6 +66,7 @@ private:
int _iteration_count;
// IB SPECIFIC:
Block<S> _WP;
int _nb_direction_kept;
int _old_nb_direction_kept;
int _nb_direction_discarded;
......@@ -363,6 +364,7 @@ public:
_size_krylov_space{0},
_iteration_count{0},
_nb_direction_kept{nbRHS},
_WP{dim, _nbRHS},
_old_nb_direction_kept{nbRHS},
_nb_direction_discarded{0},
_IB_happened{false},
......
......@@ -303,10 +303,10 @@ private:
* \param max_krylov_space_size maximum space (nb_vect) of krylov space base
* \return whether max_krylov_space_size was not reached when adding new vector into base
*/
bool do_IB_update(int max_krylov_space_size)
void do_IB_update()
{
if (!_ib_update_scheduled) {
return false;
return;
}
Block<S> &U = _U;
......@@ -331,8 +331,6 @@ private:
_base.increase(_nb_direction_kept);
_size_krylov_space = _base.get_nb_vect();
if (_size_krylov_space > max_krylov_space_size)
return true;
// PART 2
_WP.update_P( W2 ); // Pj = [P_{j-1},~Wj] * \W_2
......@@ -342,7 +340,6 @@ private:
// Gj = | \W_2^{H} | * H_j
_F.update_bottom_line(W1_W2);
_ib_update_scheduled = false;
return false;
}
template<class P>
......@@ -442,23 +439,19 @@ public:
assert( _dim == B.get_nb_row() );
// IS_INVARIANT( _nb_direction_discarded + _nb_direction_kept == nbRHS );
const auto inv_epsilon = array_inverse(epsilon);
initial_iteration(A, X, B, restarter, epsilon, inv_epsilon);
// Should not be true (because we already tested it):
convergence = (_nb_direction_discarded == _nbRHS);
while (!convergence && // Main Loop
_size_krylov_space < max_krylov_space_size)
_size_krylov_space+_nb_direction_kept < max_krylov_space_size)
{
++ _iteration_count;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (do_IB_update(max_krylov_space_size)) {
break;
}
do_IB_update();
Block<S> Vj = _base.get_Vj();
W_AMV(A, Vj, _WP); // W <- A*M*Vj: perform matrix multiplication
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment