MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
LineSearchSolver Class Reference

#include <linesearch.h>

Inheritance diagram for LineSearchSolver:
Collaboration diagram for LineSearchSolver:

Classes

class  MaxIterations
 
class  ResidualNotDecreased
 
class  SlopeNotNegative
 
class  ZeroGradient
 

Public Member Functions

 LineSearchSolver (DataManager *pDM, const NonlinearSolverOptions &options, const struct LineSearchParameters &param)
 
 ~LineSearchSolver (void)
 
void Solve (const NonlinearProblem *pNLP, Solver *pS, const integer iMaxIter, const doublereal &Tol, integer &iIterCnt, doublereal &dErr, const doublereal &SolTol, doublereal &dSolErr)
 
- Public Member Functions inherited from NonlinearSolver
 NonlinearSolver (const NonlinearSolverOptions &options)
 
virtual void SetTest (NonlinearSolverTest *pr, NonlinearSolverTest *ps)
 
virtual ~NonlinearSolver (void)
 
virtual bool MakeResTest (Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
 
virtual integer TotalAssembledJacobian (void)
 
virtual NonlinearSolverTestpGetResTest (void)
 
virtual NonlinearSolverTestpGetSolTest (void)
 
- Public Member Functions inherited from SolverDiagnostics
 SolverDiagnostics (unsigned OF=OUTPUT_DEFAULT, DriveCaller *pOM=0)
 
virtual ~SolverDiagnostics (void)
 
void SetNoOutput (void)
 
void SetOutputMeter (DriveCaller *pOM)
 
void SetOutputDriveHandler (const DriveHandler *pDH)
 
void SetOutputFlags (unsigned OF)
 
void AddOutputFlags (unsigned OF)
 
void DelOutputFlags (unsigned OF)
 
MatrixHandler::Norm_t GetCondMatNorm (void) const
 
bool outputMeter (void) const
 
bool outputIters (void) const
 
bool outputRes (void) const
 
bool outputSol (void) const
 
bool outputJac (void) const
 
bool outputStep (void) const
 
bool outputBailout (void) const
 
bool outputCounter (void) const
 
bool outputMatrixConditionNumber (void) const
 
bool outputSolverConditionNumber (void) const
 
bool outputSolverConditionStat (void) const
 
bool outputCPUTime (void) const
 
bool outputMsg (void) const
 

Private Member Functions

void LineSearch (doublereal stpmax, doublereal fold, doublereal &f, bool &check, const integer &iIterCnt)
 
void Residual (doublereal &f, integer iIterCnt)
 
void Jacobian ()
 
- Private Member Functions inherited from LineSearchParameters
 LineSearchParameters ()
 

Private Attributes

VectorHandlerpRes
 
VectorHandlerpSol
 
MyVectorHandler p
 
MyVectorHandler g
 
MyVectorHandler dXneg
 
const NonlinearProblempNLP
 
SolverpS
 
const DataManager *const pDM
 
- Private Attributes inherited from LineSearchParameters
doublereal dTolX
 
doublereal dTolMin
 
integer iMaxIterations
 
doublereal dMaxStep
 
doublereal dAlpha
 
doublereal dLambdaMin
 
doublereal dLambdaFactMin
 
doublereal dDivergenceCheck
 
doublereal dMinStepScale
 
unsigned uFlags
 

Additional Inherited Members

- Public Types inherited from NonlinearSolver
enum  Type {
  UNKNOWN = -1, NEWTONRAPHSON, MATRIXFREE, LINESEARCH,
  DEFAULT = NEWTONRAPHSON, LASTSOLVERTYPE
}
 
- Protected Types inherited from NonlinearSolver
enum  CPUTimeType { CPU_RESIDUAL, CPU_JACOBIAN, CPU_LINEAR_SOLVER, CPU_LAST_TYPE }
 
- Protected Types inherited from SolverDiagnostics
enum  {
  OUTPUT_NONE = 0x0000, OUTPUT_ITERS = 0x0001, OUTPUT_RES = 0x0002, OUTPUT_SOL = 0x0004,
  OUTPUT_JAC = 0x0008, OUTPUT_BAILOUT = 0x0010, OUTPUT_MSG = 0x0020, OUTPUT_COUNTER = 0x0040,
  OUTPUT_MAT_COND_NUM_1 = 0x0080, OUTPUT_MAT_COND_NUM_INF = 0x0100, OUTPUT_SOLVER_COND_NUM = 0x0200, OUTPUT_SOLVER_COND_STAT = 0x400,
  OUTPUT_CPU_TIME = 0x800, OUTPUT_MAT_COND_NUM = OUTPUT_MAT_COND_NUM_1 | OUTPUT_MAT_COND_NUM_INF, OUTPUT_DEFAULT = OUTPUT_MSG, OUTPUT_STEP = (OUTPUT_ITERS | OUTPUT_RES | OUTPUT_SOL | OUTPUT_JAC | OUTPUT_MAT_COND_NUM | OUTPUT_SOLVER_COND_NUM),
  OUTPUT_MASK = 0x07FF
}
 
- Protected Types inherited from NonlinearSolverOptions
enum  ScaleFlags { SCALE_ALGEBRAIC_EQUATIONS_NO = 0, SCALE_ALGEBRAIC_EQUATIONS_YES = 1 }
 
- Protected Member Functions inherited from NonlinearSolver
virtual bool MakeSolTest (Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
 
doublereal dGetCondMax () const
 
doublereal dGetCondMin () const
 
doublereal dGetCondAvg () const
 
doublereal dGetTimeCPU (CPUTimeType eType) const
 
void AddCond (doublereal dCond)
 
void AddTimeCPU (doublereal dTime, CPUTimeType eType)
 
- Protected Member Functions inherited from NonlinearSolverOptions
 NonlinearSolverOptions (bool bHonorJacRequest=false, enum ScaleFlags eScaleFlags=SCALE_ALGEBRAIC_EQUATIONS_NO, doublereal dScaleAlgebraic=1.)
 
- Protected Attributes inherited from NonlinearSolver
integer Size
 
integer TotJac
 
NonlinearSolverTestpResTest
 
NonlinearSolverTestpSolTest
 
- Protected Attributes inherited from SolverDiagnostics
unsigned OutputFlags
 
DriveCallerpOutputMeter
 
- Protected Attributes inherited from NonlinearSolverOptions
bool bHonorJacRequest
 
enum
NonlinearSolverOptions::ScaleFlags 
eScaleFlags
 
doublereal dScaleAlgebraic
 
- Private Types inherited from LineSearchParameters
enum  {
  ZERO_GRADIENT_CONTINUE = 0x001, DIVERGENCE_CHECK = 0x002, ALGORITHM_CUBIC = 0x004, ALGORITHM_FACTOR = 0x008,
  PRINT_CONVERGENCE_INFO = 0x010, SCALE_NEWTON_STEP = 0x020, RELATIVE_LAMBDA_MIN = 0x040, ABORT_AT_LAMBDA_MIN = 0x080,
  VERBOSE_MODE = 0x100, NON_NEGATIVE_SLOPE_CONTINUE = 0x200, ALGORITHM = ALGORITHM_CUBIC | ALGORITHM_FACTOR
}
 

Detailed Description

Definition at line 81 of file linesearch.h.

Constructor & Destructor Documentation

LineSearchSolver::LineSearchSolver ( DataManager pDM,
const NonlinearSolverOptions options,
const struct LineSearchParameters param 
)

Definition at line 102 of file linesearch.cc.

References LineSearchParameters::ALGORITHM_CUBIC, LineSearchParameters::ALGORITHM_FACTOR, NonlinearSolverOptions::bHonorJacRequest, LineSearchParameters::dAlpha, LineSearchParameters::dDivergenceCheck, LineSearchParameters::DIVERGENCE_CHECK, LineSearchParameters::dLambdaFactMin, LineSearchParameters::dLambdaMin, LineSearchParameters::dMaxStep, LineSearchParameters::dMinStepScale, LineSearchParameters::dTolMin, LineSearchParameters::dTolX, LineSearchParameters::iMaxIterations, LineSearchParameters::PRINT_CONVERGENCE_INFO, LineSearchParameters::RELATIVE_LAMBDA_MIN, LineSearchParameters::SCALE_NEWTON_STEP, TRACE_FLAG, TRACE_VAR, LineSearchParameters::uFlags, and LineSearchParameters::ZERO_GRADIENT_CONTINUE.

105 : NonlinearSolver(options),
106 LineSearchParameters(param),
107 pRes(0),
108 pSol(0),
109 pNLP(0),
110 pS(0),
111 pDM(pDM)
112 {
113  TRACE_VAR(dTolX);
114  // TRACE_VAR(dTolF);
118  TRACE_VAR(dAlpha);
124  TRACE_VAR(uFlags);
132 }
VectorHandler * pRes
Definition: linesearch.h:83
doublereal dMaxStep
Definition: linesearch.h:60
doublereal dTolMin
Definition: linesearch.h:58
VectorHandler * pSol
Definition: linesearch.h:84
const NonlinearProblem * pNLP
Definition: linesearch.h:88
const DataManager *const pDM
Definition: linesearch.h:90
doublereal dLambdaMin
Definition: linesearch.h:62
doublereal dTolX
Definition: linesearch.h:57
doublereal dDivergenceCheck
Definition: linesearch.h:64
#define TRACE_FLAG(varname, flag)
Definition: linesearch.cc:81
NonlinearSolver(const NonlinearSolverOptions &options)
Definition: nonlin.cc:434
#define TRACE_VAR(varname)
Definition: linesearch.cc:80
doublereal dMinStepScale
Definition: linesearch.h:65
doublereal dAlpha
Definition: linesearch.h:61
integer iMaxIterations
Definition: linesearch.h:59
doublereal dLambdaFactMin
Definition: linesearch.h:63
LineSearchSolver::~LineSearchSolver ( void  )

Definition at line 134 of file linesearch.cc.

References NO_OP.

135 {
136  NO_OP;
137 }
#define NO_OP
Definition: myassert.h:74

Member Function Documentation

void LineSearchSolver::Jacobian ( )
private

Definition at line 163 of file linesearch.cc.

References MatrixHandler::ConditionNumber(), SolverDiagnostics::GetCondMatNorm(), NonlinearProblem::Jacobian(), SolutionManager::MatrInitialize(), SolutionManager::MatrReset(), MBDYN_EXCEPT_ARGS, SolverDiagnostics::outputJac(), SolverDiagnostics::outputMatrixConditionNumber(), Solver::pGetSolutionManager(), SolutionManager::pMatHdl(), pNLP, pS, and NonlinearSolver::TotJac.

Referenced by Solve().

164 {
165  SolutionManager *const pSM = pS->pGetSolutionManager();
166 
167  const integer iMaxIterRebuild = 10;
168 
169  pSM->MatrReset();
170  integer iIter = 0;
171 
172  do {
173  try {
174  pNLP->Jacobian(pSM->pMatHdl());
175  } catch (const MatrixHandler::ErrRebuildMatrix&) {
176  silent_cout("NewtonRaphsonSolver: "
177  "rebuilding matrix..."
178  << std::endl);
179 
180  /* need to rebuild the matrix... */
181  pSM->MatrInitialize();
182  continue;
183  }
184  break;
185  } while (++iIter < iMaxIterRebuild);
186 
187  if (iIter >= iMaxIterRebuild) {
188  silent_cerr("Maximum number of iterations exceeded when rebuilding the Jacobian matrix" << std::endl);
190  }
191 
192  TotJac++;
193 
194 #ifdef USE_MPI
195  if (!bParallel || MBDynComm.Get_rank() == 0)
196 #endif /* USE_MPI */
197  {
198  if (outputJac()) {
199  silent_cout("Jacobian:" << std::endl
200  << *(pSM->pMatHdl()));
201  }
202 
204  silent_cout(" cond=" << pSM->pMatHdl()->ConditionNumber(GetCondMatNorm()) << std::endl);
205  }
206  }
207 }
bool outputMatrixConditionNumber(void) const
virtual SolutionManager * pGetSolutionManager(void) const
Definition: solver.h:398
virtual void Jacobian(MatrixHandler *pJac) const =0
integer TotJac
Definition: nonlin.h:252
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
MatrixHandler::Norm_t GetCondMatNorm(void) const
bool outputJac(void) const
const NonlinearProblem * pNLP
Definition: linesearch.h:88
virtual doublereal ConditionNumber(enum Norm_t eNorm=NORM_1) const
Definition: mh.cc:451
virtual MatrixHandler * pMatHdl(void) const =0
virtual void MatrReset(void)=0
virtual void MatrInitialize(void)
Definition: solman.cc:72
long int integer
Definition: colamd.c:51

Here is the call graph for this function:

void LineSearchSolver::LineSearch ( doublereal  stpmax,
doublereal  fold,
doublereal f,
bool check,
const integer iIterCnt 
)
private

Definition at line 210 of file linesearch.cc.

References a, LineSearchParameters::ALGORITHM_CUBIC, ASSERT, Solver::CheckTimeStepLimit(), LineSearchParameters::dAlpha, DataManager::dGetTime(), LineSearchParameters::dLambdaFactMin, LineSearchParameters::dLambdaMin, LineSearchParameters::dMinStepScale, LineSearchParameters::dTolX, dXneg, g, DataManager::GetpXCurr(), DataManager::GetpXPCurr(), LineSearchParameters::iMaxIterations, VectorHandler::InnerProd(), NonlinearSolver::MakeResTest(), MBDYN_EXCEPT_ARGS, mbdyn_stop_at_end_of_iteration(), LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE, VectorHandler::Norm(), SolverDiagnostics::outputIters(), p, pDM, pNLP, pRes, LineSearchParameters::PRINT_CONVERGENCE_INFO, pS, pSol, LineSearchParameters::RELATIVE_LAMBDA_MIN, VectorHandler::Reset(), Residual(), VectorHandler::ScalarMul(), LineSearchParameters::SCALE_NEWTON_STEP, NonlinearSolver::Size, grad::sqrt(), TRACE, TRACE_VAR, LineSearchParameters::uFlags, NonlinearProblem::Update(), and LineSearchParameters::VERBOSE_MODE.

Referenced by Solve().

212 {
213  check = false;
214 
215  if (uFlags & SCALE_NEWTON_STEP) {
216  const doublereal dNormSol = pSol->Norm();
217 
218  if (dNormSol > stpmax) {
219  const doublereal dScale = std::max(dMinStepScale, stpmax / dNormSol);
220  if (uFlags & VERBOSE_MODE) {
221  silent_cerr("line search warning: "
222  "Newton increment is reduced by factor " << dScale
223  << " at time step " << pDM->dGetTime()
224  << " at iteration " << iIterCnt
225  << " The time step is probably too large" << std::endl);
226  }
227  ASSERT(stpmax >= 0);
228  ASSERT(dScale <= 1.);
229  ASSERT(dScale > 0.);
230  *pSol *= dScale;
231  }
232  }
233 
234  p = *pSol; // save the Newton increment
235 
236  doublereal slope = g.InnerProd(p);
237 
238  TRACE_VAR(slope);
239 
240  doublereal dLambdaMinEff = -1.;
241 
242  if (slope >= 0.) {
243  if (uFlags & VERBOSE_MODE) {
244  silent_cerr("line search warning: slope=" << slope << " >= 0"
245  " at time step " << pDM->dGetTime()
246  << " at iteration " << iIterCnt << std::endl
247  << "This could be a roundoff problem" << std::endl);
248  }
249 
251  // It seems to be a numerical problem.
252  // Line search may not work in this situation.
253  // Resort to the ordinary Newton Raphson algorithm
254  slope = 0.;
255  dLambdaMinEff = 1.;
256  } else {
257  throw SlopeNotNegative(MBDYN_EXCEPT_ARGS);
258  }
259  }
260 
261  if (dLambdaMinEff < 0) {
262  // dLambdaMinEff has to be detected
263  if (uFlags & RELATIVE_LAMBDA_MIN) {
264  doublereal test = 0.;
265 
266  for (integer i = 1; i <= Size; ++i) {
267  const doublereal temp = std::abs(p(i)) / std::max(std::max(std::abs((*pDM->GetpXCurr())(i)), std::abs((*pDM->GetpXPCurr())(i))), 1.);
268  if (temp > test) {
269  test = temp;
270  }
271  }
272 
273  dLambdaMinEff = std::max(dLambdaMin, std::min(dTolX / test, 1.));
274  } else {
275  dLambdaMinEff = dLambdaMin;
276  }
277  }
278 
279  doublereal lambda = 1.;
280  doublereal tmplam;
281  doublereal lambda2 = 0; //f2, lambda2 initialized to silence g++ possible uninitialized warning
282  doublereal f2 = 0.;
283 
284 #ifdef DEBUG
285  lambda2 = f2 = NAN;
286 #endif
287 
288  integer iIter = 0;
289 
290  TRACE_VAR(dLambdaMinEff);
291 
292  do {
293  // FIXME: dLambdaMax not defined
294  // ASSERT(lambda <= dLambdaMax);
295 
296  if (iIter > 0) {
297  TRACE("Start new step from Xold, XPold with lambda = " << lambda << " ..." << std::endl);
298  pNLP->Update(&dXneg); // restore the previous state
299 
300  pSol->Reset();
301  pSol->ScalarMul(p, lambda); // scale the Newton increment by lambda
302  }
303 
304  for (integer i = 1; i <= Size; ++i)
305  dXneg(i) = -(*pSol)(i); // save the increment; so we can restore the previous state
306 
307  TRACE("Update the nonlinear problem ... pSol->Norm()=" << pSol->Norm() << std::endl);
308 
309  pNLP->Update(pSol);
310  Residual(f, iIter);
311 
312  ++iIter;
313 
314  TRACE("New value for f:" << f << std::endl);
315 
316  doublereal dErr = 0., dErrDiff = 0.;
317  MakeResTest(pS, pNLP, *pRes, 0., dErr, dErrDiff);
318 
320 #ifdef USE_MPI
321  if (!bParallel || MBDynComm.Get_rank() == 0)
322 #endif /* USE_MPI */
323  {
324  silent_cout("\t\tf(" << iIterCnt << ":" << iIter << ")=" << std::setw(12) << f
325  << "\tErr=" << std::setw(12) << dErr
326  << "\tlambda=" << std::setw(12) << lambda
327  <<"\tslope=" << slope << std::endl);
328  }
329  }
330 
331  pS->CheckTimeStepLimit(dErr, dErrDiff);
332 
333  if (f <= fold + dAlpha * lambda * slope) {
334  TRACE("Sufficient decrease in f: backtrack" << std::endl);
335  return;
336  } else if (lambda <= dLambdaMinEff) {
337  TRACE("Checking for convergence: lambda=" << lambda << " < lambdaMin=" << dLambdaMinEff << std::endl);
338  check = true; // check for convergence
339  return;
340  } else {
341  if (uFlags & ALGORITHM_CUBIC) {
342  TRACE("Calculate new value for alam ..." << std::endl);
343  if (lambda == 1.) {
344  tmplam = -slope / (2 * (f - fold - slope));
345  TRACE_VAR(tmplam);
346  } else {
347  const doublereal rhs1 = f - fold - lambda * slope;
348  const doublereal rhs2 = f2 - fold - lambda2 * slope;
349  const doublereal a = (rhs1 / (lambda * lambda) - rhs2 / (lambda2 * lambda2)) / (lambda - lambda2);
350  const doublereal b = (-lambda2 * rhs1 / (lambda * lambda) + lambda * rhs2 / (lambda2 * lambda2)) / (lambda - lambda2);
351 
352  if (a == 0.) {
353  tmplam = -slope / (2. * b);
354  TRACE_VAR(tmplam);
355  } else {
356  const doublereal disc = b * b - 3. * a * slope;
357 
358  if (disc < 0.) {
359  tmplam = 0.5 * lambda;
360  TRACE_VAR(tmplam);
361  } else if (b <= 0.) {
362  tmplam = (-b + sqrt(disc)) / (3. * a);
363  TRACE_VAR(tmplam);
364  } else {
365  tmplam = -slope / (b + sqrt(disc));
366  TRACE_VAR(tmplam);
367  }
368 
369  if (tmplam > 0.5 * lambda) {
370  tmplam = 0.5 * lambda;
371  TRACE_VAR(tmplam);
372  }
373  }
374  }
375  lambda2 = lambda;
376  f2 = f;
377  TRACE_VAR(tmplam);
378  lambda = std::max(tmplam, dLambdaFactMin * lambda);
379  } else {
380  lambda *= dLambdaFactMin;
381  }
382  }
383 
386  }
387  } while (iIter < iMaxIterations);
388 
389  if (uFlags & VERBOSE_MODE) {
390  silent_cerr("line search warning: Maximum number of line search iterations="
391  << iMaxIterations
392  << " exceeded in line search at time step " << pDM->dGetTime()
393  << " at iteration " << iIterCnt << std::endl);
394  }
395 
396  throw MaxIterations(MBDYN_EXCEPT_ARGS);
397 }
virtual void Reset(void)=0
bool outputIters(void) const
MyVectorHandler p
Definition: linesearch.h:85
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
VectorHandler * pRes
Definition: linesearch.h:83
#define TRACE(expr)
Definition: linesearch.cc:77
virtual bool MakeResTest(Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
Definition: nonlin.cc:474
doublereal dGetTime(void) const
Definition: dataman2.cc:165
virtual void CheckTimeStepLimit(doublereal dErr, doublereal dErrDiff) const
Definition: solver.cc:5205
virtual doublereal InnerProd(const VectorHandler &VH) const
Definition: vh.cc:269
MyVectorHandler g
Definition: linesearch.h:86
integer Size
Definition: nonlin.h:251
VectorHandler * pSol
Definition: linesearch.h:84
virtual doublereal Norm(void) const
Definition: vh.cc:262
const NonlinearProblem * pNLP
Definition: linesearch.h:88
const DataManager *const pDM
Definition: linesearch.h:90
int mbdyn_stop_at_end_of_iteration(void)
Definition: solver.cc:172
doublereal dLambdaMin
Definition: linesearch.h:62
doublereal dTolX
Definition: linesearch.h:57
void Residual(doublereal &f, integer iIterCnt)
Definition: linesearch.cc:140
const VectorHandler * GetpXPCurr(void) const
Definition: dataman.h:857
virtual VectorHandler & ScalarMul(const VectorHandler &VH, const doublereal &d)
Definition: vh.cc:145
MyVectorHandler dXneg
Definition: linesearch.h:87
#define ASSERT(expression)
Definition: colamd.c:977
const VectorHandler * GetpXCurr(void) const
Definition: dataman.h:853
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
#define TRACE_VAR(varname)
Definition: linesearch.cc:80
doublereal dMinStepScale
Definition: linesearch.h:65
doublereal dAlpha
Definition: linesearch.h:61
integer iMaxIterations
Definition: linesearch.h:59
static const doublereal a
Definition: hfluid_.h:289
doublereal dLambdaFactMin
Definition: linesearch.h:63
virtual void Update(const VectorHandler *pSol) const =0
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51

Here is the call graph for this function:

void LineSearchSolver::Residual ( doublereal f,
integer  iIterCnt 
)
private

Definition at line 140 of file linesearch.cc.

References VectorHandler::Dot(), SolverDiagnostics::outputRes(), pNLP, pRes, Solver::PrintResidual(), pS, VectorHandler::Reset(), and NonlinearProblem::Residual().

Referenced by LineSearch(), and Solve().

141 {
142 #ifdef USE_EXTERNAL
143  SendExternal();
144 #endif /* USE_EXTERNAL */
145 
146  pRes->Reset();
147 
148  try {
149  pNLP->Residual(pRes);
150  }
152  // The Jacobian will be rebuilt anyway
153  }
154 
155  if (outputRes()) {
156  pS->PrintResidual(*pRes, iIterCnt);
157  }
158 
159  f = 0.5 * pRes->Dot();
160 }
virtual void Reset(void)=0
VectorHandler * pRes
Definition: linesearch.h:83
virtual void PrintResidual(const VectorHandler &Res, integer iIterCnt) const
Definition: solver.cc:5194
bool outputRes(void) const
virtual void Residual(VectorHandler *pRes) const =0
virtual doublereal Dot(void) const
Definition: vh.cc:244
const NonlinearProblem * pNLP
Definition: linesearch.h:88

