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