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