From: Martin Kronbichler Date: Mon, 12 Sep 2022 06:40:20 +0000 (+0200) Subject: SolverGMRES: Add 'batched' mode for reduced overhead at small sizes X-Git-Tag: v9.5.0-rc1~980^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=56ce56b1c9dec5d6f40cfa64df04676a70cea30f;p=dealii.git SolverGMRES: Add 'batched' mode for reduced overhead at small sizes --- diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 3b4e3f19ed..79b7bd82e9 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -187,12 +187,14 @@ public: * Constructor. By default, set the number of temporary vectors to 30, * i.e. do a restart every 28 iterations. Also set preconditioning from * left, the residual of the stopping criterion to the default residual, - * and re-orthogonalization only if necessary. + * and re-orthogonalization only if necessary. Also, the batched mode with + * reduced functionality to track information is disabled by default. */ explicit AdditionalData(const unsigned int max_n_tmp_vectors = 30, const bool right_preconditioning = false, const bool use_default_residual = true, - const bool force_re_orthogonalization = false); + const bool force_re_orthogonalization = false, + const bool batched_mode = false); /** * Maximum number of temporary vectors. This parameter controls the size @@ -223,6 +225,14 @@ public: * if necessary. */ bool force_re_orthogonalization; + + /** + * Flag to control whether a reduced mode of the solver should be + * run. This is necessary when running (several) SolverGMRES instances + * involving very small and cheap linear systems where the feedback from + * all signals, eigenvalue computations, and log stream are disabled. + */ + bool batched_mode; }; /** @@ -372,6 +382,13 @@ protected: */ boost::signals2::signal re_orthogonalize_signal; + /** + * A reference to the underlying SolverControl object. In the regular case, + * this is not needed, as the signal from the base class is used, but the + * batched variant cannot use those mechanisms due to the high costs. + */ + SolverControl &solver_control; + /** * Implementation of the computation of the norm of the residual. */ @@ -433,11 +450,28 @@ protected: FullMatrix H; /** - * Auxiliary matrix for inverting @p H + * Auxiliary vector for orthogonalization */ - FullMatrix H1; + Vector gamma; + + /** + * Auxiliary vector for orthogonalization + */ + Vector ci; + + /** + * Auxiliary vector for orthogonalization + */ + Vector si; + + /** + * Auxiliary vector for orthogonalization + */ + Vector h; }; + + /** * Implementation of the Generalized minimal residual method with flexible * preconditioning (flexible GMRES or FGMRES). @@ -560,7 +594,7 @@ namespace internal if (data[i] == nullptr) { data[i] = std::move(typename VectorMemory::Pointer(mem)); - data[i]->reinit(temp); + data[i]->reinit(temp, true); } return *data[i]; } @@ -584,6 +618,24 @@ namespace internal return x.real() < y.real() || (x.real() == y.real() && x.imag() < y.imag()); } + + // A function to solve the (upper) triangular system after Givens + // rotations on a matrix that has possibly unused rows and columns + void + solve_triangular(const unsigned int dim, + const FullMatrix &H, + const Vector & rhs, + Vector & solution) + { + for (int i = dim - 1; i >= 0; --i) + { + double s = rhs(i); + for (unsigned int j = i + 1; j < dim; ++j) + s -= solution(j) * H(i, j); + solution(i) = s / H(i, i); + AssertIsFinite(solution(i)); + } + } } // namespace SolverGMRESImplementation } // namespace internal @@ -594,11 +646,13 @@ inline SolverGMRES::AdditionalData::AdditionalData( const unsigned int max_n_tmp_vectors, const bool right_preconditioning, const bool use_default_residual, - const bool force_re_orthogonalization) + const bool force_re_orthogonalization, + const bool batched_mode) : max_n_tmp_vectors(max_n_tmp_vectors) , right_preconditioning(right_preconditioning) , use_default_residual(use_default_residual) , force_re_orthogonalization(force_re_orthogonalization) + , batched_mode(batched_mode) { Assert(3 <= max_n_tmp_vectors, ExcMessage("SolverGMRES needs at least three " @@ -613,6 +667,7 @@ SolverGMRES::SolverGMRES(SolverControl & cn, const AdditionalData & data) : SolverBase(cn, mem) , additional_data(data) + , solver_control(cn) {} @@ -622,6 +677,7 @@ SolverGMRES::SolverGMRES(SolverControl & cn, const AdditionalData &data) : SolverBase(cn) , additional_data(data) + , solver_control(cn) {} @@ -791,7 +847,9 @@ SolverGMRES::solve(const MatrixType & A, // TODO:[GK] Make sure the parameter in the constructor means maximum basis // size - LogStream::Prefix prefix("GMRES"); + std::unique_ptr prefix; + if (!additional_data.batched_mode) + prefix = std::make_unique("GMRES"); // extra call to std::max to placate static analyzers: coverity rightfully // complains that data.max_n_tmp_vectors - 2 may overflow @@ -808,9 +866,11 @@ SolverGMRES::solve(const MatrixType & A, unsigned int accumulated_iterations = 0; const bool do_eigenvalues = - !condition_number_signal.empty() || !all_condition_numbers_signal.empty() || - !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() || - !hessenberg_signal.empty() || !all_hessenberg_signal.empty(); + !additional_data.batched_mode && + (!condition_number_signal.empty() || + !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() || + !all_eigenvalues_signal.empty() || !hessenberg_signal.empty() || + !all_hessenberg_signal.empty()); // for eigenvalue computation, need to collect the Hessenberg matrix (before // applying Givens rotations) FullMatrix H_orig; @@ -821,9 +881,10 @@ SolverGMRES::solve(const MatrixType & A, H.reinit(n_tmp_vectors, n_tmp_vectors - 1); // some additional vectors, also used in the orthogonalization - dealii::Vector gamma(n_tmp_vectors), ci(n_tmp_vectors - 1), - si(n_tmp_vectors - 1), h(n_tmp_vectors - 1); - + gamma.reinit(n_tmp_vectors); + ci.reinit(n_tmp_vectors - 1); + si.reinit(n_tmp_vectors - 1); + h.reinit(n_tmp_vectors - 1); unsigned int dim = 0; @@ -889,8 +950,11 @@ SolverGMRES::solve(const MatrixType & A, if (use_default_residual) { last_res = rho; - iteration_state = - this->iteration_status(accumulated_iterations, rho, x); + if (additional_data.batched_mode) + iteration_state = solver_control.check(accumulated_iterations, rho); + else + iteration_state = + this->iteration_status(accumulated_iterations, rho, x); if (iteration_state != SolverControl::iterate) break; @@ -909,8 +973,11 @@ SolverGMRES::solve(const MatrixType & A, double res = r->l2_norm(); last_res = res; - iteration_state = - this->iteration_status(accumulated_iterations, res, x); + if (additional_data.batched_mode) + iteration_state = solver_control.check(accumulated_iterations, rho); + else + iteration_state = + this->iteration_status(accumulated_iterations, res, x); if (iteration_state != SolverControl::iterate) break; @@ -978,32 +1045,33 @@ SolverGMRES::solve(const MatrixType & A, if (use_default_residual) { last_res = rho; - iteration_state = - this->iteration_status(accumulated_iterations, rho, x); + if (additional_data.batched_mode) + iteration_state = + solver_control.check(accumulated_iterations, rho); + else + iteration_state = + this->iteration_status(accumulated_iterations, rho, x); } else { - deallog << "default_res=" << rho << std::endl; + if (!additional_data.batched_mode) + deallog << "default_res=" << rho << std::endl; - dealii::Vector h_(dim); *x_ = x; *gamma_ = gamma; - H1.reinit(dim + 1, dim); - - for (unsigned int i = 0; i < dim + 1; ++i) - for (unsigned int j = 0; j < dim; ++j) - H1(i, j) = H(i, j); - - H1.backward(h_, *gamma_); + internal::SolverGMRESImplementation::solve_triangular(dim, + H, + *gamma_, + h); if (left_precondition) for (unsigned int i = 0; i < dim; ++i) - x_->add(h_(i), tmp_vectors[i]); + x_->add(h(i), tmp_vectors[i]); else { p = 0.; for (unsigned int i = 0; i < dim; ++i) - p.add(h_(i), tmp_vectors[i]); + p.add(h(i), tmp_vectors[i]); preconditioner.vmult(*r, p); x_->add(1., *r); }; @@ -1024,37 +1092,36 @@ SolverGMRES::solve(const MatrixType & A, const double preconditioned_res = x_->l2_norm(); last_res = preconditioned_res; - iteration_state = - this->iteration_status(accumulated_iterations, - preconditioned_res, - x); + if (additional_data.batched_mode) + iteration_state = + solver_control.check(accumulated_iterations, rho); + else + iteration_state = + this->iteration_status(accumulated_iterations, + preconditioned_res, + x); } } - }; + } + // end of inner iteration. now calculate the solution from the temporary // vectors - h.reinit(dim); - H1.reinit(dim + 1, dim); + internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h); - for (unsigned int i = 0; i < dim + 1; ++i) - for (unsigned int j = 0; j < dim; ++j) - H1(i, j) = H(i, j); - - compute_eigs_and_cond(H_orig, - dim, - all_eigenvalues_signal, - all_hessenberg_signal, - condition_number_signal); - - H1.backward(h, gamma); + if (do_eigenvalues) + compute_eigs_and_cond(H_orig, + dim, + all_eigenvalues_signal, + all_hessenberg_signal, + condition_number_signal); if (left_precondition) for (unsigned int i = 0; i < dim; ++i) x.add(h(i), tmp_vectors[i]); else { - p = 0.; - for (unsigned int i = 0; i < dim; ++i) + p.equ(h(0), tmp_vectors[0]); + for (unsigned int i = 1; i < dim; ++i) p.add(h(i), tmp_vectors[i]); preconditioner.vmult(v, p); x.add(1., v); @@ -1064,13 +1131,14 @@ SolverGMRES::solve(const MatrixType & A, } while (iteration_state == SolverControl::iterate); - compute_eigs_and_cond(H_orig, - dim, - eigenvalues_signal, - hessenberg_signal, - condition_number_signal); + if (do_eigenvalues) + compute_eigs_and_cond(H_orig, + dim, + eigenvalues_signal, + hessenberg_signal, + condition_number_signal); - if (!krylov_space_signal.empty()) + if (!additional_data.batched_mode && !krylov_space_signal.empty()) krylov_space_signal(tmp_vectors); // in case of failure: throw exception