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(0); 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(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w)); 101 if (changed_y) { 102 PetscCall(VecAXPY(X, -lambda, Y)); 103 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, X)); 104 } else { 105 PetscCall(VecCopy(W, X)); 106 } 107 PetscCall((*linesearch->ops->snesfunc)(snes,X,F)); 108 109 PetscCall(SNESLineSearchComputeNorms(linesearch)); 110 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 111 112 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 113 114 if (monitor) { 115 PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 116 PetscCall(PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm)); 117 PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 118 } 119 if (lambda <= steptol) { 120 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 /*MC 126 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 127 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 128 to find roots of dot(F, Y) via a secant method. 129 130 Options Database Keys: 131 + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda 132 . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value 133 . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0 134 - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed. 135 136 Notes: 137 This method does NOT use the objective function if it is provided with SNESSetObjective(). 138 139 This method is the preferred line search for SNESQN and SNESNCG. 140 141 Level: advanced 142 143 .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 144 M*/ 145 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 146 { 147 PetscFunctionBegin; 148 linesearch->ops->apply = SNESLineSearchApply_CP; 149 linesearch->ops->destroy = NULL; 150 linesearch->ops->setfromoptions = NULL; 151 linesearch->ops->reset = NULL; 152 linesearch->ops->view = NULL; 153 linesearch->ops->setup = NULL; 154 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 155 156 linesearch->max_its = 1; 157 PetscFunctionReturn(0); 158 } 159