Skip to content
Snippets Groups Projects
Commit ae72026e authored by graeser's avatar graeser Committed by graeser
Browse files

Simplify branch structure in damping code 

[[Imported from SVN: r13079]]
parent 7135d642
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
+ 78
− 78
View file @ ae72026e
......@@ -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;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment