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