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 = (*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