Here is the call graph for this function:

void LineSearchSolver::Solve ( const NonlinearProblem pNLP,
Solver pS,
const integer  iMaxIter,
const doublereal Tol,
integer iIterCnt,
doublereal dErr,
const doublereal SolTol,
doublereal dSolErr 
)
virtual

Implements NonlinearSolver.

Definition at line 400 of file linesearch.cc.

References LineSearchParameters::ABORT_AT_LAMBDA_MIN, NonlinearSolver::AddCond(), NonlinearSolver::AddTimeCPU(), ASSERT, SolutionManager::bGetConditionNumber(), Solver::CheckTimeStepLimit(), NonlinearSolver::CPU_JACOBIAN, NonlinearSolver::CPU_LINEAR_SOLVER, NonlinearSolver::CPU_RESIDUAL, LineSearchParameters::dDivergenceCheck, NonlinearSolver::dGetCondAvg(), NonlinearSolver::dGetCondMax(), NonlinearSolver::dGetCondMin(), DataManager::dGetTime(), NonlinearSolver::dGetTimeCPU(), LineSearchParameters::DIVERGENCE_CHECK, LineSearchParameters::dMaxStep, VectorHandler::Dot(), LineSearchParameters::dTolMin, dXneg, g, DataManager::GetpXCurr(), DataManager::GetpXPCurr(), MatrixHandler::iGetNumCols(), MatrixHandler::iGetNumRows(), VectorHandler::iGetSize(), MyVectorHandler::iGetSize(), Jacobian(), LineSearch(), NonlinearSolver::MakeResTest(), NonlinearSolver::MakeSolTest(), MatrixHandler::MatTVecDecMul(), mbdyn_clock_time(), MBDYN_EXCEPT_ARGS, SolverDiagnostics::outputBailout(), SolverDiagnostics::outputCPUTime(), SolverDiagnostics::outputIters(), SolverDiagnostics::outputSol(), SolverDiagnostics::outputSolverConditionNumber(), SolverDiagnostics::outputSolverConditionStat(), p, pDM, Solver::pGetSolutionManager(), SolutionManager::pMatHdl(), pNLP, pRes, SolutionManager::pResHdl(), LineSearchParameters::PRINT_CONVERGENCE_INFO, Solver::PrintResidual(), Solver::PrintSolution(), pS, pSol, SolutionManager::pSolHdl(), MyVectorHandler::Reset(), Residual(), MyVectorHandler::Resize(), LineSearchParameters::SCALE_NEWTON_STEP, NonlinearSolver::Size, SolutionManager::Solve(), grad::sqrt(), TRACE, TRACE_VAR, LineSearchParameters::uFlags, LineSearchParameters::VERBOSE_MODE, and LineSearchParameters::ZERO_GRADIENT_CONTINUE.

