1 #include <petsc/private/linesearchimpl.h> 2 #include <petscsnes.h> 3 4 static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch) { 5 PetscBool changed_y, changed_w; 6 Vec X, Y, F, W; 7 SNES snes; 8 PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep; 9 PetscReal lambda, lambda_old, lambda_update, delLambda; 10 PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s; 11 PetscInt i, max_its; 12 PetscViewer monitor; 13 14 PetscFunctionBegin; 15 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL)); 16 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 17 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 18 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 19 PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its)); 20 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 21 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 22 23 /* precheck */ 24 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 25 lambda_old = 0.0; 26 27 PetscCall(VecDot(F, Y, &fty_old)); 28 if (PetscAbsScalar(fty_old) < atol * ynorm) { 29 if (monitor) { 30 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 31 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm))); 32 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 33 } 34 PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS)); 35 PetscFunctionReturn(0); 36 } 37 38 fty_init = fty_old; 39 40 for (i = 0; i < max_its; i++) { 41 /* compute the norm at lambda */ 42 PetscCall(VecWAXPY(W, -lambda, Y, X)); 43 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 44 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 45 PetscCall(VecDot(F, Y, &fty)); 46 47 delLambda = lambda - lambda_old; 48 49 /* check for convergence */ 50 if (PetscAbsReal(delLambda) < steptol * lambda) break; 51 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break; 52 if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break; 53 if (monitor) { 54 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 55 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old))); 56 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 57 } 58 59 /* compute the search direction */ 60 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { 61 s = (fty - fty_old) / delLambda; 62 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 63 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 64 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 65 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 66 PetscCall(VecDot(F, Y, &fty_mid1)); 67 s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda; 68 } else { 69 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 70 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 71 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 72 PetscCall(VecDot(F, Y, &fty_mid1)); 73 PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X)); 74 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 75 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 76 PetscCall(VecDot(F, Y, &fty_mid2)); 77 s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda); 78 } 79 /* if the solve is going in the wrong direction, fix it */ 80 if (PetscRealPart(s) > 0.) s = -s; 81 if (s == 0.0) break; 82 lambda_update = lambda - PetscRealPart(fty / s); 83 84 /* switch directions if we stepped out of bounds */ 85 if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s); 86 87 if (PetscIsInfOrNanReal(lambda_update)) break; 88 if (lambda_update > maxstep) break; 89 90 /* compute the new state of the line search */ 91 lambda_old = lambda; 92 lambda = lambda_update; 93 fty_old = fty; 94 } 95 /* construct the solution */ 96 PetscCall(VecWAXPY(W, -lambda, Y, X)); 97 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 98 /* postcheck */ 99 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 100 if (changed_y) { 101 PetscCall(VecAXPY(X, -lambda, Y)); 102 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, X)); 103 } else { 104 PetscCall(VecCopy(W, X)); 105 } 106 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 107 108 PetscCall(SNESLineSearchComputeNorms(linesearch)); 109 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 110 111 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 112 113 if (monitor) { 114 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 115 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm)); 116 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 117 } 118 if (lambda <= steptol) { PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); } 119 PetscFunctionReturn(0); 120 } 121 122 /*MC 123 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 124 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 125 to find roots of dot(F, Y) via a secant method. 126 127 Options Database Keys: 128 + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda 129 . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value 130 . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0 131 - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed. 132 133 Notes: 134 This method does NOT use the objective function if it is provided with SNESSetObjective(). 135 136 This method is the preferred line search for SNESQN and SNESNCG. 137 138 Level: advanced 139 140 .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 141 M*/ 142 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) { 143 PetscFunctionBegin; 144 linesearch->ops->apply = SNESLineSearchApply_CP; 145 linesearch->ops->destroy = NULL; 146 linesearch->ops->setfromoptions = NULL; 147 linesearch->ops->reset = NULL; 148 linesearch->ops->view = NULL; 149 linesearch->ops->setup = NULL; 150 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 151 152 linesearch->max_its = 1; 153 PetscFunctionReturn(0); 154 } 155