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: `SNESLineSearch`, `SNESLineSearchType`, `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