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