From: Martin Kronbichler Date: Mon, 11 Mar 2024 14:50:56 +0000 (+0100) Subject: Cleanup of SolverGMRES and SolverFGMRES implementations X-Git-Tag: v9.6.0-rc1~482^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b70ca60cd41aafa079d87f96ccd333dc6eacd2c6;p=dealii.git Cleanup of SolverGMRES and SolverFGMRES implementations --- diff --git a/examples/step-12/step-12.cc b/examples/step-12/step-12.cc index 511acf5dab..6c7de2cbde 100644 --- a/examples/step-12/step-12.cc +++ b/examples/step-12/step-12.cc @@ -481,7 +481,7 @@ namespace Step12 SolverControl solver_control(1000, 1e-6 * right_hand_side.l2_norm()); SolverGMRES>::AdditionalData additional_data; - additional_data.max_n_tmp_vectors = 100; + additional_data.max_basis_size = 100; SolverGMRES> solver(solver_control, additional_data); // Here we create the preconditioner, diff --git a/examples/step-22/doc/results.dox b/examples/step-22/doc/results.dox index 32f7be0f72..30473dc225 100644 --- a/examples/step-22/doc/results.dox +++ b/examples/step-22/doc/results.dox @@ -544,7 +544,7 @@ $k=100$ temporary vectors: 1e-6*system_rhs.l2_norm()); GrowingVectorMemory > vector_memory; SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; + gmres_data.max_basis_size = 100; SolverGMRES > gmres(solver_control, vector_memory, gmres_data); diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index f722f144b4..78cf154529 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -25,7 +25,6 @@ #include #include -#include #include #include #include @@ -129,11 +128,14 @@ namespace internal * Implementation of the Restarted Preconditioned Direct Generalized Minimal * Residual Method. The stopping criterion is the norm of the residual. * - * The AdditionalData structure contains the number of temporary vectors used. - * The size of the Arnoldi basis is this number minus three. Additionally, it - * allows you to choose between right or left preconditioning. The default is - * left preconditioning. Finally it includes a flag indicating whether or not - * the default residual is used as stopping criterion. + * The AdditionalData structure contains the size of the Arnoldi basis used + * for orthogonalization. It is related to the number of temporary vectors + * used, which is the basis size plus two. Additionally, it allows you to + * choose between right or left preconditioning. The default is left + * preconditioning. Furthermore, it includes a flag indicating whether or not + * the default residual is used as stopping criterion and an option for the + * orthogonalization algorithm, see LinearAlgebra::OrthogonalizationStrategy + * for available options. * * *

Left versus right preconditioning

@@ -142,7 +144,7 @@ namespace internal * preconditioning. As expected, this switches between solving for the systems * P-1A and AP-1, respectively. * - * A second consequence is the type of residual which is used to measure + * A second consequence is the type of residual, which is used to measure * convergence. With left preconditioning, this is the preconditioned * residual, while with right preconditioning, it is the residual of the * unpreconditioned system. @@ -156,13 +158,13 @@ namespace internal * *

The size of the Arnoldi basis

