From d56e6e10aaaccdc9ca8cdad887e4e63670b9e1f0 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 4 Dec 2020 12:26:50 +0100 Subject: [PATCH] Reshuffle code in SoverCG --- include/deal.II/lac/solver_cg.h | 128 ++++++++++++++++---------------- 1 file changed, 65 insertions(+), 63 deletions(-) diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index 5aaf20a3e6..e2e420eed1 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -341,28 +341,26 @@ SolverCG::solve(const MatrixType & A, !condition_number_signal.empty() || !all_condition_numbers_signal.empty() || !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty(); - // vectors used for eigenvalue - // computations + // vectors used for eigenvalue computations std::vector diagonal; std::vector offdiagonal; - int it = 0; - double res = -std::numeric_limits::max(); - typename VectorType::value_type eigen_beta_alpha = 0; - // resize the vectors, but do not set - // the values since they'd be overwritten + // resize the vectors, but do not set the values since they'd be overwritten // soon anyway. g.reinit(x, true); d.reinit(x, true); h.reinit(x, true); - number gh, beta; + int it = 0; + number gh = number(); + number beta = number(); + number alpha = number(); + number old_alpha = number(); - // compute residual. if vector is - // zero, then short-circuit the - // full computation + // compute residual. if vector is zero, then short-circuit the full + // computation if (!x.all_zero()) { A.vmult(g, x); @@ -370,32 +368,56 @@ SolverCG::solve(const MatrixType & A, } else g.equ(-1., b); - res = g.l2_norm(); - conv = this->iteration_status(0, res, x); + double res = g.l2_norm(); + conv = this->iteration_status(0, res, x); if (conv != SolverControl::iterate) return; - if (std::is_same::value == false) - { - preconditioner.vmult(h, g); - - d.equ(-1., h); - - gh = g * h; - } - else - { - d.equ(-1., g); - gh = res * res; - } - while (conv == SolverControl::iterate) { it++; + old_alpha = alpha; + + if (it > 1) + { + if (std::is_same::value == + false) + { + preconditioner.vmult(h, g); + beta = gh; + Assert(std::abs(beta) != 0., ExcDivideByZero()); + gh = g * h; + beta = gh / beta; + d.sadd(beta, -1., h); + } + else + { + beta = gh; + gh = res * res; + beta = gh / beta; + d.sadd(beta, -1., g); + } + } + else + { + if (std::is_same::value == + false) + { + preconditioner.vmult(h, g); + d.equ(-1., h); + gh = g * h; + } + else + { + d.equ(-1., g); + gh = res * res; + } + } + A.vmult(h, d); - number alpha = d * h; + alpha = d * h; Assert(std::abs(alpha) != 0., ExcDivideByZero()); alpha = gh / alpha; @@ -404,44 +426,24 @@ SolverCG::solve(const MatrixType & A, print_vectors(it, x, g, d); - conv = this->iteration_status(it, res, x); - if (conv != SolverControl::iterate) - break; - - if (std::is_same::value == - false) + if (it > 1) { - preconditioner.vmult(h, g); - - beta = gh; - Assert(std::abs(beta) != 0., ExcDivideByZero()); - gh = g * h; - beta = gh / beta; - d.sadd(beta, -1., h); - } - else - { - beta = gh; - gh = res * res; - beta = gh / beta; - d.sadd(beta, -1., g); + this->coefficients_signal(old_alpha, beta); + // set up the vectors containing the diagonal and the off diagonal of + // the projected matrix. + if (do_eigenvalues) + { + diagonal.push_back(number(1.) / old_alpha + eigen_beta_alpha); + eigen_beta_alpha = beta / old_alpha; + offdiagonal.push_back(std::sqrt(beta) / old_alpha); + } + compute_eigs_and_cond(diagonal, + offdiagonal, + all_eigenvalues_signal, + all_condition_numbers_signal); } - this->coefficients_signal(alpha, beta); - // set up the vectors - // containing the diagonal - // and the off diagonal of - // the projected matrix. - if (do_eigenvalues) - { - diagonal.push_back(number(1.) / alpha + eigen_beta_alpha); - eigen_beta_alpha = beta / alpha; - offdiagonal.push_back(std::sqrt(beta) / alpha); - } - compute_eigs_and_cond(diagonal, - offdiagonal, - all_eigenvalues_signal, - all_condition_numbers_signal); + conv = this->iteration_status(it, res, x); } compute_eigs_and_cond(diagonal, -- 2.39.5