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