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