00001
00002
00003
00004
00005 #ifndef PLAYA_NEWTON_ARMIJO_SOLVER_IMPL_HPP
00006 #define PLAYA_NEWTON_ARMIJO_SOLVER_IMPL_HPP
00007
00008 #include "PlayaNewtonArmijoSolverDecl.hpp"
00009 #include "PlayaNonlinearOperator.hpp"
00010 #include "PlayaTabs.hpp"
00011 #include "PlayaOut.hpp"
00012 #include "Teuchos_ParameterList.hpp"
00013
00014 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00015 #include "PlayaLinearCombinationImpl.hpp"
00016 #include "PlayaLinearSolverImpl.hpp"
00017 #include "PlayaLinearOperatorImpl.hpp"
00018 #endif
00019
00020
00021 using std::setw;
00022
00023 namespace Playa
00024 {
00025
00026
00027
00028 template <class Scalar> inline
00029 NewtonArmijoSolver<Scalar>::NewtonArmijoSolver(
00030 const ParameterList& params,
00031 const LinearSolver<Scalar>& linSolver)
00032 : NonlinearSolverBase<Scalar>(params),
00033 linSolver_(linSolver),
00034 tauR_(10.0*Teuchos::ScalarTraits<Scalar>::eps()),
00035 tauA_(10.0*Teuchos::ScalarTraits<Scalar>::eps()),
00036 alpha_(1.0e-4),
00037 stepReduction_(0.5),
00038 maxIters_(20),
00039 maxLineSearch_(20),
00040 verb_(0)
00041 {
00042 if (params.isParameter("Tau Relative")) tauR_ = params.get<Scalar>("Tau Relative");
00043 if (params.isParameter("Tau Absolute")) tauA_ = params.get<Scalar>("Tau Absolute");
00044 if (params.isParameter("Alpha")) alpha_ = params.get<Scalar>("Alpha");
00045 if (params.isParameter("Step Reduction")) stepReduction_ = params.get<Scalar>("Step Reduction");
00046 if (params.isParameter("Max Iterations")) maxIters_ = params.get<int>("Max Iterations");
00047 if (params.isParameter("Max Backtracks")) maxLineSearch_ = params.get<int>("Max Backtracks");
00048 if (params.isParameter("Verbosity")) verb_ = params.get<int>("Verbosity");
00049 }
00050
00051 template <class Scalar> inline
00052 SolverState<Scalar> NewtonArmijoSolver<Scalar>::solve(const NonlinearOperator<Scalar>& F,
00053 Vector<Scalar>& soln) const
00054 {
00055 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00056 typedef typename Teuchos::ScalarTraits<Scalar> ST;
00057
00058 Tabs tab0(0);
00059 PLAYA_MSG1(verb_, tab0 << " begin Playa::NewtonArmijoSolver::solve()");
00060
00061 soln = F.getInitialGuess().copy();
00062 Vector<Scalar> newtonStep = soln.copy();
00063
00064 F.setEvalPt(soln);
00065 Vector<Scalar> resid = F.getFunctionValue();
00066
00067 ScalarMag r0 = resid.norm2();
00068 ScalarMag normF0 = r0;
00069
00070 for (int i=0; i<maxIters_; i++)
00071 {
00072 Tabs tab1;
00073 PLAYA_MSG2(verb_, tab1 << "Newton iter #" << setw(6) << i << " |F|=" << setw(12) << normF0 << " |F|/|F0|="
00074 << setw(12) << normF0/r0);
00075
00076 if (normF0 < r0*tauR_ + tauA_)
00077 {
00078 PLAYA_MSG3(verb_, tab1 << "|F|=" << setw(12) << normF0);
00079 PLAYA_MSG3(verb_, tab1 << "Relative tolerance tauR=" << setw(12) << tauR_);
00080 PLAYA_MSG3(verb_, tab1 << "Absolute tolerance tauA=" << setw(12) << tauA_);
00081 PLAYA_MSG3(verb_, tab1 << " F0*tauR+tauA=" << setw(12) << r0*tauR_ + tauA_);
00082 PLAYA_MSG2(verb_, tab1 << "converged!");
00083 PLAYA_MSG1(verb_, tab0 << " done Playa::NewtonArmijoSolver::solve()");
00084 soln = F.currentEvalPt().copy();
00085 return SolverState<Scalar>(SolveConverged, "NewtonArmijoSolver::solve converged",
00086 i, normF0);
00087 }
00088 LinearOperator<Scalar> J = F.getJacobian();
00089
00090
00091 SolverState<Scalar> linSolverState = linSolver_.solve(J, resid, newtonStep);
00092 if (linSolverState.finalState() != SolveConverged)
00093 {
00094 PLAYA_MSG1(verb_, tab0 << " done Playa::NewtonArmijoSolver::solve()");
00095 return SolverState<Scalar>(SolveCrashed,
00096 "NewtonArmijoSolver::solve: linear solve failed with message ["
00097 + linSolverState.finalMsg() + "]", i, normF0);
00098 }
00099
00100
00101 Scalar t = ST::one();
00102
00103 bool stepAccepted = false;
00104 soln = F.currentEvalPt().copy();
00105
00106 for (int j=0; j<maxLineSearch_; j++)
00107 {
00108 Tabs tab2;
00109 Vector<Scalar> tmp = soln - t*newtonStep;
00110 F.setEvalPt( tmp );
00111 resid = F.getFunctionValue();
00112 ScalarMag normF1 = resid.norm2();
00113 PLAYA_MSG2(verb_, tab2 << "step t=" << setw(12) << t << " |F|=" << setw(12) << normF1);
00114 if (normF1 < (ST::one() - alpha_*t)*normF0)
00115 {
00116 stepAccepted = true;
00117 normF0 = normF1;
00118 break;
00119 }
00120 t = stepReduction_*t;
00121 }
00122
00123 if (!stepAccepted)
00124 {
00125 PLAYA_MSG1(verb_, tab0 << " done Playa::NewtonArmijoSolver::solve()");
00126 return SolverState<Scalar>(SolveCrashed,
00127 "NewtonArmijoSolver: line search failed",i, normF0);
00128 }
00129 }
00130
00131 PLAYA_MSG1(verb_, tab0 << " done Playa::NewtonArmijoSolver::solve()");
00132 return SolverState<Scalar>(SolveFailedToConverge, "NewtonArmijoSolver: convergence failure after "
00133 + Teuchos::toString(maxIters_) + " steps.", maxIters_, normF0);
00134 }
00135
00136 }
00137
00138
00139 #endif