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, Fold; 191 PetscReal fnorm, gnorm, dummyNorm, beta = 0.0; 192 PetscScalar dXdotF, dXolddotFold, dXdotFold, lXdotF, lXdotFold; 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 Fold = 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(F, lX);CHKERRQ(ierr); 244 ierr = VecScale(lX, -1.0);CHKERRQ(ierr); 245 } else { 246 ierr = VecCopy(X, lX);CHKERRQ(ierr); 247 ierr = SNESSolve(snes->pc, B, lX);CHKERRQ(ierr); 248 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 249 ierr = VecAXPY(lX, -1.0, X);CHKERRQ(ierr); 250 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 251 snes->reason = SNES_DIVERGED_INNER; 252 PetscFunctionReturn(0); 253 } 254 } 255 ierr = VecDot(lX, F, &dXdotF);CHKERRQ(ierr); 256 for(i = 1; i < maxits; i++) { 257 lsSuccess = PETSC_TRUE; 258 /* some update types require the old residual */ 259 ierr = VecCopy(F, Fold);CHKERRQ(ierr); 260 ierr = (*snes->ops->linesearch)(snes, snes->lsP, X, F, lX, fnorm, 0.0, G, W, &dummyNorm, &gnorm, &lsSuccess);CHKERRQ(ierr); 261 if (!lsSuccess) { 262 if (++snes->numFailures >= snes->maxFailures) { 263 snes->reason = SNES_DIVERGED_LINE_SEARCH; 264 break; 265 } 266 } 267 if (snes->nfuncs >= snes->max_funcs) { 268 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 269 break; 270 } 271 if (snes->domainerror) { 272 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 273 PetscFunctionReturn(0); 274 } 275 276 /* copy over the solution */ 277 ierr = VecCopy(G, F);CHKERRQ(ierr); 278 ierr = VecCopy(W, X);CHKERRQ(ierr); 279 fnorm = gnorm; 280 281 /* Monitor convergence */ 282 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 283 snes->iter = i; 284 snes->norm = fnorm; 285 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 286 SNESLogConvHistory(snes,snes->norm,0); 287 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 288 289 /* Test for convergence */ 290 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 291 if (snes->reason) break; 292 293 /* Call general purpose update function */ 294 if (snes->ops->update) { 295 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 296 } 297 if (!snes->pc) { 298 ierr = VecCopy(F,dX);CHKERRQ(ierr); 299 ierr = VecScale(dX,-1.0);CHKERRQ(ierr); 300 } else { 301 ierr = VecCopy(X,dX);CHKERRQ(ierr); 302 ierr = SNESSolve(snes->pc, B, dX);CHKERRQ(ierr); 303 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 304 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 305 snes->reason = SNES_DIVERGED_INNER; 306 PetscFunctionReturn(0); 307 } 308 ierr = VecAXPY(dX,-1.0,X);CHKERRQ(ierr); 309 } 310 311 /* compute the conjugate direction lX = dX + beta*lX with beta = ((F, dX) / (F_old, dX_old) (Fletcher-Reeves update)*/ 312 switch(ncg->betatype) { 313 case 0: /* Fletcher-Reeves */ 314 dXolddotFold = dXdotF; 315 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 316 beta = PetscRealPart(dXdotF / dXolddotFold); 317 break; 318 case 1: /* Polak-Ribiere-Poylak */ 319 dXolddotFold = dXdotF; 320 ierr = VecDot(dX, Fold, &dXdotFold);CHKERRQ(ierr); 321 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 322 beta = PetscRealPart(((dXdotF - dXdotFold) / dXolddotFold)); 323 if (beta < 0.0) beta = 0.0; /* restart */ 324 break; 325 case 2: /* Hestenes-Stiefel */ 326 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 327 ierr = VecDot(dX, Fold, &dXdotFold);CHKERRQ(ierr); 328 ierr = VecDot(lX, F, &lXdotF);CHKERRQ(ierr); 329 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 330 beta = PetscRealPart((dXdotF - dXdotFold) / (lXdotF - lXdotFold)); 331 break; 332 case 3: /* Dai-Yuan */ 333 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 334 ierr = VecDot(lX, F, &lXdotF);CHKERRQ(ierr); 335 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 336 beta = PetscRealPart(dXdotF / (lXdotFold - lXdotF));CHKERRQ(ierr); 337 break; 338 case 4: /* Conjugate Descent */ 339 ierr = VecDot(dX, F, &dXdotF);CHKERRQ(ierr); 340 ierr = VecDot(lX, Fold, &lXdotFold);CHKERRQ(ierr); 341 beta = PetscRealPart(dXdotF / lXdotFold);CHKERRQ(ierr); 342 break; 343 } 344 if (ncg->monitor) { 345 ierr = PetscPrintf(PETSC_COMM_WORLD, "beta = %e\n", beta);CHKERRQ(ierr); 346 } 347 ierr = VecAYPX(lX, beta, dX);CHKERRQ(ierr); 348 /* line search for the proper contribution of lX to the solution */ 349 } 350 if (i == maxits) { 351 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);CHKERRQ(ierr); 352 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 353 } 354 PetscFunctionReturn(0); 355 } 356 357 /*MC 358 SNESNCG - Nonlinear Conjugate-Gradient method for the solution of general nonlinear systems. 359 360 Level: beginner 361 362 Options Database: 363 + -snes_ls_damping - damping factor to apply to F(x) (used only if -snes_ls is basic or basicnonorms) 364 - -snes_ls <basic,basicnormnorms,quadratic> 365 366 Notes: This solves the nonlinear system of equations F(x) = 0 using the nonlinear generalization of the conjugate 367 gradient method. This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise 368 chooses the initial search direction as F(x) for the initial guess x. 369 370 371 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR, SNESNGMRES, SNESNQN 372 M*/ 373 EXTERN_C_BEGIN 374 #undef __FUNCT__ 375 #define __FUNCT__ "SNESCreate_NCG" 376 PetscErrorCode SNESCreate_NCG(SNES snes) 377 { 378 PetscErrorCode ierr; 379 SNES_NCG * neP; 380 381 PetscFunctionBegin; 382 snes->ops->destroy = SNESDestroy_NCG; 383 snes->ops->setup = SNESSetUp_NCG; 384 snes->ops->setfromoptions = SNESSetFromOptions_NCG; 385 snes->ops->view = SNESView_NCG; 386 snes->ops->solve = SNESSolve_NCG; 387 snes->ops->reset = SNESReset_NCG; 388 389 snes->usesksp = PETSC_FALSE; 390 snes->usespc = PETSC_TRUE; 391 392 ierr = PetscNewLog(snes, SNES_NCG, &neP);CHKERRQ(ierr); 393 snes->data = (void*) neP; 394 neP->monitor = PETSC_NULL; 395 neP->betatype = 0; 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