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 if (linesearch->ops->vidirderiv) { 29 PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old)); 30 } else { 31 PetscCall(VecDot(F, Y, &fty_old)); 32 } 33 if (PetscAbsScalar(fty_old) < atol * ynorm) { 34 if (monitor) { 35 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 36 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))); 37 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 38 } 39 PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS)); 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 fty_init = fty_old; 44 45 for (i = 0; i < max_its; i++) { 46 /* compute the norm at lambda */ 47 PetscCall(VecWAXPY(W, -lambda, Y, X)); 48 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 49 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 50 if (linesearch->ops->vidirderiv) { 51 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty)); 52 } else { 53 PetscCall(VecDot(F, Y, &fty)); 54 } 55 56 delLambda = lambda - lambda_old; 57 58 /* check for convergence */ 59 if (PetscAbsReal(delLambda) < steptol * lambda) break; 60 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break; 61 if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break; 62 if (monitor) { 63 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 64 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old))); 65 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 66 } 67 68 /* compute the search direction */ 69 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { 70 s = (fty - fty_old) / delLambda; 71 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 72 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 73 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 74 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 75 if (linesearch->ops->vidirderiv) { 76 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1)); 77 } else { 78 PetscCall(VecDot(F, Y, &fty_mid1)); 79 } 80 s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda; 81 } else { 82 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 83 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 84 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 85 if (linesearch->ops->vidirderiv) { 86 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1)); 87 } else { 88 PetscCall(VecDot(F, Y, &fty_mid1)); 89 } 90 PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X)); 91 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 92 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 93 if (linesearch->ops->vidirderiv) { 94 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2)); 95 } else { 96 PetscCall(VecDot(F, Y, &fty_mid2)); 97 } 98 s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda); 99 } 100 /* if the solve is going in the wrong direction, fix it */ 101 if (PetscRealPart(s) > 0.) s = -s; 102 if (s == 0.0) break; 103 lambda_update = lambda - PetscRealPart(fty / s); 104 105 /* switch directions if we stepped out of bounds */ 106 if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s); 107 108 if (PetscIsInfOrNanReal(lambda_update)) break; 109 if (lambda_update > maxstep) break; 110 111 /* compute the new state of the line search */ 112 lambda_old = lambda; 113 lambda = lambda_update; 114 fty_old = fty; 115 } 116 /* construct the solution */ 117 PetscCall(VecWAXPY(W, -lambda, Y, X)); 118 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 119 /* postcheck */ 120 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 121 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 122 if (changed_y) { 123 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 124 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 125 } 126 PetscCall(VecCopy(W, X)); 127 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 128 129 PetscCall(SNESLineSearchComputeNorms(linesearch)); 130 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 131 132 if (monitor) { 133 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 134 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm)); 135 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 136 } 137 if (lambda <= steptol) PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /*MC 142 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 143 artificial $G(x)$ for which the `SNESFunctionFn` $ F(x) = grad G(x)$. Therefore, this line search seeks 144 to find roots of $ F^T Y$ via a secant method. 145 146 Options Database Keys: 147 + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda 148 . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value 149 . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor on entry to the line search, default is 1.0 150 - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed. 151 152 Level: advanced 153 154 Notes: 155 This method does NOT use the objective function if it is provided with `SNESSetObjective()`. 156 157 This method is the preferred line search for `SNESQN` and `SNESNCG`. 158 159 .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION` 160 M*/ 161 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 162 { 163 PetscFunctionBegin; 164 linesearch->ops->apply = SNESLineSearchApply_CP; 165 linesearch->ops->destroy = NULL; 166 linesearch->ops->setfromoptions = NULL; 167 linesearch->ops->reset = NULL; 168 linesearch->ops->view = NULL; 169 linesearch->ops->setup = NULL; 170 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 171 172 linesearch->max_its = 1; 173 PetscFunctionReturn(PETSC_SUCCESS); 174 } 175