1 #include <../src/snes/impls/ncg/snesncgimpl.h> 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "SNESReset_NCG" 5 PetscErrorCode SNESReset_NCG(SNES snes) 6 { 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 if (snes->work) {ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr);} 11 PetscFunctionReturn(0); 12 } 13 14 /* 15 SNESDestroy_NCG - Destroys the private SNES_NCG context that was created with SNESCreate_NCG(). 16 17 Input Parameter: 18 . snes - the SNES context 19 20 Application Interface Routine: SNESDestroy() 21 */ 22 #undef __FUNCT__ 23 #define __FUNCT__ "SNESDestroy_NCG" 24 PetscErrorCode SNESDestroy_NCG(SNES snes) 25 { 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 ierr = SNESReset_NCG(snes);CHKERRQ(ierr); 30 ierr = PetscFree(snes->data);CHKERRQ(ierr); 31 PetscFunctionReturn(0); 32 } 33 34 /* 35 SNESSetUp_NCG - Sets up the internal data structures for the later use 36 of the SNESNCG nonlinear solver. 37 38 Input Parameters: 39 + snes - the SNES context 40 - x - the solution vector 41 42 Application Interface Routine: SNESSetUp() 43 */ 44 #undef __FUNCT__ 45 #define __FUNCT__ "SNESSetUp_NCG" 46 PetscErrorCode SNESSetUp_NCG(SNES snes) 47 { 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 ierr = SNESDefaultGetWork(snes,4);CHKERRQ(ierr); 52 PetscFunctionReturn(0); 53 } 54 /* 55 SNESSetFromOptions_NCG - Sets various parameters for the SNESLS method. 56 57 Input Parameter: 58 . snes - the SNES context 59 60 Application Interface Routine: SNESSetFromOptions() 61 */ 62 #undef __FUNCT__ 63 #define __FUNCT__ "SNESSetFromOptions_NCG" 64 static PetscErrorCode SNESSetFromOptions_NCG(SNES snes) 65 { 66 SNES_NCG *ncg = (SNES_NCG *)snes->data; 67 const char *betas[] = {"FR","PRP","HS", "DY", "CD"}; 68 PetscErrorCode ierr; 69 PetscBool debug, flg; 70 PetscInt indx; 71 72 PetscFunctionBegin; 73 ierr = PetscOptionsHead("SNES NCG options");CHKERRQ(ierr); 74 ierr = PetscOptionsBool("-snes_ncg_monitor", "Monitor NCG iterations", "SNES", ncg->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 75 ierr = PetscOptionsEList("-snes_ncg_type","Nonlinear CG update used","", betas, 5, "FR",&indx, &flg);CHKERRQ(ierr); 76 if (flg) { 77 ncg->betatype = indx; 78 } 79 if (debug) { 80 ncg->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 81 } 82 ierr = PetscOptionsTail();CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 /* 87 SNESView_NCG - Prints info from the SNESNCG data structure. 88 89 Input Parameters: 90 + SNES - the SNES context 91 - viewer - visualization context 92 93 Application Interface Routine: SNESView() 94 */ 95 #undef __FUNCT__ 96 #define __FUNCT__ "SNESView_NCG" 97 static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer) 98 { 99 PetscBool iascii; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 104 if (iascii) { 105 ierr = PetscViewerASCIIPrintf(viewer," line search type variant: %s\n", SNESLineSearchTypeName(snes->ls_type));CHKERRQ(ierr); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "SNESLineSearchExactLinear_NCG" 112 PetscErrorCode SNESLineSearchExactLinear_NCG(SNES snes,void *lsctx,Vec X,Vec F,Vec Y,PetscReal fnorm,PetscReal dummyXnorm,Vec G,Vec W,PetscReal *dummyYnorm,PetscReal *gnorm,PetscBool *flag) 113 { 114 PetscScalar alpha, ptAp; 115 PetscErrorCode ierr; 116 /* SNES_NCG *ncg = (SNES_NCG *) snes->data; */ 117 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 118 119 PetscFunctionBegin; 120 /* 121 122 The exact step size for unpreconditioned linear CG is just: 123 124 alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy) 125 126 */ 127 128 ierr = SNESComputeJacobian(snes, X, &snes->jacobian, &snes->jacobian_pre, &flg);CHKERRQ(ierr); 129 ierr = VecDot(F, F, &alpha);CHKERRQ(ierr); 130 ierr = MatMult(snes->jacobian, Y, W);CHKERRQ(ierr); 131 ierr = VecDot(Y, W, &ptAp);CHKERRQ(ierr); 132 alpha = alpha / ptAp; 133 ierr = PetscPrintf(((PetscObject)snes)->comm, "alpha: %G\n", PetscRealPart(alpha));CHKERRQ(ierr); 134 ierr = VecCopy(X, W);CHKERRQ(ierr); 135 ierr = VecAXPY(W, alpha, Y);CHKERRQ(ierr); 136 ierr = SNESComputeFunction(snes, W, G);CHKERRQ(ierr); 137 ierr = VecNorm(G, NORM_2, gnorm);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 EXTERN_C_BEGIN 142 #undef __FUNCT__ 143 #define __FUNCT__ "SNESLineSearchSetType_NCG" 144 PetscErrorCode SNESLineSearchSetType_NCG(SNES snes, SNESLineSearchType type) 145 { 146 PetscErrorCode ierr; 147 PetscFunctionBegin; 148 149 switch (type) { 150 case SNES_LS_BASIC: 151 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 152 break; 153 case SNES_LS_BASIC_NONORMS: 154 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 155 break; 156 case SNES_LS_QUADRATIC: 157 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 158 break; 159 case SNES_LS_TEST: 160 ierr = SNESLineSearchSet(snes,SNESLineSearchExactLinear_NCG,PETSC_NULL);CHKERRQ(ierr); 161 break; 162 case SNES_LS_SECANT: 163 ierr = SNESLineSearchSet(snes,SNESLineSearchSecant,PETSC_NULL);CHKERRQ(ierr); 164 break; 165 default: 166 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type"); 167 break; 168 } 169 snes->ls_type = type; 170 PetscFunctionReturn(0); 171 } 172 EXTERN_C_END 173 174 /* 175 SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method. 176 177 Input Parameters: 178 . snes - the SNES context 179 180 Output Parameter: 181 . outits - number of iterations until termination 182 183 Application Interface Routine: SNESSolve() 184 */ 185 #undef __FUNCT__ 186 #define __FUNCT__ "SNESSolve_NCG" 187 PetscErrorCode SNESSolve_NCG(SNES snes) 188 { 189 SNES_NCG *ncg = (SNES_NCG *)snes->data; 190 Vec X, dX, lX, F, W, B, G, dXold; 191 PetscReal fnorm, gnorm, dummyNorm, beta = 0.0; 192 PetscScalar dXdotdX, dXolddotdXold, dXdotdXold, lXdotdX, lXdotdXold; 193 PetscInt maxits, i; 194 PetscErrorCode ierr; 195 SNESConvergedReason reason; 196 PetscBool lsSuccess = PETSC_TRUE; 197 PetscFunctionBegin; 198 snes->reason = SNES_CONVERGED_ITERATING; 199 200 maxits = snes->max_its; /* maximum number of iterations */ 201 X = snes->vec_sol; /* X^n */ 202 dXold = snes->work[0]; /* The previous iterate of X */ 203 dX = snes->work[1]; /* the steepest direction */ 204 lX = snes->vec_sol_update; /* the conjugate direction */ 205 F = snes->vec_func; /* residual vector */ 206 W = snes->work[2]; /* work vector and solution output for the line search */ 207 B = snes->vec_rhs; /* the right hand side */ 208 G = snes->work[3]; /* the line search result function evaluation */ 209 210 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 211 snes->iter = 0; 212 snes->norm = 0.; 213 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 214 215 /* compute the initial function and preconditioned update delX */ 216 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 217 if (snes->domainerror) { 218 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 219 PetscFunctionReturn(0); 220 } 221 /* convergence test */ 222 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 223 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 224 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 225 snes->norm = fnorm; 226 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 227 SNESLogConvHistory(snes,fnorm,0); 228 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 229 230 /* set parameter for default relative tolerance convergence test */ 231 snes->ttol = fnorm*snes->rtol; 232 /* test convergence */ 233 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 234 if (snes->reason) PetscFunctionReturn(0); 235 236 /* Call general purpose update function */ 237 if (snes->ops->update) { 238 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 239 } 240 241 /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */ 242 if (snes->pc) { 243 ierr = VecCopy(X, dX);CHKERRQ(ierr); 244 ierr = SNESSolve(snes->pc, B, dX);CHKERRQ(ierr); 245 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 246 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 247 snes->reason = SNES_DIVERGED_INNER; 248 PetscFunctionReturn(0); 249 } 250 ierr = VecAYPX(dX,-1.0,X);CHKERRQ(ierr); 251 } else { 252 ierr = VecCopy(F, dX);CHKERRQ(ierr); 253 } 254 ierr = VecCopy(dX, lX);CHKERRQ(ierr); 255 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 256 for(i = 1; i < maxits + 1; i++) { 257 lsSuccess = PETSC_TRUE; 258 /* some update types require the old update direction or conjugate direction */ 259 if (ncg->betatype == 1 || ncg->betatype == 2 || ncg->betatype == 3 || ncg->betatype == 4) { /* prp, hs, dy, cd need dXold*/ 260 ierr = VecCopy(dX, dXold);CHKERRQ(ierr); 261 } 262 ierr = (*snes->ops->linesearch)(snes, snes->lsP, X, F, lX, fnorm, 0.0, G, W, &dummyNorm, &gnorm, &lsSuccess);CHKERRQ(ierr); 263 if (!lsSuccess) { 264 if (++snes->numFailures >= snes->maxFailures) { 265 snes->reason = SNES_DIVERGED_LINE_SEARCH; 266 PetscFunctionReturn(0); 267 } 268 } 269 if (snes->nfuncs >= snes->max_funcs) { 270 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 271 PetscFunctionReturn(0); 272 } 273 if (snes->domainerror) { 274 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 275 PetscFunctionReturn(0); 276 } 277 278 /* copy over the solution */ 279 ierr = VecCopy(G, F);CHKERRQ(ierr); 280 ierr = VecCopy(W, X);CHKERRQ(ierr); 281 fnorm = gnorm; 282 283 /* Monitor convergence */ 284 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 285 snes->iter = i; 286 snes->norm = fnorm; 287 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 288 SNESLogConvHistory(snes,snes->norm,0); 289 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 290 291 /* Test for convergence */ 292 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 293 if (snes->reason) PetscFunctionReturn(0); 294 295 /* Call general purpose update function */ 296 if (snes->ops->update) { 297 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 298 } 299 if (snes->pc) { 300 ierr = VecCopy(X,dX);CHKERRQ(ierr); 301 ierr = SNESSolve(snes->pc, B, dX);CHKERRQ(ierr); 302 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 303 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 304 snes->reason = SNES_DIVERGED_INNER; 305 PetscFunctionReturn(0); 306 } 307 ierr = VecAYPX(dX,-1.0,X);CHKERRQ(ierr); 308 } else { 309 ierr = VecCopy(F,dX);CHKERRQ(ierr); 310 } 311 312 /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/ 313 switch(ncg->betatype) { 314 case 0: /* Fletcher-Reeves */ 315 dXolddotdXold = dXdotdX; 316 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 317 beta = PetscRealPart(dXdotdX / dXolddotdXold); 318 break; 319 case 1: /* Polak-Ribiere-Poylak */ 320 dXolddotdXold = dXdotdX; 321 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 322 ierr = VecDot(dXold, dX, &dXdotdXold);CHKERRQ(ierr); 323 beta = PetscRealPart(((dXdotdX - dXdotdXold) / dXolddotdXold)); 324 if (beta < 0.0) beta = 0.0; /* restart */ 325 break; 326 case 2: /* Hestenes-Stiefel */ 327 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 328 ierr = VecDot(dX, dXold, &dXdotdXold);CHKERRQ(ierr); 329 ierr = VecDot(lX, dX, &lXdotdX);CHKERRQ(ierr); 330 ierr = VecDot(lX, dXold, &lXdotdXold);CHKERRQ(ierr); 331 beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold)); 332 break; 333 case 3: /* Dai-Yuan */ 334 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 335 ierr = VecDot(lX, dX, &lXdotdXold);CHKERRQ(ierr); 336 ierr = VecDot(lX, dXold, &lXdotdXold);CHKERRQ(ierr); 337 beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX));CHKERRQ(ierr); 338 break; 339 case 4: /* Conjugate Descent */ 340 ierr = VecDot(lX, dXold, &lXdotdXold);CHKERRQ(ierr); 341 ierr = VecDot(dX, dX, &dXdotdX);CHKERRQ(ierr); 342 beta = PetscRealPart(dXdotdX / lXdotdXold);CHKERRQ(ierr); 343 break; 344 } 345 if (ncg->monitor) { 346 ierr = PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", beta);CHKERRQ(ierr); 347 } 348 ierr = VecAYPX(lX, beta, dX);CHKERRQ(ierr); 349 } 350 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr); 351 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 352 PetscFunctionReturn(0); 353 } 354 355 /*MC 356 SNESNCG - Nonlinear Conjugate-Gradient method for the solution of general nonlinear systems. 357 358 Level: beginner 359 360 Options Database: 361 + -snes_ls_damping - damping factor to apply to F(x) (used only if -snes_ls is basic or basicnonorms) 362 - -snes_ls <basic,basicnormnorms,quadratic> 363 364 Notes: This solves the nonlinear system of equations F(x) = 0 using the nonlinear generalization of the conjugate 365 gradient method. This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise 366 chooses the initial search direction as F(x) for the initial guess x. 367 368 369 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR, SNESNGMRES, SNESNQN 370 M*/ 371 EXTERN_C_BEGIN 372 #undef __FUNCT__ 373 #define __FUNCT__ "SNESCreate_NCG" 374 PetscErrorCode SNESCreate_NCG(SNES snes) 375 { 376 PetscErrorCode ierr; 377 SNES_NCG * neP; 378 379 PetscFunctionBegin; 380 snes->ops->destroy = SNESDestroy_NCG; 381 snes->ops->setup = SNESSetUp_NCG; 382 snes->ops->setfromoptions = SNESSetFromOptions_NCG; 383 snes->ops->view = SNESView_NCG; 384 snes->ops->solve = SNESSolve_NCG; 385 snes->ops->reset = SNESReset_NCG; 386 387 snes->usesksp = PETSC_FALSE; 388 snes->usespc = PETSC_TRUE; 389 390 ierr = PetscNewLog(snes, SNES_NCG, &neP);CHKERRQ(ierr); 391 snes->data = (void*) neP; 392 neP->monitor = PETSC_NULL; 393 neP->betatype = 4; 394 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NCG",SNESLineSearchSetType_NCG);CHKERRQ(ierr); 395 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398 EXTERN_C_END 399