408 {
409  ASSERT(pS != NULL);
410  SolutionManager *const pSM = pS->pGetSolutionManager();
411  this->pS = pS;
412  pNLP = pNonLinProblem;
413  pRes = pSM->pResHdl();
414  pSol = pSM->pSolHdl();
415  Size = pRes->iGetSize();
416 
417  if (g.iGetSize() != Size) {
418  TRACE("Resize temporary vectors ..." << std::endl);
419  g.Resize(Size);
420  p.Resize(Size);
421  dXneg.Resize(Size);
422  }
423 
424  bool check = false;
425  iIterCnt = 0;
426  dSolErr = 0.;
427  dErr = 0.;
428  doublereal dErrDiff = 0.;
429  doublereal f;
430  Residual(f, iIterCnt);
431 
432  TRACE("\t\tf(0) = " << f << std::endl);
433 
434  bool bTest = MakeResTest(pS, pNLP, *pRes, 1e-2 * Tol, dErr, dErrDiff); // use a more stringent test for the first iteration
435  doublereal dPrevErr = std::numeric_limits<doublereal>::max(); // disable error test for the first iteration
436  doublereal dErrFactor = 1.;
437 
438  if (outputIters()) {
439 #ifdef USE_MPI
440  if (!bParallel || MBDynComm.Get_rank() == 0)
441 #endif /* USE_MPI */
442  {
443  silent_cout("\tIteration(" << iIterCnt << ") " << std::setw(12) << dErr);
444 
446  silent_cout(" f=" << std::setw(12) << f);
447  }
448 
449  silent_cout(std::endl);
450  }
451  }
452 
453  pS->CheckTimeStepLimit(dErr, dErrDiff);
454 
455  if (bTest) {
456  return;
457  }
458 
459  const doublereal stpmax = (uFlags & SCALE_NEWTON_STEP)
460  ? dMaxStep * std::max(sqrt(pDM->GetpXPCurr()->Dot() + pDM->GetpXCurr()->Dot()), static_cast<doublereal>(Size))
461  : std::numeric_limits<doublereal>::max();
462 
463  TRACE_VAR(stpmax);
464 
465  doublereal dStartTimeCPU, dEndTimeCPU, dJacobianCPU, dResidualCPU, dLinSolveCPU;
466 
467  while (true) {
468  if (outputCPUTime()) {
469  dStartTimeCPU = mbdyn_clock_time();
470  }
471 
472  Jacobian();
473 
474  ASSERT(pSM->pMatHdl()->iGetNumCols() == Size);
475  ASSERT(pSM->pMatHdl()->iGetNumCols() == pRes->iGetSize());
476  ASSERT(pSM->pMatHdl()->iGetNumRows() == pRes->iGetSize());
477 
478  g.Reset();
479  pSM->pMatHdl()->MatTVecDecMul(g, *pRes); // compute gradient g = \nabla f = fjac^T \, fvec = -Jac^T \, pRes
480 
481  if (outputCPUTime()) {
482  dEndTimeCPU = mbdyn_clock_time();
483  dJacobianCPU = dEndTimeCPU - dStartTimeCPU;
484  AddTimeCPU(dJacobianCPU, CPU_JACOBIAN);
485  dStartTimeCPU = dEndTimeCPU;
486  }
487 
488  const doublereal fold = f;
489 
490  pSM->Solve();
491 
492  if (outputCPUTime()) {
493  dEndTimeCPU = mbdyn_clock_time();
494  dLinSolveCPU = dEndTimeCPU - dStartTimeCPU;
495  AddTimeCPU(dLinSolveCPU, CPU_LINEAR_SOLVER);
496  }
497 
498  if (outputSol()) {
499  pS->PrintSolution(*pSol, iIterCnt);
500  }
501 
502  if (outputCPUTime()) {
503  dStartTimeCPU = mbdyn_clock_time();
504  }
505 
506  LineSearch(stpmax, fold, f, check, iIterCnt);
507 
508  bTest = MakeResTest(pS, pNLP, *pRes, Tol, dErr, dErrDiff);
509 
510  if (outputCPUTime()) {
511  dEndTimeCPU = mbdyn_clock_time();
512  dResidualCPU = dEndTimeCPU - dStartTimeCPU;
513  AddTimeCPU(dResidualCPU, CPU_RESIDUAL);
514  }
515 
516  if (iIterCnt > 0) {
517  dErrFactor *= dErr / dPrevErr;
518  }
519 
520  iIterCnt++;
521 
522 #ifdef USE_MPI
523  if (!bParallel || MBDynComm.Get_rank() == 0)
524 #endif /* USE_MPI */
525  {
527  if (outputIters()) {
528  silent_cout("\tIteration(" << iIterCnt << ") " << std::setw(12) << dErr << " J");
529  }
530 
532  silent_cout(" cond=");
533  doublereal dCond;
534  if (pSM->bGetConditionNumber(dCond)) {
535  silent_cout(dCond);
537  AddCond(dCond);
538  silent_cout(" " << dGetCondMin() << " " << dGetCondMax() << " " << dGetCondAvg());
539  }
540  } else {
541  silent_cout("NA");
542  }
543  }
544 
545  if (outputCPUTime()) {
546  silent_cout(" CPU:" << dResidualCPU << "/" << dGetTimeCPU(CPU_RESIDUAL)
547  << "+" << dJacobianCPU << "/" << dGetTimeCPU(CPU_JACOBIAN)
548  << "+" << dLinSolveCPU << "/" << dGetTimeCPU(CPU_LINEAR_SOLVER));
549  }
550 
552  silent_cout(" f=" << std::setw(12) << f);
553 
554  if (iIterCnt > 1) {
555  silent_cout(" Err(n-1)=" << std::setw(12) << dPrevErr
556  << " Err(n)/Err(n-1)=" << std::setw(12) << dErr / dPrevErr
557  << " Err(n)/Err(1)=" << std::setw(12) << dErrFactor);
558  }
559  }
560 
561  if (outputIters()
563  || (uFlags & PRINT_CONVERGENCE_INFO)) {
564  silent_cout(std::endl);
565  }
566  }
567  }
568 
569  if (bTest) {
570  return;
571  }
572 
573  bTest = MakeSolTest(pS, *pSol, SolTol, dSolErr);
574 
575  if (outputIters()) {
576 #ifdef USE_MPI
577  if (!bParallel || MBDynComm.Get_rank() == 0)
578 #endif /* USE_MPI */
579  {
580  silent_cout("\t\tSolErr "
581  << dSolErr << std::endl);
582  }
583  }
584 
585  if (bTest) {
586  if (uFlags & VERBOSE_MODE) {
587  silent_cerr("line search warning: Convergence on solution "
588  "at time step " << pDM->dGetTime()
589  << "at iteration " << iIterCnt
590  << "\tErr(n)=" << dErr);
591 
592  if (iIterCnt >= 1) {
593  silent_cerr("\tErr(n-1)=" << dPrevErr
594  << "\tErr(n)/Err(n-1)=" << dErr/dPrevErr
595  << "\tErr(n)/Err(1) " << dErrFactor);
596  }
597 
598  silent_cerr(std::endl);
599  }
600  throw ConvergenceOnSolution(MBDYN_EXCEPT_ARGS);
601  }
602 
603  if (check) { // lambda <= dLambdaMinEff: check for gradient zero
604  doublereal test = 0.;
605  const doublereal den = std::max(f, 0.5 * Size);
606 
607  for (integer i = 1; i <= Size; ++i) {
608  const doublereal absX = std::max(std::abs((*pDM->GetpXCurr())(i)),
609  std::abs((*pDM->GetpXPCurr())(i)));
610  const doublereal temp = std::abs(g(i)) * std::max(absX, 1.) / den;
611  if (temp > test) {
612  test = temp;
613  }
614  }
615 
616  if (test < dTolMin) { // we are at a local minimum of the function f
617  if (uFlags & VERBOSE_MODE) {
618  silent_cerr("line search warning: Zero gradient detected at time step "
619  << pDM->dGetTime() << " at iteration "
620  << iIterCnt << " (spurious convergence) test=" << test
621  << " < tol=" << dTolMin
622  << "\tErr(n)=" << dErr
623  << " > Tol = " << Tol << std::endl);
624  }
625 
627  return;
628  } else {
629  throw ZeroGradient(MBDYN_EXCEPT_ARGS);
630  }
631  } else {
632  if (uFlags & VERBOSE_MODE) {
633  silent_cerr("line search warning: lambda min"
634  " has been reached at time step " << pDM->dGetTime()
635  << " at iteration " << iIterCnt
636  << " but the gradient is not zero" << std::endl);
637  }
638 
639  if (uFlags & ABORT_AT_LAMBDA_MIN) {
640  throw NoConvergence(MBDYN_EXCEPT_ARGS);
641  }
642  }
643  } else if (dErr > dPrevErr) { // f <= fold + dAlpha * lambda * slope
644  // We should not get here since f was decreased also Err should be decreased
645  // If not it could be a numerical problem (e.g. the tolerance could be too small)
646  if (uFlags & VERBOSE_MODE) {
647  silent_cerr("line search warning: f has been reduced during line search but the residual could not be reduced at time step "
648  << pDM->dGetTime()
649  << " at iteration " << iIterCnt << std::endl);
650  }
651 
652  // Do not throw an exception here because if we specify for example
653  // tolerance: <<Tol>>, test, minmax;
654  // it could happen that the the norm of the residual vector is decreased
655  // but the maximum residual is not!
656 
657  // In any case we will check for divergence if DIVERGENCE_CHECK is enabled.
658  }
659 
660  if (dErrFactor > dDivergenceCheck) {
661  if (uFlags & DIVERGENCE_CHECK) {
662  if (uFlags & VERBOSE_MODE) {
663  silent_cerr("line search warning: The residual could not be decreased"
664  " at time step " << pDM->dGetTime()
665  << " at iteration " << iIterCnt << std::endl);
666 
667  if (iIterCnt > 1) {
668  silent_cerr("Err(n-1)=" << dPrevErr
669  << "\tErr(n)=" << dErr
670  << "\tErr(n)/Err(n-1)=" << dErr/dPrevErr
671  << "\tErr(n)/Err(1)=" << dErrFactor << std::endl);
672  }
673  }
674  throw ResidualNotDecreased(MBDYN_EXCEPT_ARGS);
675  }
676  }
677 
678  if (iIterCnt >= std::abs(iMaxIter)) {
679  if (iMaxIter < 0 && dErrFactor < 1.) {
680  return;
681  }
682  if (outputBailout()) {
683  pS->PrintResidual(*pRes, iIterCnt);
684  }
685  throw NoConvergence(MBDYN_EXCEPT_ARGS);
686  }
687 
688  dPrevErr = dErr;
689  }
690 }
virtual VectorHandler * pResHdl(void) const =0
virtual SolutionManager * pGetSolutionManager(void) const
Definition: solver.h:398
void AddCond(doublereal dCond)
Definition: nonlin.h:331
doublereal dGetCondMax() const
Definition: nonlin.h:275
virtual integer iGetNumCols(void) const =0
bool outputIters(void) const
MyVectorHandler p
Definition: linesearch.h:85
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
VectorHandler * pRes
Definition: linesearch.h:83
bool outputBailout(void) const
#define TRACE(expr)
Definition: linesearch.cc:77
doublereal dMaxStep
Definition: linesearch.h:60
virtual bool MakeSolTest(Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
Definition: nonlin.cc:497
virtual void PrintResidual(const VectorHandler &Res, integer iIterCnt) const
Definition: solver.cc:5194
virtual bool MakeResTest(Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
Definition: nonlin.cc:474
doublereal dGetTime(void) const
Definition: dataman2.cc:165
double mbdyn_clock_time()
Definition: clock_time.cc:51
virtual VectorHandler & MatTVecDecMul(VectorHandler &out, const VectorHandler &in) const
Definition: mh.cc:368
virtual void CheckTimeStepLimit(doublereal dErr, doublereal dErrDiff) const
Definition: solver.cc:5205
MyVectorHandler g
Definition: linesearch.h:86
doublereal dGetCondAvg() const
Definition: nonlin.h:277
doublereal dTolMin
Definition: linesearch.h:58
void LineSearch(doublereal stpmax, doublereal fold, doublereal &f, bool &check, const integer &iIterCnt)
Definition: linesearch.cc:210
integer Size
Definition: nonlin.h:251
VectorHandler * pSol
Definition: linesearch.h:84
doublereal dGetCondMin() const
Definition: nonlin.h:276
virtual doublereal Dot(void) const
Definition: vh.cc:244
virtual integer iGetSize(void) const =0
doublereal dGetTimeCPU(CPUTimeType eType) const
Definition: nonlin.h:346
bool outputSolverConditionNumber(void) const
virtual void Resize(integer iSize)
Definition: vh.cc:347
const DataManager *const pDM
Definition: linesearch.h:90
virtual integer iGetSize(void) const
Definition: vh.h:255
bool outputSol(void) const
doublereal dDivergenceCheck
Definition: linesearch.h:64
void Residual(doublereal &f, integer iIterCnt)
Definition: linesearch.cc:140
virtual MatrixHandler * pMatHdl(void) const =0
const VectorHandler * GetpXPCurr(void) const
Definition: dataman.h:857
void AddTimeCPU(doublereal dTime, CPUTimeType eType)
Definition: nonlin.h:355
virtual void Solve(void)=0
MyVectorHandler dXneg
Definition: linesearch.h:87
#define ASSERT(expression)
Definition: colamd.c:977
const VectorHandler * GetpXCurr(void) const
Definition: dataman.h:853
GradientExpression< UnaryExpr< FuncSqrt, Expr > > sqrt(const GradientExpression< Expr > &u)
Definition: gradient.h:2974
virtual void Reset(void)
Definition: vh.cc:459
#define TRACE_VAR(varname)
Definition: linesearch.cc:80
bool bGetConditionNumber(doublereal &dCond) const
Definition: solman.cc:107
virtual void PrintSolution(const VectorHandler &Sol, integer iIterCnt) const
Definition: solver.cc:5200
bool outputCPUTime(void) const
virtual VectorHandler * pSolHdl(void) const =0
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
bool outputSolverConditionStat(void) const
virtual integer iGetNumRows(void) const =0

Here is the call graph for this function:

Member Data Documentation

MyVectorHandler LineSearchSolver::dXneg
private

Definition at line 87 of file linesearch.h.

Referenced by LineSearch(), and Solve().

MyVectorHandler LineSearchSolver::g
private

Definition at line 86 of file linesearch.h.

Referenced by LineSearch(), and Solve().

MyVectorHandler LineSearchSolver::p
private

Definition at line 85 of file linesearch.h.

Referenced by LineSearch(), and Solve().

const DataManager* const LineSearchSolver::pDM
private

Definition at line 90 of file linesearch.h.

Referenced by LineSearch(), and Solve().

const NonlinearProblem* LineSearchSolver::pNLP
private

Definition at line 88 of file linesearch.h.

Referenced by Jacobian(), LineSearch(), Residual(), and Solve().

VectorHandler* LineSearchSolver::pRes
private

Definition at line 83 of file linesearch.h.

Referenced by LineSearch(), Residual(), and Solve().

Solver* LineSearchSolver::pS
private

Definition at line 89 of file linesearch.h.

Referenced by Jacobian(), LineSearch(), Residual(), and Solve().

VectorHandler* LineSearchSolver::pSol
private

Definition at line 84 of file linesearch.h.

Referenced by LineSearch(), and Solve().


The documentation for this class was generated from the following files: