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, domainerror; 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 22 ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, PETSC_NULL);CHKERRQ(ierr); 23 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 24 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 25 ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 26 ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);CHKERRQ(ierr); 27 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 28 ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 29 30 /* precheck */ 31 ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr); 32 lambda_old = 0.0; 33 ierr = VecDot(F, Y, &fty_old);CHKERRQ(ierr); 34 fty_init = fty_old; 35 36 for (i = 0; i < max_its; i++) { 37 38 /* compute the norm at lambda */ 39 ierr = VecCopy(X, W);CHKERRQ(ierr); 40 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 41 if (linesearch->ops->viproject) { 42 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 43 } 44 ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr); 45 46 ierr = VecDot(F, Y, &fty);CHKERRQ(ierr); 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 && i > 0) break; 54 if (monitor) { 55 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 56 ierr = PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", 57 lambda, lambda_old, PetscRealPart(fty), PetscRealPart(fty_old));CHKERRQ(ierr); 58 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 59 } 60 61 /* compute the search direction */ 62 if (linesearch->order == SNES_LINESEARCH_LINEAR) { 63 s = (fty - fty_old) / delLambda; 64 } else if (linesearch->order == SNES_LINESEARCH_QUADRATIC) { 65 ierr = VecCopy(X, W);CHKERRQ(ierr); 66 ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr); 67 if (linesearch->ops->viproject) { 68 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 69 } 70 ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr); 71 ierr = VecDot(F, Y, &fty_mid1);CHKERRQ(ierr); 72 s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda; 73 } else { 74 ierr = VecCopy(X, W);CHKERRQ(ierr); 75 ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr); 76 if (linesearch->ops->viproject) { 77 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 78 } 79 ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr); 80 ierr = VecDot(F, Y, &fty_mid1);CHKERRQ(ierr); 81 ierr = VecCopy(X, W);CHKERRQ(ierr); 82 ierr = VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);CHKERRQ(ierr); 83 if (linesearch->ops->viproject) { 84 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 85 } 86 ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr); 87 ierr = VecDot(F, Y, &fty_mid2);CHKERRQ(ierr); 88 s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda); 89 } 90 /* if the solve is going in the wrong direction, fix it */ 91 if (PetscRealPart(s) > 0.) s = -s; 92 lambda_update = lambda - PetscRealPart(fty / s); 93 94 /* switch directions if we stepped out of bounds */ 95 if (lambda_update < steptol) { 96 lambda_update = lambda + PetscRealPart(fty / s); 97 } 98 99 if (PetscIsInfOrNanScalar(lambda_update)) break; 100 if (lambda_update > maxstep) { 101 break; 102 } 103 104 /* compute the new state of the line search */ 105 lambda_old = lambda; 106 lambda = lambda_update; 107 fty_old = fty; 108 } 109 /* construct the solution */ 110 ierr = VecCopy(X, W);CHKERRQ(ierr); 111 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 112 if (linesearch->ops->viproject) { 113 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 114 } 115 /* postcheck */ 116 ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 117 if (changed_y) { 118 ierr = VecAXPY(X, -lambda, Y);CHKERRQ(ierr); 119 if (linesearch->ops->viproject) { 120 ierr = (*linesearch->ops->viproject)(snes, X);CHKERRQ(ierr); 121 } 122 } else { 123 ierr = VecCopy(W, X);CHKERRQ(ierr); 124 } 125 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 126 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 127 if (domainerror) { 128 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 132 ierr = SNESLineSearchComputeNorms(linesearch);CHKERRQ(ierr); 133 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 134 135 ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 136 137 if (monitor) { 138 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 139 ierr = PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", lambda, gnorm);CHKERRQ(ierr); 140 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 141 } 142 if (lambda <= steptol) { 143 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 144 } 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNCT__ 149 #define __FUNCT__ "SNESLineSearchCreate_CP" 150 /*MC 151 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 152 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 153 to find roots of f^ty via a secant method. 154 155 Options Database Keys: 156 . -snes_linesearch_damping - initial trial step length 157 158 Notes: 159 This method is the preferred line search for SNESQN and SNESNCG. 160 161 Level: advanced 162 163 .keywords: SNES, SNESLineSearch, damping 164 165 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 166 M*/ 167 PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 168 { 169 PetscFunctionBegin; 170 linesearch->ops->apply = SNESLineSearchApply_CP; 171 linesearch->ops->destroy = PETSC_NULL; 172 linesearch->ops->setfromoptions = PETSC_NULL; 173 linesearch->ops->reset = PETSC_NULL; 174 linesearch->ops->view = PETSC_NULL; 175 linesearch->ops->setup = PETSC_NULL; 176 linesearch->order = SNES_LINESEARCH_LINEAR; 177 PetscFunctionReturn(0); 178 } 179