From ae72026ef303315ff3ca656ed5c588be1d35b65c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@mi.fu-berlin.de> Date: Fri, 28 Mar 2014 12:25:07 +0000 Subject: [PATCH] Simplify branch structure in damping code [[Imported from SVN: r13079]] --- .../blocknonlineartnnmgproblem.hh | 156 +++++++++--------- 1 file changed, 78 insertions(+), 78 deletions(-) diff --git a/dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh b/dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh index abef0a8..26446e6 100644 --- a/dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh +++ b/dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh @@ -447,110 +447,110 @@ class BlockNonlinearTNNMGProblem */ double computeDampingParameter(const VectorType& u, const VectorType& projected_v) const { + if (not(linesearch.enable_flag)) + return 1.0; + double alpha = 1.0; - if (linesearch.enable_flag) + // normalize vector v to increase stability + VectorType v = projected_v; + double scalingFactor = 1.0; + if (linesearch.normed_flag) { - // normalize vector v to increase stability - VectorType v = projected_v; - double scalingFactor = 1.0; - if (linesearch.normed_flag) - { - scalingFactor = v.two_norm(); - if (scalingFactor > 0.0) - v /= scalingFactor; - else - return 0.0; - } + scalingFactor = v.two_norm(); + if (scalingFactor > 0.0) + v /= scalingFactor; + else + return 0.0; + } - // initialize alpha by 1 (up to rescaling) - alpha = scalingFactor; + // initialize alpha by 1 (up to rescaling) + alpha = scalingFactor; - // setup directional convex functional - double Avv; - double resv; - { - VectorType temp(u.size()); + // setup directional convex functional + double Avv; + double resv; + { + VectorType temp(u.size()); - // compute quadratic part a*<Av,v> - problem_.A.mv(v, temp); - Avv = problem_.a*(temp*v); + // compute quadratic part a*<Av,v> + problem_.A.mv(v, temp); + Avv = problem_.a*(temp*v); - // compute linear part <f-a*Au,v> - temp = problem_.f; - problem_.A.template usmv<VectorType, VectorType>(-problem_.a, u, temp); - resv = temp*v; + // compute linear part <f-a*Au,v> + temp = problem_.f; + problem_.A.template usmv<VectorType, VectorType>(-problem_.a, u, temp); + resv = temp*v; - // only touch the lowRankFactor if its coefficient is not zero; it might not be initialized otherwise - if (problem_.am != 0.0) - { + // only touch the lowRankFactor if its coefficient is not zero; it might not be initialized otherwise + if (problem_.am != 0.0) + { - // compute quadratic part from low rank factor am*<Am'*Am v,v> = am*<Am v, Am v> - temp.resize(problem_.lowRankFactor_.N()); - problem_.lowRankFactor_.mv(v,temp); - Avv += problem_.am * (temp * temp); + // compute quadratic part from low rank factor am*<Am'*Am v,v> = am*<Am v, Am v> + temp.resize(problem_.lowRankFactor_.N()); + problem_.lowRankFactor_.mv(v,temp); + Avv += problem_.am * (temp * temp); - // compute linear part -am*<Am'*Am u,v> = -am*a<Am u, Am v> - VectorType temp2(problem_.lowRankFactor_.N()); - problem_.lowRankFactor_.mv(u,temp2); - resv -= problem_.am * (temp * temp2); - } + // compute linear part -am*<Am'*Am u,v> = -am*a<Am u, Am v> + VectorType temp2(problem_.lowRankFactor_.N()); + problem_.lowRankFactor_.mv(u,temp2); + resv -= problem_.am * (temp * temp2); } - DirectionalConvexFunction<NonlinearityType> psi(Avv, resv, problem_.phi, u, v); - - // if neccessary compute energy before coarse correction - double oldEnergy; - if (linesearch.descentonly_flag or linesearch.try_one_flag) - oldEnergy = psi(0); + } + DirectionalConvexFunction<NonlinearityType> psi(Avv, resv, problem_.phi, u, v); - // if neccessary compute energy at alpha - double energy=0; - if (linesearch.descentonly_flag or linesearch.try_one_flag or config.get("logging.energy", true)) - energy = psi(alpha); + // if neccessary compute energy before coarse correction + double oldEnergy; + if (linesearch.descentonly_flag or linesearch.try_one_flag) + oldEnergy = psi(0); - if (linesearch.descentonly_flag) - { - while (energy>oldEnergy) - { - alpha *= linesearch.descentonly_reduction; - oldEnergy = psi(0); - energy = psi(alpha); - } + // if neccessary compute energy at alpha + double energy=0; + if (linesearch.descentonly_flag or linesearch.try_one_flag or config.get("logging.energy", true)) + energy = psi(alpha); - } - else + if (linesearch.descentonly_flag) + { + while (energy>oldEnergy) { - if (not(linesearch.try_one_flag) or (energy > oldEnergy)) - { - int bisectionsteps = 0; - Bisection bisection(0.0, linesearch.acceptance, linesearch.tol, linesearch.fast_quadratic_minimize_flag); - alpha = bisection.minimize(psi, scalingFactor, 0.0, bisectionsteps); - if (config.get("logging.energy", true)) - energy = psi(alpha); - } + alpha *= linesearch.descentonly_reduction; + oldEnergy = psi(0); + energy = psi(alpha); } - // restrict alpha to [0,1] - if (linesearch.damponly_flag) + } + else + { + if (not(linesearch.try_one_flag) or (energy > oldEnergy)) { - if (alpha < 0.0) - alpha = 0.0; - if (alpha > scalingFactor) - alpha = scalingFactor; + int bisectionsteps = 0; + Bisection bisection(0.0, linesearch.acceptance, linesearch.tol, linesearch.fast_quadratic_minimize_flag); + alpha = bisection.minimize(psi, scalingFactor, 0.0, bisectionsteps); if (config.get("logging.energy", true)) energy = psi(alpha); } + } + // restrict alpha to [0,1] + if (linesearch.damponly_flag) + { + if (alpha < 0.0) + alpha = 0.0; + if (alpha > scalingFactor) + alpha = scalingFactor; if (config.get("logging.energy", true)) - { - outStream.setf(std::ios::scientific); - outStream << std::setw(15) << energy; - } + energy = psi(alpha); + } - // rescale alpha - alpha = alpha/scalingFactor; + if (config.get("logging.energy", true)) + { + outStream.setf(std::ios::scientific); + outStream << std::setw(15) << energy; } + // rescale alpha + alpha = alpha/scalingFactor; + return alpha; }; -- GitLab