1 /* Defines the basic SNES object */ 2 #include <../src/snes/impls/ngmres/snesngmres.h> 3 #include <petscblaslapack.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "SNESReset_NGMRES" 7 PetscErrorCode SNESReset_NGMRES(SNES snes) 8 { 9 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = VecDestroyVecs(ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr); 14 ierr = VecDestroyVecs(ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "SNESDestroy_NGMRES" 20 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 21 { 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 26 if (snes->data) { 27 SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data; 28 ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->fnorms, ngmres->q);CHKERRQ(ierr); 29 ierr = PetscFree(ngmres->s);CHKERRQ(ierr); 30 #if PETSC_USE_COMPLEX 31 ierr = PetscFree(ngmres->rwork); 32 #endif 33 ierr = PetscFree(ngmres->work); 34 } 35 ierr = PetscFree(snes->data); 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "SNESSetUp_NGMRES" 41 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 42 { 43 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 44 PetscInt msize,hsize; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 msize = ngmres->msize; /* restart size */ 49 hsize = msize * msize; 50 51 52 /* explicit least squares minimization solve */ 53 ierr = PetscMalloc5(hsize,PetscScalar,&ngmres->h, 54 msize,PetscScalar,&ngmres->beta, 55 msize,PetscScalar,&ngmres->xi, 56 msize,PetscReal, &ngmres->fnorms, 57 hsize,PetscScalar,&ngmres->q);CHKERRQ(ierr); 58 ngmres->nrhs = 1; 59 ngmres->lda = msize; 60 ngmres->ldb = msize; 61 ierr = PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);CHKERRQ(ierr); 62 ierr = PetscMemzero(ngmres->h, hsize*sizeof(PetscScalar));CHKERRQ(ierr); 63 ierr = PetscMemzero(ngmres->q, hsize*sizeof(PetscScalar));CHKERRQ(ierr); 64 ierr = PetscMemzero(ngmres->xi, msize*sizeof(PetscScalar));CHKERRQ(ierr); 65 ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr); 66 ngmres->lwork = 12*msize; 67 #if PETSC_USE_COMPLEX 68 ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork); 69 #endif 70 ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work); 71 72 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr); 73 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr); 74 ierr = SNESDefaultGetWork(snes, 5);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 80 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 81 { 82 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 83 PetscErrorCode ierr; 84 PetscBool debug; 85 86 PetscFunctionBegin; 87 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 88 ierr = PetscOptionsBool("-snes_ngmres_additive", "Use additive variant vs. choice", "SNES", ngmres->additive, &ngmres->additive, PETSC_NULL);CHKERRQ(ierr); 89 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 90 ierr = PetscOptionsInt("-snes_ngmres_restart_it","Tolerance iterations before restart","SNES", ngmres->restart_it, &ngmres->restart_it, PETSC_NULL);CHKERRQ(ierr); 91 ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 92 if (debug) { 93 ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 94 } 95 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscOptionsReal("-snes_ngmres_gammaC", " Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 98 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 99 ierr = PetscOptionsTail();CHKERRQ(ierr); 100 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "SNESView_NGMRES" 106 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 107 { 108 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 109 PetscBool iascii; 110 const char *cstr; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 115 if (iascii) { 116 cstr = SNESLineSearchTypeName(snes->ls_type); 117 ierr = PetscViewerASCIIPrintf(viewer, " line search variant: %s\n",cstr);CHKERRQ(ierr); 118 ierr = PetscViewerASCIIPrintf(viewer, " Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer, " Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(viewer, " Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 EXTERN_C_BEGIN 126 #undef __FUNCT__ 127 #define __FUNCT__ "SNESLineSearchSetType_NGMRES" 128 PetscErrorCode SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type) 129 { 130 PetscErrorCode ierr; 131 PetscFunctionBegin; 132 133 switch (type) { 134 case SNES_LS_BASIC: 135 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 136 break; 137 case SNES_LS_BASIC_NONORMS: 138 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 139 break; 140 case SNES_LS_QUADRATIC: 141 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 142 break; 143 default: 144 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type."); 145 break; 146 } 147 snes->ls_type = type; 148 PetscFunctionReturn(0); 149 } 150 EXTERN_C_END 151 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "SNESSolve_NGMRES" 155 156 PetscErrorCode SNESSolve_NGMRES(SNES snes) 157 { 158 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 159 /* present solution, residual, and preconditioned residual */ 160 Vec X, F, B, D, Y; 161 162 /* candidate linear combination answers */ 163 Vec XA, FA, XM, FM; 164 165 /* previous iterations to construct the subspace */ 166 Vec *Fdot = ngmres->Fdot; 167 Vec *Xdot = ngmres->Xdot; 168 169 /* coefficients and RHS to the minimization problem */ 170 PetscScalar *beta = ngmres->beta; 171 PetscScalar *xi = ngmres->xi; 172 PetscReal fnorm, fMnorm, fAnorm, ynorm, xnorm = 0.0; 173 PetscReal nu; 174 PetscScalar alph_total = 0.; 175 PetscScalar qentry; 176 PetscInt i, j, k, k_restart, l, ivec, restart_count = 0; 177 178 /* solution selection data */ 179 PetscBool selectA, selectRestart; 180 PetscReal dnorm, dminnorm, dcurnorm; 181 PetscReal fminnorm; 182 183 SNESConvergedReason reason; 184 PetscBool lssucceed; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 /* variable initialization */ 189 snes->reason = SNES_CONVERGED_ITERATING; 190 X = snes->vec_sol; 191 F = snes->vec_func; 192 B = snes->vec_rhs; 193 XA = snes->vec_sol_update; 194 FA = snes->work[0]; 195 D = snes->work[1]; 196 197 /* work for the line search */ 198 Y = snes->work[2]; 199 XM = snes->work[3]; 200 FM = snes->work[4]; 201 202 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 203 snes->iter = 0; 204 snes->norm = 0.; 205 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 206 207 /* initialization */ 208 209 /* r = F(x) */ 210 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 211 if (snes->domainerror) { 212 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 213 PetscFunctionReturn(0); 214 } 215 216 /* nu = (r, r) */ 217 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 218 fminnorm = fnorm; 219 nu = fnorm*fnorm; 220 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 221 222 /* q_{00} = nu */ 223 Q(0,0) = nu; 224 ngmres->fnorms[0] = fnorm; 225 /* Fdot[0] = F */ 226 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 227 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 228 229 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 230 snes->norm = fnorm; 231 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 232 SNESLogConvHistory(snes, fnorm, 0); 233 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 234 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 235 if (snes->reason) PetscFunctionReturn(0); 236 237 k_restart = 1; 238 l = 1; 239 for (k=1; k < snes->max_its+1; k++) { 240 /* select which vector of the stored subspace will be updated */ 241 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 242 243 /* Computation of x^M */ 244 if (snes->pc) { 245 ierr = VecCopy(X, XM);CHKERRQ(ierr); 246 ierr = SNESSolve(snes->pc, B, XM);CHKERRQ(ierr); 247 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 248 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 249 snes->reason = SNES_DIVERGED_INNER; 250 PetscFunctionReturn(0); 251 } 252 ierr = SNESComputeFunction(snes, XM, FM);CHKERRQ(ierr); 253 ierr = VecNorm(FM, NORM_2, &fMnorm);CHKERRQ(ierr); 254 } else { 255 /* no preconditioner -- just take gradient descent with line search */ 256 ierr = VecCopy(F, Y);CHKERRQ(ierr); 257 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FM,XM,&ynorm,&fMnorm,&lssucceed);CHKERRQ(ierr); 258 if (!lssucceed) { 259 if (++snes->numFailures >= snes->maxFailures) { 260 snes->reason = SNES_DIVERGED_LINE_SEARCH; 261 PetscFunctionReturn(0); 262 } 263 } 264 } 265 266 /* r = F(x) */ 267 nu = fMnorm*fMnorm; 268 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */ 269 270 /* construct the right hand side and xi factors */ 271 ierr = VecMDot(FM, l, Fdot, xi);CHKERRQ(ierr); 272 273 for (i = 0; i < l; i++) { 274 beta[i] = nu - xi[i]; 275 } 276 277 /* construct h */ 278 for (j = 0; j < l; j++) { 279 for (i = 0; i < l; i++) { 280 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 281 } 282 } 283 284 if(l == 1) { 285 /* simply set alpha[0] = beta[0] / H[0, 0] */ 286 beta[0] = beta[0] / H(0, 0); 287 } else { 288 #ifdef PETSC_MISSING_LAPACK_GELSS 289 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 290 #else 291 ngmres->m = PetscBLASIntCast(l); 292 ngmres->n = PetscBLASIntCast(l); 293 ngmres->info = PetscBLASIntCast(0); 294 ngmres->rcond = -1.; 295 #ifdef PETSC_USE_COMPLEX 296 LAPACKgelss_(&ngmres->m, 297 &ngmres->n, 298 &ngmres->nrhs, 299 ngmres->h, 300 &ngmres->lda, 301 ngmres->beta, 302 &ngmres->ldb, 303 ngmres->s, 304 &ngmres->rcond, 305 &ngmres->rank, 306 ngmres->work, 307 &ngmres->lwork, 308 ngmres->rwork, 309 &ngmres->info); 310 #else 311 LAPACKgelss_(&ngmres->m, 312 &ngmres->n, 313 &ngmres->nrhs, 314 ngmres->h, 315 &ngmres->lda, 316 ngmres->beta, 317 &ngmres->ldb, 318 ngmres->s, 319 &ngmres->rcond, 320 &ngmres->rank, 321 ngmres->work, 322 &ngmres->lwork, 323 &ngmres->info); 324 #endif 325 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 326 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 327 #endif 328 } 329 330 alph_total = 0.; 331 for (i = 0; i < l; i++) { 332 alph_total += beta[i]; 333 } 334 335 ierr = VecCopy(XM, XA);CHKERRQ(ierr); 336 ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr); 337 338 ierr = VecMAXPY(XA, l, beta, Xdot);CHKERRQ(ierr); 339 ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr); 340 ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr); 341 342 /* differences for selection and restart */ 343 ierr=VecCopy(XA, D);CHKERRQ(ierr); 344 ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr); 345 ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr); 346 dminnorm = -1.0; 347 for(i=0;i<l;i++) { 348 ierr=VecCopy(XA, D);CHKERRQ(ierr); 349 ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr); 350 ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr); 351 if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm; 352 } 353 354 /* combination (additive) or selection (multiplicative) of the N-GMRES solution */ 355 if (ngmres->additive) { 356 /* X = X + \lambda(XA - X) */ 357 if (ngmres->monitor) { 358 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);CHKERRQ(ierr); 359 } 360 ierr = VecCopy(XA, Y);CHKERRQ(ierr); 361 ierr = VecAYPX(Y, -1.0, X);CHKERRQ(ierr); 362 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FA,XA,&ynorm,&fnorm,&lssucceed);CHKERRQ(ierr); 363 if (!lssucceed) { 364 if (++snes->numFailures >= snes->maxFailures) { 365 snes->reason = SNES_DIVERGED_LINE_SEARCH; 366 PetscFunctionReturn(0); 367 } 368 } 369 if (ngmres->monitor) { 370 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", fnorm);CHKERRQ(ierr); 371 } 372 ierr = VecCopy(FA, F);CHKERRQ(ierr); 373 ierr = VecCopy(XA, X);CHKERRQ(ierr); 374 } else { 375 selectA = PETSC_TRUE; 376 /* Conditions for choosing the accelerated answer */ 377 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 378 if (fAnorm >= ngmres->gammaA*fminnorm) { 379 selectA = PETSC_FALSE; 380 } 381 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 382 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 383 } else { 384 selectA=PETSC_FALSE; 385 } 386 if (selectA) { 387 if (ngmres->monitor) { 388 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 389 } 390 /* copy it over */ 391 fnorm = fAnorm; 392 nu = fnorm*fnorm; 393 ierr = VecCopy(FA, F);CHKERRQ(ierr); 394 ierr = VecCopy(XA, X);CHKERRQ(ierr); 395 } else { 396 if (ngmres->monitor) { 397 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 398 } 399 fnorm = fMnorm; 400 nu = fnorm*fnorm; 401 ierr = VecCopy(FM, F);CHKERRQ(ierr); 402 ierr = VecCopy(XM, X);CHKERRQ(ierr); 403 } 404 } 405 406 selectRestart = PETSC_FALSE; 407 408 /* difference stagnation restart */ 409 if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) { 410 if (ngmres->monitor) { 411 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr); 412 } 413 selectRestart = PETSC_TRUE; 414 } 415 /* residual stagnation restart */ 416 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 417 if (ngmres->monitor) { 418 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 419 } 420 selectRestart = PETSC_TRUE; 421 } 422 423 /* if the restart conditions persist for more than restart_it iterations, restart. */ 424 if (selectRestart) { 425 restart_count++; 426 } else { 427 restart_count = 0; 428 } 429 430 /* restart after restart conditions have persisted for a fixed number of iterations */ 431 if (restart_count >= ngmres->restart_it) { 432 if (ngmres->monitor){ 433 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 434 } 435 restart_count = 0; 436 k_restart = 1; 437 l = 1; 438 /* q_{00} = nu */ 439 ngmres->fnorms[0] = fnorm; 440 nu = fnorm*fnorm; 441 Q(0,0) = nu; 442 /* Fdot[0] = F */ 443 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 444 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 445 } else { 446 /* select the current size of the subspace */ 447 if (l < ngmres->msize) l++; 448 k_restart++; 449 /* place the current entry in the list of previous entries */ 450 ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr); 451 ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr); 452 ngmres->fnorms[ivec] = fnorm; 453 if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of r^A */ 454 for (i = 0; i < l; i++) { 455 ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr); 456 Q(i, ivec) = qentry; 457 Q(ivec, i) = qentry; 458 } 459 } 460 461 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 462 snes->iter = k; 463 snes->norm = fnorm; 464 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 465 SNESLogConvHistory(snes, snes->norm, snes->iter); 466 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 467 468 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 469 if (snes->reason) PetscFunctionReturn(0); 470 } 471 snes->reason = SNES_DIVERGED_MAX_IT; 472 PetscFunctionReturn(0); 473 } 474 475 /*MC 476 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 477 478 Level: beginner 479 480 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 481 SIAM Journal on Scientific Computing, 21(5), 2000. 482 483 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 484 J. Assoc. Comput. Mach., 12:547–560, 1965." 485 486 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 487 M*/ 488 489 EXTERN_C_BEGIN 490 #undef __FUNCT__ 491 #define __FUNCT__ "SNESCreate_NGMRES" 492 PetscErrorCode SNESCreate_NGMRES(SNES snes) 493 { 494 SNES_NGMRES *ngmres; 495 PetscErrorCode ierr; 496 497 PetscFunctionBegin; 498 snes->ops->destroy = SNESDestroy_NGMRES; 499 snes->ops->setup = SNESSetUp_NGMRES; 500 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 501 snes->ops->view = SNESView_NGMRES; 502 snes->ops->solve = SNESSolve_NGMRES; 503 snes->ops->reset = SNESReset_NGMRES; 504 505 snes->usespc = PETSC_TRUE; 506 snes->usesksp = PETSC_FALSE; 507 508 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 509 snes->data = (void*) ngmres; 510 ngmres->msize = 10; 511 512 ngmres->restart_it = 2; 513 ngmres->gammaA = 2.0; 514 ngmres->gammaC = 2.0; 515 ngmres->deltaB = 0.9; 516 ngmres->epsilonB = 0.1; 517 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 518 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 519 PetscFunctionReturn(0); 520 } 521 EXTERN_C_END 522