* - * The maximal basis size is controlled by AdditionalData::max_n_tmp_vectors, - * and it is this number minus 2. If the number of iteration steps exceeds - * this number, all basis vectors are discarded and the iteration starts anew - * from the approximation obtained so far. + * The maximal basis size is controlled by AdditionalData::max_basis_size. If + * the number of iteration steps exceeds this number, all basis vectors are + * discarded and the iteration starts anew from the approximation obtained so + * far. * - * Note that the minimizing property of GMRes only pertains to the Krylov - * space spanned by the Arnoldi basis. Therefore, restarted GMRes is + * Note that the minimizing property of GMRES only pertains to the Krylov + * space spanned by the Arnoldi basis. Therefore, restarted GMRES is * not minimizing anymore. The choice of the basis length is a * trade-off between memory consumption and convergence speed, since a longer * basis means minimization over a larger space. @@ -199,14 +201,15 @@ public: struct AdditionalData { /** - * 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. Also, the batched mode with - * reduced functionality to track information is disabled by default. + * Constructor. By default, set the size of the Arnoldi basis. Also set + * preconditioning from left, the residual of the stopping criterion to + * the default residual, and re-orthogonalization only if necessary. Also, + * the batched mode with reduced functionality to track information is + * disabled by default. Finally, the default orthogonalization algorithm + * is the modified Gram-Schmidt method. */ explicit AdditionalData( - const unsigned int max_n_tmp_vectors = 30, + const unsigned int max_basis_size = 30, const bool right_preconditioning = false, const bool use_default_residual = true, const bool force_re_orthogonalization = false, @@ -216,13 +219,27 @@ public: LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt); /** - * Maximum number of temporary vectors. This parameter controls the size - * of the Arnoldi basis, which for historical reasons is - * #max_n_tmp_vectors-2. SolverGMRES assumes that there are at least three - * temporary vectors, so this value must be greater than or equal to three. + * Maximum number of temporary vectors. Together with #max_basis_size, + * this parameter controls the size of the Arnoldi basis, which for + * historical reasons is #max_n_tmp_vectors-2. SolverGMRES assumes that + * there are at least three temporary vectors, so this value must be + * greater than or equal to three. If both this variable and + * #max_basis_size are set to a non-zero value, the solver uses the latter + * variable. + * + * @deprecated Use #max_basis_size instead. */ unsigned int max_n_tmp_vectors; + /** + * Maximum size of the Arnoldi basis. SolverGMRES assumes that there is at + * least one vector in the Arnoldi basis, so this value must be greater + * than or equal to one. Note that whenever this variable is set to a + * non-zero, including the value set by the default constructor, this + * variable takes precedence. + */ + unsigned int max_basis_size; + /** * Flag for right preconditioning. * @@ -419,17 +436,6 @@ protected: virtual double criterion(); - /** - * Transformation of an upper Hessenberg matrix into tridiagonal structure - * by givens rotation of the last column - */ - void - givens_rotation(Vector &h, - Vector &b, - Vector &ci, - Vector &si, - int col) const; - /** * Estimates the eigenvalues from the Hessenberg matrix, H_orig, generated * during the inner iterations. Uses these estimate to compute the condition @@ -454,17 +460,12 @@ protected: /** * Auxiliary vector for orthogonalization */ - Vector gamma; - - /** - * Auxiliary vector for orthogonalization - */ - Vector ci; + Vector projected_rhs; /** * Auxiliary vector for orthogonalization */ - Vector si; + std::vector> givens_rotations; /** * Auxiliary vector for orthogonalization @@ -572,6 +573,51 @@ private: #ifndef DOXYGEN + +template +inline SolverGMRES::AdditionalData::AdditionalData( + const unsigned int max_basis_size, + const bool right_preconditioning, + const bool use_default_residual, + const bool force_re_orthogonalization, + const bool batched_mode, + const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy) + : max_n_tmp_vectors(0) + , max_basis_size(max_basis_size) + , right_preconditioning(right_preconditioning) + , use_default_residual(use_default_residual) + , force_re_orthogonalization(force_re_orthogonalization) + , batched_mode(batched_mode) + , orthogonalization_strategy(orthogonalization_strategy) +{ + Assert(max_basis_size >= 1, + ExcMessage("SolverGMRES needs at least one vector in the " + "Arnoldi basis.")); +} + + + +template +SolverGMRES::SolverGMRES(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data) + : SolverBase(cn, mem) + , additional_data(data) + , solver_control(cn) +{} + + + +template +SolverGMRES::SolverGMRES(SolverControl &cn, + const AdditionalData &data) + : SolverBase(cn) + , additional_data(data) + , solver_control(cn) +{} + + + namespace internal { namespace SolverGMRESImplementation @@ -622,111 +668,6 @@ namespace internal - // A comparator for better printing eigenvalues - inline bool - complex_less_pred(const std::complex &x, - const std::complex &y) - { - 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 - inline 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 - - - -template -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 batched_mode, - const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy) - : 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) - , orthogonalization_strategy(orthogonalization_strategy) -{ - Assert(3 <= max_n_tmp_vectors, - ExcMessage("SolverGMRES needs at least three " - "temporary vectors.")); -} - - - -template -SolverGMRES::SolverGMRES(SolverControl &cn, - VectorMemory &mem, - const AdditionalData &data) - : SolverBase(cn, mem) - , additional_data(data) - , solver_control(cn) -{} - - - -template -SolverGMRES::SolverGMRES(SolverControl &cn, - const AdditionalData &data) - : SolverBase(cn) - , additional_data(data) - , solver_control(cn) -{} - - - -template -inline void -SolverGMRES::givens_rotation(Vector &h, - Vector &b, - Vector &ci, - Vector &si, - int col) const -{ - for (int i = 0; i < col; ++i) - { - const double s = si(i); - const double c = ci(i); - const double dummy = h(i); - h(i) = c * dummy + s * h(i + 1); - h(i + 1) = -s * dummy + c * h(i + 1); - }; - - const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1)); - si(col) = h(col + 1) * r; - ci(col) = h(col) * r; - h(col) = ci(col) * h(col) + si(col) * h(col + 1); - b(col + 1) = -si(col) * b(col); - b(col) *= ci(col); -} - - - -namespace internal -{ - namespace SolverGMRESImplementation - { template struct is_dealii_compatible_distributed_vector; @@ -1225,6 +1166,65 @@ namespace internal return 0.0; } + + // A comparator for better printing eigenvalues + inline bool + complex_less_pred(const std::complex &x, + const std::complex &y) + { + return x.real() < y.real() || + (x.real() == y.real() && x.imag() < y.imag()); + } + + + + // A function to compute the Givens rotation for the QR factorization of + // the Hessenberg matrix involved in the Arnoldi process, transforming it + // into an upper triangular matrix. + inline void + givens_rotation(Vector &h, + Vector &b, + std::vector> &rotations, + const int col) + { + for (int i = 0; i < col; ++i) + { + const double c = rotations[i].first; + const double s = rotations[i].second; + const double tmp = h(i); + h(i) = c * tmp + s * h(i + 1); + h(i + 1) = -s * tmp + c * h(i + 1); + } + + const double r = + 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1)); + rotations[col].second = h(col + 1) * r; + rotations[col].first = h(col) * r; + h(col) = + rotations[col].first * h(col) + rotations[col].second * h(col + 1); + b(col + 1) = -rotations[col].second * b(col); + b(col) *= rotations[col].first; + } + + + + // A function to solve the (upper) triangular system after Givens + // rotations on a matrix that has possibly unused rows and columns + inline 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 @@ -1290,22 +1290,20 @@ SolverGMRES::solve(const MatrixType &A, const VectorType &b, const PreconditionerType &preconditioner) { - // TODO:[?] Check, why there are two different start residuals. - // TODO:[GK] Make sure the parameter in the constructor means maximum basis - // size - 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 - const unsigned int n_tmp_vectors = - std::max(additional_data.max_n_tmp_vectors, 3u); + const unsigned int basis_size = + (additional_data.max_basis_size > 0 ? + additional_data.max_basis_size : + std::max(additional_data.max_n_tmp_vectors, 3u) - 2); // Generate an object where basis vectors are stored. - internal::SolverGMRESImplementation::TmpVectors tmp_vectors( - n_tmp_vectors, this->memory); + internal::SolverGMRESImplementation::TmpVectors basis_vectors( + basis_size + 2, this->memory); // number of the present iteration; this // number is not reset to zero upon a @@ -1322,21 +1320,18 @@ SolverGMRES::solve(const MatrixType &A, // applying Givens rotations) FullMatrix H_orig; if (do_eigenvalues) - H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1); + H_orig.reinit(basis_size + 1, basis_size); // matrix used for the orthogonalization process later - H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true); + H.reinit(basis_size + 1, basis_size, /* omit_initialization */ true); // some additional vectors, also used in the orthogonalization - 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; + projected_rhs.reinit(basis_size + 1); + givens_rotations.resize(basis_size); + h.reinit(basis_size + 1); SolverControl::State iteration_state = SolverControl::iterate; - double last_res = std::numeric_limits::lowest(); + double res = std::numeric_limits::lowest(); // switch to determine whether we want a left or a right preconditioner. at // present, left is default, but both ways are implemented @@ -1347,15 +1342,14 @@ SolverGMRES::solve(const MatrixType &A, // residual as stopping criterion. const bool use_default_residual = additional_data.use_default_residual; - // define two aliases - VectorType &v = tmp_vectors(0, x); - VectorType &p = tmp_vectors(n_tmp_vectors - 1, x); + // define an alias + VectorType &p = basis_vectors(basis_size + 1, x); // Following vectors are needed when we are not using the default residuals // as stopping criterion typename VectorMemory::Pointer r; typename VectorMemory::Pointer x_; - std::unique_ptr> gamma_; + std::unique_ptr> gamma; if (!use_default_residual) { r = std::move(typename VectorMemory::Pointer(this->memory)); @@ -1363,7 +1357,7 @@ SolverGMRES::solve(const MatrixType &A, r->reinit(x); x_->reinit(x); - gamma_ = std::make_unique>(gamma.size()); + gamma = std::make_unique>(projected_rhs.size()); } bool re_orthogonalize = additional_data.force_re_orthogonalization; @@ -1374,42 +1368,45 @@ SolverGMRES::solve(const MatrixType &A, // restart do { - double rho = 0.0; + VectorType &v = basis_vectors(0, x); + double norm_v = 0.; if (left_precondition) { A.vmult(p, x); p.sadd(-1., 1., b); preconditioner.vmult(v, p); - rho = v.l2_norm(); + norm_v = v.l2_norm(); } else { A.vmult(v, x); - rho = dealii::internal::SolverGMRESImplementation::sadd_and_norm(v, - -1, - b, - 1.0); + norm_v = dealii::internal::SolverGMRESImplementation::sadd_and_norm( + v, -1, b, 1.0); } + projected_rhs(0) = norm_v; + if (norm_v != 0) + v /= norm_v; + // check the residual here as well since it may be that we got the exact // (or an almost exact) solution vector at the outset. if we wouldn't // check here, the next scaling operation would produce garbage if (use_default_residual) { - last_res = rho; + res = norm_v; if (additional_data.batched_mode) - iteration_state = solver_control.check(accumulated_iterations, rho); + iteration_state = solver_control.check(accumulated_iterations, res); else iteration_state = - this->iteration_status(accumulated_iterations, rho, x); + this->iteration_status(accumulated_iterations, res, x); if (iteration_state != SolverControl::iterate) break; } else { - deallog << "default_res=" << rho << std::endl; + deallog << "default_res=" << norm_v << std::endl; if (left_precondition) { @@ -1419,10 +1416,9 @@ SolverGMRES::solve(const MatrixType &A, else preconditioner.vmult(*r, v); - double res = r->l2_norm(); - last_res = res; + res = r->l2_norm(); if (additional_data.batched_mode) - iteration_state = solver_control.check(accumulated_iterations, rho); + iteration_state = solver_control.check(accumulated_iterations, res); else iteration_state = this->iteration_status(accumulated_iterations, res, x); @@ -1431,98 +1427,91 @@ SolverGMRES::solve(const MatrixType &A, break; } - gamma(0) = rho; - - v *= 1. / rho; - - // inner iteration doing at most as many steps as there are temporary - // vectors. the number of steps actually been done is propagated outside - // through the @p dim variable - for (unsigned int inner_iteration = 0; - ((inner_iteration < n_tmp_vectors - 2) && - (iteration_state == SolverControl::iterate)); + // inner iteration doing at most as many steps as the size of the + // Arnoldi basis + unsigned int inner_iteration = 0; + for (; (inner_iteration < basis_size && + iteration_state == SolverControl::iterate); ++inner_iteration) { ++accumulated_iterations; // yet another alias - VectorType &vv = tmp_vectors(inner_iteration + 1, x); + VectorType &vv = basis_vectors(inner_iteration + 1, x); if (left_precondition) { - A.vmult(p, tmp_vectors[inner_iteration]); + A.vmult(p, basis_vectors[inner_iteration]); preconditioner.vmult(vv, p); } else { - preconditioner.vmult(p, tmp_vectors[inner_iteration]); + preconditioner.vmult(p, basis_vectors[inner_iteration]); A.vmult(vv, p); } - dim = inner_iteration + 1; - - const double s = - internal::SolverGMRESImplementation::iterated_gram_schmidt( - additional_data.orthogonalization_strategy, - tmp_vectors, - dim, - accumulated_iterations, - vv, - h, - re_orthogonalize, - re_orthogonalize_signal); - h(inner_iteration + 1) = s; - - // s=0 is a lucky breakdown, the solver will reach convergence, - // but we must not divide by zero here. - if (s != 0) - vv *= 1. / s; + norm_v = internal::SolverGMRESImplementation::iterated_gram_schmidt( + additional_data.orthogonalization_strategy, + basis_vectors, + inner_iteration + 1, + accumulated_iterations, + vv, + h, + re_orthogonalize, + re_orthogonalize_signal); + + // norm_v = 0 is a lucky breakdown, the solver will reach + // convergence, but we must not divide by zero here. + if (norm_v != 0) + vv /= norm_v; + + h(inner_iteration + 1) = norm_v; // for eigenvalues, get the resulting coefficients from the // orthogonalization process if (do_eigenvalues) - for (unsigned int i = 0; i < dim + 1; ++i) + for (unsigned int i = 0; i < inner_iteration + 2; ++i) H_orig(i, inner_iteration) = h(i); // Transformation into tridiagonal structure - givens_rotation(h, gamma, ci, si, inner_iteration); + internal::SolverGMRESImplementation::givens_rotation(h, + projected_rhs, + givens_rotations, + inner_iteration); // append vector on matrix - for (unsigned int i = 0; i < dim; ++i) + for (unsigned int i = 0; i < inner_iteration + 1; ++i) H(i, inner_iteration) = h(i); // default residual - rho = std::fabs(gamma(dim)); + res = std::fabs(projected_rhs(inner_iteration + 1)); if (use_default_residual) { - last_res = rho; if (additional_data.batched_mode) iteration_state = - solver_control.check(accumulated_iterations, rho); + solver_control.check(accumulated_iterations, res); else iteration_state = - this->iteration_status(accumulated_iterations, rho, x); + this->iteration_status(accumulated_iterations, res, x); } else { if (!additional_data.batched_mode) - deallog << "default_res=" << rho << std::endl; + deallog << "default_res=" << res << std::endl; - *x_ = x; - *gamma_ = gamma; - internal::SolverGMRESImplementation::solve_triangular(dim, - H, - *gamma_, - h); + *x_ = x; + *gamma = projected_rhs; + internal::SolverGMRESImplementation::solve_triangular( + inner_iteration + 1, H, *gamma, h); if (left_precondition) - for (unsigned int i = 0; i < dim; ++i) - x_->add(h(i), tmp_vectors[i]); + for (unsigned int i = 0; i < inner_iteration + 1; ++i) + x_->add(h(i), basis_vectors[i]); else { p = 0.; - for (unsigned int i = 0; i < dim; ++i) - p.add(h(i), tmp_vectors[i]); + for (unsigned int i = 0; i < inner_iteration + 1; ++i) + p.add(h(i), basis_vectors[i]); preconditioner.vmult(*r, p); x_->add(1., *r); }; @@ -1531,69 +1520,69 @@ SolverGMRES::solve(const MatrixType &A, // Now *r contains the unpreconditioned residual!! if (left_precondition) { - const double res = r->l2_norm(); - last_res = res; - + res = r->l2_norm(); iteration_state = this->iteration_status(accumulated_iterations, res, x); } else { preconditioner.vmult(*x_, *r); - const double preconditioned_res = x_->l2_norm(); - last_res = preconditioned_res; + res = x_->l2_norm(); if (additional_data.batched_mode) iteration_state = - solver_control.check(accumulated_iterations, rho); + solver_control.check(accumulated_iterations, res); else iteration_state = - this->iteration_status(accumulated_iterations, - preconditioned_res, - x); + this->iteration_status(accumulated_iterations, res, x); } } } // end of inner iteration. now calculate the solution from the temporary // vectors - internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h); + internal::SolverGMRESImplementation::solve_triangular(inner_iteration, + H, + projected_rhs, + h); if (do_eigenvalues) compute_eigs_and_cond(H_orig, - dim, + inner_iteration, all_eigenvalues_signal, all_hessenberg_signal, condition_number_signal); if (left_precondition) dealii::internal::SolverGMRESImplementation::add( - x, dim, h, tmp_vectors, false); + x, inner_iteration, h, basis_vectors, false); else { dealii::internal::SolverGMRESImplementation::add( - p, dim, h, tmp_vectors, true); + p, inner_iteration, h, basis_vectors, true); preconditioner.vmult(v, p); x.add(1., v); } + + // in the last round, print the eigenvalues from the last Arnoldi step + if (iteration_state != SolverControl::iterate && do_eigenvalues) + compute_eigs_and_cond(H_orig, + inner_iteration, + eigenvalues_signal, + hessenberg_signal, + condition_number_signal); + // end of outer iteration. restart if no convergence and the number of // iterations is not exceeded } while (iteration_state == SolverControl::iterate); - if (do_eigenvalues) - compute_eigs_and_cond(H_orig, - dim, - eigenvalues_signal, - hessenberg_signal, - condition_number_signal); - if (!additional_data.batched_mode && !krylov_space_signal.empty()) - krylov_space_signal(tmp_vectors); + krylov_space_signal(basis_vectors); // in case of failure: throw exception AssertThrow(iteration_state == SolverControl::success, - SolverControl::NoConvergence(accumulated_iterations, last_res)); + SolverControl::NoConvergence(accumulated_iterations, res)); } @@ -1721,7 +1710,7 @@ SolverFGMRES::solve(const MatrixType &A, // Generate an object where basis vectors are stored. typename internal::SolverGMRESImplementation::TmpVectors v( - basis_size, this->memory); + basis_size + 1, this->memory); typename internal::SolverGMRESImplementation::TmpVectors z( basis_size, this->memory); @@ -1731,84 +1720,85 @@ SolverFGMRES::solve(const MatrixType &A, // matrix used for the orthogonalization process later H.reinit(basis_size + 1, basis_size); + std::vector> givens_rotations(basis_size); Vector h(basis_size + 1); // Vectors for projected system - Vector projected_rhs; - Vector y; + Vector projected_rhs(basis_size + 1); + Vector y(basis_size); // Iteration starts here double res = std::numeric_limits::lowest(); - typename VectorMemory::Pointer aux(this->memory); - aux->reinit(x); do { - A.vmult(*aux, x); - aux->sadd(-1., 1., b); + A.vmult(v(0, x), x); + v[0].sadd(-1., 1., b); - double beta = aux->l2_norm(); - res = beta; + double norm_v = v[0].l2_norm(); + res = norm_v; iteration_state = this->iteration_status(accumulated_iterations, res, x); if (iteration_state == SolverControl::success) break; H.reinit(basis_size + 1, basis_size); - double a = beta; - for (unsigned int j = 0; j < basis_size; ++j) + projected_rhs(0) = norm_v; + unsigned int inner_iteration = 0; + for (; (inner_iteration < basis_size && + iteration_state == SolverControl::iterate); + ++inner_iteration) { - if (a != 0) // treat lucky breakdown - v(j, x).equ(1. / a, *aux); - else - v(j, x) = 0.; - + // norm_v = 0 is a lucky breakdown, the solver will reach + // convergence, but we must not divide by zero here. + if (norm_v != 0) + v[inner_iteration] /= norm_v; - preconditioner.vmult(z(j, x), v[j]); - A.vmult(*aux, z[j]); + preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]); + A.vmult(v(inner_iteration + 1, x), z[inner_iteration]); // Gram-Schmidt - bool re_orthogonalize = false; - const double s = - internal::SolverGMRESImplementation::iterated_gram_schmidt< - VectorType>(additional_data.orthogonalization_strategy, - v, - j + 1, - 0, - *aux, - h, - re_orthogonalize); - for (unsigned int i = 0; i <= j; ++i) - H(i, j) = h(i); - H(j + 1, j) = a = s; + bool re_orthogonalize = false; + norm_v = internal::SolverGMRESImplementation::iterated_gram_schmidt< + VectorType>(additional_data.orthogonalization_strategy, + v, + inner_iteration + 1, + 0, + v[inner_iteration + 1], + h, + re_orthogonalize); // Compute projected solution + h(inner_iteration + 1) = norm_v; + internal::SolverGMRESImplementation::givens_rotation(h, + projected_rhs, + givens_rotations, + inner_iteration); + + // append vector on Hessenberg matrix + for (unsigned int i = 0; i < inner_iteration + 1; ++i) + H(i, inner_iteration) = h(i); - if (j > 0) - { - H1.reinit(j + 1, j); - projected_rhs.reinit(j + 1); - y.reinit(j); - projected_rhs(0) = beta; - H1.fill(H); - - // check convergence. note that the vector 'x' we pass to the - // criterion is not the final solution we compute if we - // decide to jump out of the iteration (we update 'x' again - // right after the current loop) - Householder house(H1); - res = house.least_squares(y, projected_rhs); - iteration_state = - this->iteration_status(++accumulated_iterations, res, x); - if (iteration_state != SolverControl::iterate) - break; - } + // default residual + res = std::fabs(projected_rhs(inner_iteration + 1)); + + // check convergence. note that the vector 'x' we pass to the + // criterion is not the final solution we compute if we + // decide to jump out of the iteration (we update 'x' again + // right after the current loop) + iteration_state = + this->iteration_status(++accumulated_iterations, res, x); } - // Update solution vector - for (unsigned int j = 0; j < y.size(); ++j) - x.add(y(j), z[j]); + // Solve triangular system with projected quantities and update solution + // vector + internal::SolverGMRESImplementation::solve_triangular(inner_iteration, + H, + projected_rhs, + y); + dealii::internal::SolverGMRESImplementation::add( + x, inner_iteration, y, z, false); } while (iteration_state == SolverControl::iterate); diff --git a/tests/examples/step-56.output b/tests/examples/step-56.output index 5edf0cac08..f72504944d 100644 --- a/tests/examples/step-56.output +++ b/tests/examples/step-56.output @@ -7,10 +7,10 @@ Refinement cycle 0 Assembling Multigrid... Solving... Number of FGMRES iterations: 14 - Total number of iterations used for approximation of A inverse: 45 - Total number of iterations used for approximation of S inverse: 15 + Total number of iterations used for approximation of A inverse: 42 + Total number of iterations used for approximation of S inverse: 14 - Note: The mean value was adjusted by 8.38229e-15 + Note: The mean value was adjusted by 8.58753e-15 Velocity L2 Error: 0.00340778 Pressure L2 Error: 0.017624 @@ -23,10 +23,10 @@ Refinement cycle 1 Assembling Multigrid... Solving... Number of FGMRES iterations: 16 - Total number of iterations used for approximation of A inverse: 52 - Total number of iterations used for approximation of S inverse: 17 + Total number of iterations used for approximation of A inverse: 49 + Total number of iterations used for approximation of S inverse: 16 - Note: The mean value was adjusted by -1.75663e-14 + Note: The mean value was adjusted by -1.77281e-14 Velocity L2 Error: 0.000426317 Pressure L2 Error: 0.00414007 @@ -39,10 +39,10 @@ Refinement cycle 2 Assembling Multigrid... Solving... Number of FGMRES iterations: 16 - Total number of iterations used for approximation of A inverse: 52 - Total number of iterations used for approximation of S inverse: 17 + Total number of iterations used for approximation of A inverse: 49 + Total number of iterations used for approximation of S inverse: 16 - Note: The mean value was adjusted by -2.40243e-14 + Note: The mean value was adjusted by -2.39441e-14 Velocity L2 Error: 5.33246e-05 Pressure L2 Error: 0.00102063 @@ -56,10 +56,10 @@ Refinement cycle 0 Solving... Computing preconditioner... Number of FGMRES iterations: 14 - Total number of iterations used for approximation of A inverse: 75 - Total number of iterations used for approximation of S inverse: 15 + Total number of iterations used for approximation of A inverse: 70 + Total number of iterations used for approximation of S inverse: 14 - Note: The mean value was adjusted by 3.10407e-15 + Note: The mean value was adjusted by 3.04715e-15 Velocity L2 Error: 0.00340778 Pressure L2 Error: 0.017624 @@ -72,10 +72,10 @@ Refinement cycle 1 Solving... Computing preconditioner... Number of FGMRES iterations: 16 - Total number of iterations used for approximation of A inverse: 134 - Total number of iterations used for approximation of S inverse: 17 + Total number of iterations used for approximation of A inverse: 126 + Total number of iterations used for approximation of S inverse: 16 - Note: The mean value was adjusted by -3.36177e-15 + Note: The mean value was adjusted by -3.40459e-15 Velocity L2 Error: 0.000426317 Pressure L2 Error: 0.00414007 @@ -88,10 +88,10 @@ Refinement cycle 2 Solving... Computing preconditioner... Number of FGMRES iterations: 16 - Total number of iterations used for approximation of A inverse: 229 - Total number of iterations used for approximation of S inverse: 17 + Total number of iterations used for approximation of A inverse: 216 + Total number of iterations used for approximation of S inverse: 16 - Note: The mean value was adjusted by -1.16485e-15 + Note: The mean value was adjusted by -1.33296e-15 Velocity L2 Error: 5.33247e-05 Pressure L2 Error: 0.00102063 diff --git a/tests/lac/gmres_reorthogonalize_01.cc b/tests/lac/gmres_reorthogonalize_01.cc index 077cae7298..78a5b995a0 100644 --- a/tests/lac/gmres_reorthogonalize_01.cc +++ b/tests/lac/gmres_reorthogonalize_01.cc @@ -69,7 +69,7 @@ test(unsigned int variant, unsigned int min_convergence_steps) SolverControl control(1000, 1e2 * std::numeric_limits::epsilon()); typename SolverGMRES>::AdditionalData data; - data.max_n_tmp_vectors = 80; + data.max_basis_size = 80; SolverGMRES> solver(control, data); auto print_re_orthogonalization = [](int accumulated_iterations) { diff --git a/tests/lac/gmres_reorthogonalize_02.cc b/tests/lac/gmres_reorthogonalize_02.cc index dd44f34476..572485a7e6 100644 --- a/tests/lac/gmres_reorthogonalize_02.cc +++ b/tests/lac/gmres_reorthogonalize_02.cc @@ -43,7 +43,7 @@ test() SolverControl control(1000, 1e3 * std::numeric_limits::epsilon()); typename SolverGMRES>::AdditionalData data; - data.max_n_tmp_vectors = 202; + data.max_basis_size = 200; SolverGMRES> solver(control, data); auto print_re_orthogonalization = [](int accumulated_iterations) { diff --git a/tests/lac/gmres_reorthogonalize_03.cc b/tests/lac/gmres_reorthogonalize_03.cc index 399cc5c789..bccc1a3524 100644 --- a/tests/lac/gmres_reorthogonalize_03.cc +++ b/tests/lac/gmres_reorthogonalize_03.cc @@ -67,7 +67,7 @@ test(unsigned int variant) SolverControl control(1000, 1e2 * std::numeric_limits::epsilon()); typename SolverGMRES>::AdditionalData data; - data.max_n_tmp_vectors = 80; + data.max_basis_size = 80; data.force_re_orthogonalization = true; SolverGMRES> solver(control, data); diff --git a/tests/lac/gmres_reorthogonalize_04.cc b/tests/lac/gmres_reorthogonalize_04.cc index 18cb893b5d..fd41939e92 100644 --- a/tests/lac/gmres_reorthogonalize_04.cc +++ b/tests/lac/gmres_reorthogonalize_04.cc @@ -73,7 +73,7 @@ test(const unsigned int n_expected_steps) // of re-orthogonalization) SolverControl control(1000, 1e3 * std::numeric_limits::epsilon()); typename SolverGMRES>::AdditionalData data; - data.max_n_tmp_vectors = 202; + data.max_basis_size = 200; SolverGMRES> solver(control, data); auto print_re_orthogonalization = [](int accumulated_iterations) { diff --git a/tests/lac/gmres_reorthogonalize_05.cc b/tests/lac/gmres_reorthogonalize_05.cc index 1743b1910c..37774657d8 100644 --- a/tests/lac/gmres_reorthogonalize_05.cc +++ b/tests/lac/gmres_reorthogonalize_05.cc @@ -70,7 +70,7 @@ test() SolverControl control(1000, 1e2 * std::numeric_limits::epsilon()); typename SolverGMRES>::AdditionalData data; - data.max_n_tmp_vectors = 202; + data.max_basis_size = 200; data.force_re_orthogonalization = true; SolverGMRES> solver(control, data); diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index 5a4abcd39b..800b4bffdb 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -89,9 +89,9 @@ main() SolverControl control(100, 1.e-3); SolverControl verbose_control(100, 1.e-3, true); SolverCG<> cg(control, mem); - SolverGMRES<>::AdditionalData data1(8); + SolverGMRES<>::AdditionalData data1(6); SolverGMRES<> gmres(control, mem, data1); - SolverGMRES<>::AdditionalData data2(8, true); + SolverGMRES<>::AdditionalData data2(6, true); SolverGMRES<> gmresright(control, mem, data2); SolverMinRes<> minres(control, mem); SolverBicgstab<> bicgstab(control, mem); @@ -99,7 +99,7 @@ main() SolverQMRS<> qmrs(control, mem); SolverFIRE<> fire(control, mem); - SolverGMRES<>::AdditionalData data3(8); + SolverGMRES<>::AdditionalData data3(6); data3.orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt; SolverGMRES<> gmresclassical(control, mem, data3); diff --git a/tests/lac/solver.output b/tests/lac/solver.output index 653812ce53..4dd322bb37 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -110,13 +110,13 @@ DEAL:no-fail:Bicgstab::Failure step 10 value 0.001961 DEAL:no-fail::Exception: SolverControl::NoConvergence(state.last_step, state.last_residual) DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 -DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res) DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 -DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res) DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 -DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res) DEAL:no-fail:SQMR::Starting value 11.00 DEAL:no-fail:SQMR::Failure step 10 value 0.4215 DEAL:no-fail::Exception: SolverControl::NoConvergence(step, state.last_residual) diff --git a/tests/lac/solver_gmres_01.cc b/tests/lac/solver_gmres_01.cc index e140dc5c2f..f788c3e658 100644 --- a/tests/lac/solver_gmres_01.cc +++ b/tests/lac/solver_gmres_01.cc @@ -105,7 +105,7 @@ main() deallog << "Solve with PreconditionIdentity: " << std::endl; SolverControl control(40, 1e-4); SolverGMRES>::AdditionalData - data3(8); + data3(6); data3.orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt; SolverGMRES> solver(control, @@ -127,7 +127,7 @@ main() deallog << "Solve with PreconditionIdentity: " << std::endl; SolverControl control(40, 1e-4); SolverGMRES>::AdditionalData - data3(8); + data3(6); data3.orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt; SolverGMRES> solver(control, diff --git a/tests/lac/solver_gmres_02.cc b/tests/lac/solver_gmres_02.cc index bcf6aee22f..3dd6103bc9 100644 --- a/tests/lac/solver_gmres_02.cc +++ b/tests/lac/solver_gmres_02.cc @@ -103,7 +103,7 @@ main() { deallog << "Solve with PreconditionIdentity: " << std::endl; SolverControl control(40, 1e-4); - SolverGMRES>::AdditionalData data3(8); + SolverGMRES>::AdditionalData data3(6); data3.orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt; SolverGMRES> solver(control, data3); @@ -123,7 +123,7 @@ main() deallog << "Solve with PreconditionIdentity: " << std::endl; SolverControl control(40, 1e-4); - SolverGMRES>::AdditionalData data3(8); + SolverGMRES>::AdditionalData data3(6); data3.orthogonalization_strategy = LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt; SolverGMRES> solver(control, data3); diff --git a/tests/lac/solver_monitor.cc b/tests/lac/solver_monitor.cc index 4a828b0aec..22d0c8ed46 100644 --- a/tests/lac/solver_monitor.cc +++ b/tests/lac/solver_monitor.cc @@ -97,7 +97,7 @@ main() SolverGMRES<> gmres(control, mem, - SolverGMRES<>::AdditionalData(/*max_vecs=*/8)); + SolverGMRES<>::AdditionalData(/*basis_size=*/6)); gmres.connect(&monitor_norm); gmres.connect(&monitor_mean); diff --git a/tests/lac/solver_selector_01.output b/tests/lac/solver_selector_01.output index 8f2bd5c749..d936b867f3 100644 --- a/tests/lac/solver_selector_01.output +++ b/tests/lac/solver_selector_01.output @@ -13,6 +13,6 @@ DEAL:cg::Convergence step 39 value 7.335e-08 DEAL:Bicgstab::Starting value 36.00 DEAL:Bicgstab::Convergence step 28 value 9.087e-08 DEAL:GMRES::Starting value 34.45 -DEAL:GMRES::Convergence step 39 value 6.627e-08 +DEAL:GMRES::Convergence step 40 value 9.928e-08 DEAL:FGMRES::Starting value 36.00 -DEAL:FGMRES::Convergence step 39 value 7.890e-08 +DEAL:FGMRES::Convergence step 40 value 9.803e-08 diff --git a/tests/lac/solver_selector_03.output b/tests/lac/solver_selector_03.output index a716a5ebb1..d4c88c708a 100644 --- a/tests/lac/solver_selector_03.output +++ b/tests/lac/solver_selector_03.output @@ -3,4 +3,4 @@ DEAL::Size 37 Unknowns 1296 DEAL:cg::Starting value 36.00 DEAL:cg::Convergence step 39 value 7.335e-08 DEAL:GMRES::Starting value 34.45 -DEAL:GMRES::Convergence step 39 value 6.627e-08 +DEAL:GMRES::Convergence step 40 value 9.928e-08 diff --git a/tests/lac/solver_signals.cc b/tests/lac/solver_signals.cc index 48089657bf..761e9a5c15 100644 --- a/tests/lac/solver_signals.cc +++ b/tests/lac/solver_signals.cc @@ -154,7 +154,8 @@ main() solver_cg.solve(A, u, f, PreconditionIdentity()); u = 0; - SolverGMRES<> solver_gmres(solver_control); + SolverGMRES<> solver_gmres(solver_control, + SolverGMRES<>::AdditionalData(28)); // Attach all possible slots. solver_gmres.connect_condition_number_slot( std::bind(output_double_number, diff --git a/tests/lac/solver_signals.with_lapack=true.output b/tests/lac/solver_signals.with_lapack=true.output index f597ef56d0..ac1b0228b1 100644 --- a/tests/lac/solver_signals.with_lapack=true.output +++ b/tests/lac/solver_signals.with_lapack=true.output @@ -128,36 +128,35 @@ DEAL:GMRES::Convergence step 65 value 0.0009799 DEAL:GMRES::Eigenvalues: (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000) DEAL:GMRES::Final condition number: 288.2 DEAL:GMRES::Hessenberg matrix: -DEAL:GMRES::0.1982 0.8628 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0.8628 6.098 1.971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 1.971 4.008 1.721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 1.721 3.583 2.004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 2.004 4.008 1.910 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 1.910 3.325 1.950 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 1.950 4.036 1.744 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 1.744 3.951 1.992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 1.992 4.080 1.994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 1.920 4.007 2.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 2.002 3.981 2.022 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 2.022 4.028 1.908 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 1.908 4.238 1.925 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 1.925 4.011 2.042 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 2.042 3.827 2.111 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.111 3.742 2.046 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 4.070 1.879 0 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.879 4.168 1.933 0 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.933 4.013 1.886 0 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.886 3.963 1.957 0 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.957 3.818 2.033 0 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.033 4.006 2.201 0 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.201 3.944 1.862 0 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.862 4.039 1.827 0 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.827 3.704 1.844 0 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.844 4.073 2.046 0 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 3.592 1.647 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.647 3.431 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.155 0 -DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0.1982 0.8628 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0.8628 6.098 1.971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 1.971 4.008 1.721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 1.721 3.583 2.004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 2.004 4.008 1.910 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 1.910 3.325 1.950 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 1.950 4.036 1.744 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 1.744 3.951 1.992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 1.992 4.080 1.994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 1.920 4.007 2.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 2.002 3.981 2.022 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 2.022 4.028 1.908 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 1.908 4.238 1.925 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 1.925 4.011 2.042 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 2.042 3.827 2.111 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.111 3.742 2.046 0 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 4.070 1.879 0 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.879 4.168 1.933 0 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.933 4.013 1.886 0 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.886 3.963 1.957 0 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.957 3.818 2.033 0 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.033 4.006 2.201 0 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.201 3.944 1.862 0 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.862 4.039 1.827 0 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.827 3.704 1.844 0 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.844 4.073 2.046 0 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 3.592 1.647 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.647 3.431 +DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.155 DEAL:GMRES::Final Eigenvalues: (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000) DEAL:GMRES::Final condition number: 288.2 DEAL:GMRES::Arnoldi vectors norms: diff --git a/tests/matrix_free/stokes_computation.cc b/tests/matrix_free/stokes_computation.cc index 51b6d4c908..c685e92c4d 100644 --- a/tests/matrix_free/stokes_computation.cc +++ b/tests/matrix_free/stokes_computation.cc @@ -1386,13 +1386,19 @@ namespace StokesClass solution.block(block_vel) = distributed_stokes_solution.block(block_vel); solution.block(block_p) = distributed_stokes_solution.block(block_p); + LinearAlgebra::distributed::BlockVector vec( + distributed_stokes_solution); + stokes_matrix.vmult(vec, distributed_stokes_solution); + vec -= distributed_stokes_rhs; + deallog << "Solved-in " << (solver_control_cheap.last_step() != numbers::invalid_unsigned_int ? solver_control_cheap.last_step() : 0) << " iterations, final residual: " - << solver_control_cheap.last_value() << std::endl; + << solver_control_cheap.last_value() << " and true residual " + << vec.l2_norm() << std::endl; } diff --git a/tests/matrix_free/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output b/tests/matrix_free/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output index 1f2687b42c..ba7275b72b 100644 --- a/tests/matrix_free/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output +++ b/tests/matrix_free/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output @@ -4,16 +4,16 @@ DEAL:2d::n_sinker: 4 max/min viscosity ratio: 1000.00 DEAL:2d:: DEAL:2d::Number of active cells: 1024 (on 6 levels) DEAL:2d::Number of degrees of freedom: 9539 (8450+1089) -DEAL:2d::Solved-in 69 iterations, final residual: 5.58979e-08 +DEAL:2d::Solved-in 64 iterations, final residual: 5.31067e-08 and true residual 5.31067e-08 DEAL:2d:: DEAL:2d::Number of active cells: 4096 (on 7 levels) DEAL:2d::Number of degrees of freedom: 37507 (33282+4225) -DEAL:2d::Solved-in 69 iterations, final residual: 2.58019e-08 +DEAL:2d::Solved-in 64 iterations, final residual: 2.55967e-08 and true residual 2.55967e-08 DEAL:2d:: DEAL:3d::Sinker problem in 3D. DEAL:3d::n_sinker: 4 max/min viscosity ratio: 1000.00 DEAL:3d:: DEAL:3d::Number of active cells: 512 (on 4 levels) DEAL:3d::Number of degrees of freedom: 15468 (14739+729) -DEAL:3d::Solved-in 62 iterations, final residual: 1.99453e-08 +DEAL:3d::Solved-in 58 iterations, final residual: 2.08315e-08 and true residual 2.08315e-08 DEAL:3d:: diff --git a/tests/multigrid-global-coarsening/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output b/tests/multigrid-global-coarsening/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output index 1f2687b42c..46baf3fb41 100644 --- a/tests/multigrid-global-coarsening/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output +++ b/tests/multigrid-global-coarsening/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output @@ -4,16 +4,16 @@ DEAL:2d::n_sinker: 4 max/min viscosity ratio: 1000.00 DEAL:2d:: DEAL:2d::Number of active cells: 1024 (on 6 levels) DEAL:2d::Number of degrees of freedom: 9539 (8450+1089) -DEAL:2d::Solved-in 69 iterations, final residual: 5.58979e-08 +DEAL:2d::Solved-in 64 iterations, final residual: 5.31067e-08 DEAL:2d:: DEAL:2d::Number of active cells: 4096 (on 7 levels) DEAL:2d::Number of degrees of freedom: 37507 (33282+4225) -DEAL:2d::Solved-in 69 iterations, final residual: 2.58019e-08 +DEAL:2d::Solved-in 64 iterations, final residual: 2.55967e-08 DEAL:2d:: DEAL:3d::Sinker problem in 3D. DEAL:3d::n_sinker: 4 max/min viscosity ratio: 1000.00 DEAL:3d:: DEAL:3d::Number of active cells: 512 (on 4 levels) DEAL:3d::Number of degrees of freedom: 15468 (14739+729) -DEAL:3d::Solved-in 62 iterations, final residual: 1.99453e-08 +DEAL:3d::Solved-in 58 iterations, final residual: 2.08315e-08 DEAL:3d:: diff --git a/tests/multigrid/step-39-03.output b/tests/multigrid/step-39-03.output index 002798281b..e271e6ef08 100644 --- a/tests/multigrid/step-39-03.output +++ b/tests/multigrid/step-39-03.output @@ -8,7 +8,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 9.41508 -DEAL:GMRES::Convergence step 10 value 7.45878e-13 +DEAL:GMRES::Convergence step 10 value 7.45877e-13 DEAL::energy-error: 0.439211 DEAL::L2-error: 0.0109342 DEAL::Estimate 0.979555 @@ -22,7 +22,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 10.0041 -DEAL:GMRES::Convergence step 21 value 1.79039e-13 +DEAL:GMRES::Convergence step 21 value 1.79024e-13 DEAL::energy-error: 0.332582 DEAL::L2-error: 0.00548996 DEAL::Estimate 0.838060 @@ -36,7 +36,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 10.8018 -DEAL:GMRES::Convergence step 28 value 2.71098e-13 +DEAL:GMRES::Convergence step 28 value 2.71118e-13 DEAL::energy-error: 0.237705 DEAL::L2-error: 0.00250869 DEAL::Estimate 0.611449 @@ -50,7 +50,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 13.4185 -DEAL:GMRES::Convergence step 32 value 6.11150e-13 +DEAL:GMRES::Convergence step 31 value 9.79574e-13 DEAL::energy-error: 0.170175 DEAL::L2-error: 0.00120328 DEAL::Estimate 0.448676 @@ -64,7 +64,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 16.6112 -DEAL:GMRES::Convergence step 32 value 7.81623e-13 +DEAL:GMRES::Convergence step 32 value 7.99534e-13 DEAL::energy-error: 0.121878 DEAL::L2-error: 0.000591777 DEAL::Estimate 0.327296 @@ -78,7 +78,7 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:GMRES::Starting value 19.3095 -DEAL:GMRES::Convergence step 31 value 8.40512e-13 +DEAL:GMRES::Convergence step 31 value 9.61737e-13 DEAL::energy-error: 0.0863235 DEAL::L2-error: 0.000292807 DEAL::Estimate 0.232279 diff --git a/tests/non_matching/step-70.with_p4est=true.with_petsc=true.mpirun=2.output b/tests/non_matching/step-70.with_p4est=true.with_petsc=true.mpirun=2.output index 72980d7cfd..23504a72e7 100644 --- a/tests/non_matching/step-70.with_p4est=true.with_petsc=true.mpirun=2.output +++ b/tests/non_matching/step-70.with_p4est=true.with_petsc=true.mpirun=2.output @@ -6,29 +6,29 @@ DEAL:0:: Number of degrees of freedom: 9539 (8450+1089 -- 0+0) DEAL:0::Tracer particles: 1504 DEAL:0::Solid particles: 9216 DEAL:0:FGMRES::Starting value 30.9301 -DEAL:0:FGMRES::Convergence step 147 value 2.63430e-09 -DEAL:0:: Solved in 147 iterations. +DEAL:0:FGMRES::Convergence step 148 value 2.35433e-09 +DEAL:0:: Solved in 148 iterations. DEAL:0:: Number of degrees of freedom: 9845 (8722+1123 -- 9216+1504) DEAL:0::Cycle 1: DEAL:0::Time : 0.00250000, time step: 0.00250000 DEAL:0:FGMRES::Starting value 2.92442 -DEAL:0:FGMRES::Convergence step 141 value 5.01043e-09 -DEAL:0:: Solved in 141 iterations. +DEAL:0:FGMRES::Convergence step 146 value 5.16889e-09 +DEAL:0:: Solved in 146 iterations. DEAL:0::Cycle 2: DEAL:0::Time : 0.00500000, time step: 0.00250000 DEAL:0:FGMRES::Starting value 0.637099 -DEAL:0:FGMRES::Convergence step 121 value 4.79602e-09 -DEAL:0:: Solved in 121 iterations. +DEAL:0:FGMRES::Convergence step 117 value 5.22134e-09 +DEAL:0:: Solved in 117 iterations. DEAL:0::Cycle 3: DEAL:0::Time : 0.00750000, time step: 0.00250000 DEAL:0:FGMRES::Starting value 0.649739 -DEAL:0:FGMRES::Convergence step 126 value 5.21191e-09 -DEAL:0:: Solved in 126 iterations. +DEAL:0:FGMRES::Convergence step 125 value 5.07375e-09 +DEAL:0:: Solved in 125 iterations. DEAL:0::Cycle 4: DEAL:0::Time : 0.0100000, time step: 0.00250000 DEAL:0:FGMRES::Starting value 0.744112 -DEAL:0:FGMRES::Convergence step 131 value 4.88195e-09 -DEAL:0:: Solved in 131 iterations. +DEAL:0:FGMRES::Convergence step 130 value 5.49751e-09 +DEAL:0:: Solved in 130 iterations. DEAL:1::Running StokesImmersedProblem<2> using Trilinos. DEAL:1::Cycle 0: @@ -37,27 +37,27 @@ DEAL:1:: Number of degrees of freedom: 9539 (8450+1089 -- 0+0) DEAL:1::Tracer particles: 1504 DEAL:1::Solid particles: 9216 DEAL:1:FGMRES::Starting value 30.9301 -DEAL:1:FGMRES::Convergence step 147 value 2.63430e-09 -DEAL:1:: Solved in 147 iterations. +DEAL:1:FGMRES::Convergence step 148 value 2.35433e-09 +DEAL:1:: Solved in 148 iterations. DEAL:1:: Number of degrees of freedom: 9845 (8722+1123 -- 9216+1504) DEAL:1::Cycle 1: DEAL:1::Time : 0.00250000, time step: 0.00250000 DEAL:1:FGMRES::Starting value 2.92442 -DEAL:1:FGMRES::Convergence step 141 value 5.01043e-09 -DEAL:1:: Solved in 141 iterations. +DEAL:1:FGMRES::Convergence step 146 value 5.16889e-09 +DEAL:1:: Solved in 146 iterations. DEAL:1::Cycle 2: DEAL:1::Time : 0.00500000, time step: 0.00250000 DEAL:1:FGMRES::Starting value 0.637099 -DEAL:1:FGMRES::Convergence step 121 value 4.79602e-09 -DEAL:1:: Solved in 121 iterations. +DEAL:1:FGMRES::Convergence step 117 value 5.22134e-09 +DEAL:1:: Solved in 117 iterations. DEAL:1::Cycle 3: DEAL:1::Time : 0.00750000, time step: 0.00250000 DEAL:1:FGMRES::Starting value 0.649739 -DEAL:1:FGMRES::Convergence step 126 value 5.21191e-09 -DEAL:1:: Solved in 126 iterations. +DEAL:1:FGMRES::Convergence step 125 value 5.07375e-09 +DEAL:1:: Solved in 125 iterations. DEAL:1::Cycle 4: DEAL:1::Time : 0.0100000, time step: 0.00250000 DEAL:1:FGMRES::Starting value 0.744112 -DEAL:1:FGMRES::Convergence step 131 value 4.88195e-09 -DEAL:1:: Solved in 131 iterations. +DEAL:1:FGMRES::Convergence step 130 value 5.49751e-09 +DEAL:1:: Solved in 130 iterations. diff --git a/tests/trilinos/deal_solver_03.cc b/tests/trilinos/deal_solver_03.cc index 8570662278..ffa1065eb8 100644 --- a/tests/trilinos/deal_solver_03.cc +++ b/tests/trilinos/deal_solver_03.cc @@ -82,7 +82,7 @@ main(int argc, char **argv) check_solver_within_range(solver.solve(A, u, f, preconditioner), control.last_step(), - 74, + 71, 76); } } diff --git a/tests/trilinos/deal_solver_03.output b/tests/trilinos/deal_solver_03.output index 5b6f34323b..3e4091158b 100644 --- a/tests/trilinos/deal_solver_03.output +++ b/tests/trilinos/deal_solver_03.output @@ -1,4 +1,4 @@ DEAL::Size 32 Unknowns 961 DEAL::Solver type: N6dealii11SolverGMRESINS_16TrilinosWrappers3MPI6VectorEEE -DEAL::Solver stopped within 74 - 76 iterations +DEAL::Solver stopped within 71 - 76 iterations diff --git a/tests/trilinos_tpetra/deal_solver_03.cc b/tests/trilinos_tpetra/deal_solver_03.cc index 56085099d1..1b926dae5c 100644 --- a/tests/trilinos_tpetra/deal_solver_03.cc +++ b/tests/trilinos_tpetra/deal_solver_03.cc @@ -75,8 +75,11 @@ main(int argc, char **argv) u.compress(VectorOperation::insert); GrowingVectorMemory> mem; + SolverGMRES>::AdditionalData + sol_data(28); SolverGMRES> solver(control, - mem); + mem, + sol_data); PreconditionIdentity preconditioner; deallog << "Solver type: " << typeid(solver).name() << std::endl;