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, Fold; 189 PetscReal fnorm, gnorm, dummyNorm, beta = 0.0; 190 PetscScalar dXdotF, dXolddotFold, dXdotFold, lXdotF, lXdotFold; 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 Fold = 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, F, &dXdotF);CHKERRQ(ierr); 254 255 for(i = 1; i < maxits + 1; i++) { 256 lsSuccess = PETSC_TRUE; 257 /* some update types require the old update direction or conjugate direction */ 258 if (ncg->betatype == 1 || ncg->betatype == 2 || ncg->betatype == 3 || ncg->betatype == 4) { /* prp, hs, dy, cd need dXold*/ 259 ierr = VecCopy(F, Fold);CHKERRQ(ierr); 260 } 261 ierr = SNESLineSearchPreCheckApply(snes, X, lX, PETSC_NULL);CHKERRQ(ierr); 262 ierr = SNESLineSearchApply(snes, X, F, lX, fnorm, 0.0, W, G, &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 dXolddotFold = dXdotF; 316 ierr = VecDot(F, dX, &dXdotF);CHKERRQ(ierr); 317 beta = PetscRealPart(dXdotF / dXolddotFold); 318 break; 319 case 1: /* Polak-Ribiere-Poylak */ 320 dXolddotFold = dXdotF; 321 ierr = VecDot(F, dX, &dXdotF);CHKERRQ(ierr); 322 ierr = VecDot(Fold, dX, &dXdotFold);CHKERRQ(ierr); 323 beta = PetscRealPart(((dXdotF - dXdotFold) / dXolddotFold)); 324 if (beta < 0.0) beta = 0.0; /* restart */ 325 break; 326 case 2: /* Hestenes-Stiefel */ 327 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 328 ierr = VecDot(dX, Fold, &dXdotFold);CHKERRQ(ierr); 329 ierr = VecDot(lX, F, &lXdotF);CHKERRQ(ierr); 330 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 331 beta = PetscRealPart((dXdotF - dXdotFold) / (lXdotF - lXdotFold)); 332 break; 333 case 3: /* Dai-Yuan */ 334 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 335 ierr = VecDot(lX, F, &lXdotF);CHKERRQ(ierr); 336 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 337 beta = PetscRealPart(dXdotF / (lXdotFold - lXdotF));CHKERRQ(ierr); 338 break; 339 case 4: /* Conjugate Descent */ 340 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 341 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 342 beta = PetscRealPart(dXdotF / lXdotFold);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 snes->max_funcs = 30000; 391 snes->max_its = 10000; 392 393 ierr = PetscNewLog(snes, SNES_NCG, &neP);CHKERRQ(ierr); 394 snes->data = (void*) neP; 395 neP->monitor = PETSC_NULL; 396 neP->betatype = 4; 397 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NCG",SNESLineSearchSetType_NCG);CHKERRQ(ierr); 398 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 EXTERN_C_END 402