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; 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 = SNESComputeFunction(snes, XM, FM);CHKERRQ(ierr); 252 ierr = VecNorm(FM, NORM_2, &fMnorm);CHKERRQ(ierr); 253 } else { 254 /* no preconditioner -- just take gradient descent with line search */ 255 ierr = VecCopy(F, Y);CHKERRQ(ierr); 256 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FM,XM,&ynorm,&fMnorm,&lssucceed);CHKERRQ(ierr); 257 if (!lssucceed) { 258 if (++snes->numFailures >= snes->maxFailures) { 259 snes->reason = SNES_DIVERGED_LINE_SEARCH; 260 PetscFunctionReturn(0); 261 } 262 } 263 } 264 265 /* r = F(x) */ 266 nu = fMnorm*fMnorm; 267 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */ 268 269 /* construct the right hand side and xi factors */ 270 ierr = VecMDot(FM, l, Fdot, xi);CHKERRQ(ierr); 271 272 for (i = 0; i < l; i++) { 273 beta[i] = nu - xi[i]; 274 } 275 276 /* construct h */ 277 for (j = 0; j < l; j++) { 278 for (i = 0; i < l; i++) { 279 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 280 } 281 } 282 283 if(l == 1) { 284 /* simply set alpha[0] = beta[0] / H[0, 0] */ 285 beta[0] = beta[0] / H(0, 0); 286 } else { 287 #ifdef PETSC_MISSING_LAPACK_GELSS 288 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 289 #else 290 ngmres->m = PetscBLASIntCast(l); 291 ngmres->n = PetscBLASIntCast(l); 292 ngmres->info = PetscBLASIntCast(0); 293 ngmres->rcond = -1.; 294 #ifdef PETSC_USE_COMPLEX 295 LAPACKgelss_(&ngmres->m, 296 &ngmres->n, 297 &ngmres->nrhs, 298 ngmres->h, 299 &ngmres->lda, 300 ngmres->beta, 301 &ngmres->ldb, 302 ngmres->s, 303 &ngmres->rcond, 304 &ngmres->rank, 305 ngmres->work, 306 &ngmres->lwork, 307 ngmres->rwork, 308 &ngmres->info); 309 #else 310 LAPACKgelss_(&ngmres->m, 311 &ngmres->n, 312 &ngmres->nrhs, 313 ngmres->h, 314 &ngmres->lda, 315 ngmres->beta, 316 &ngmres->ldb, 317 ngmres->s, 318 &ngmres->rcond, 319 &ngmres->rank, 320 ngmres->work, 321 &ngmres->lwork, 322 &ngmres->info); 323 #endif 324 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 325 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 326 #endif 327 } 328 329 alph_total = 0.; 330 for (i = 0; i < l; i++) { 331 alph_total += beta[i]; 332 } 333 334 ierr = VecCopy(XM, XA);CHKERRQ(ierr); 335 ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr); 336 337 ierr = VecMAXPY(XA, l, beta, Xdot);CHKERRQ(ierr); 338 ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr); 339 ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr); 340 341 /* differences for selection and restart */ 342 ierr=VecCopy(XA, D);CHKERRQ(ierr); 343 ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr); 344 ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr); 345 dminnorm = -1.0; 346 for(i=0;i<l;i++) { 347 ierr=VecCopy(XA, D);CHKERRQ(ierr); 348 ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr); 349 ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr); 350 if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm; 351 } 352 353 /* combination (additive) or selection (multiplicative) of the N-GMRES solution */ 354 if (ngmres->additive) { 355 /* X = X + \lambda(XA - X) */ 356 if (ngmres->monitor) { 357 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);CHKERRQ(ierr); 358 } 359 ierr = VecCopy(XA, Y);CHKERRQ(ierr); 360 ierr = VecAYPX(Y, -1.0, X);CHKERRQ(ierr); 361 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FA,XA,&ynorm,&fnorm,&lssucceed);CHKERRQ(ierr); 362 if (!lssucceed) { 363 if (++snes->numFailures >= snes->maxFailures) { 364 snes->reason = SNES_DIVERGED_LINE_SEARCH; 365 PetscFunctionReturn(0); 366 } 367 } 368 if (ngmres->monitor) { 369 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", fnorm);CHKERRQ(ierr); 370 } 371 ierr = VecCopy(FA, F);CHKERRQ(ierr); 372 ierr = VecCopy(XA, X);CHKERRQ(ierr); 373 } else { 374 selectA = PETSC_TRUE; 375 /* Conditions for choosing the accelerated answer */ 376 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 377 if (fAnorm >= ngmres->gammaA*fminnorm) { 378 selectA = PETSC_FALSE; 379 } 380 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 381 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 382 } else { 383 selectA=PETSC_FALSE; 384 } 385 if (selectA) { 386 if (ngmres->monitor) { 387 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 388 } 389 /* copy it over */ 390 fnorm = fAnorm; 391 nu = fnorm*fnorm; 392 ierr = VecCopy(FA, F);CHKERRQ(ierr); 393 ierr = VecCopy(XA, X);CHKERRQ(ierr); 394 } else { 395 if (ngmres->monitor) { 396 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 397 } 398 fnorm = fMnorm; 399 nu = fnorm*fnorm; 400 ierr = VecCopy(FM, F);CHKERRQ(ierr); 401 ierr = VecCopy(XM, X);CHKERRQ(ierr); 402 } 403 } 404 405 selectRestart = PETSC_FALSE; 406 407 /* difference stagnation restart */ 408 if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) { 409 if (ngmres->monitor) { 410 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr); 411 } 412 selectRestart = PETSC_TRUE; 413 } 414 /* residual stagnation restart */ 415 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 416 if (ngmres->monitor) { 417 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 418 } 419 selectRestart = PETSC_TRUE; 420 } 421 422 /* if the restart conditions persist for more than restart_it iterations, restart. */ 423 if (selectRestart) { 424 restart_count++; 425 } else { 426 restart_count = 0; 427 } 428 429 /* restart after restart conditions have persisted for a fixed number of iterations */ 430 if (restart_count >= ngmres->restart_it) { 431 if (ngmres->monitor){ 432 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 433 } 434 restart_count = 0; 435 k_restart = 1; 436 l = 1; 437 /* q_{00} = nu */ 438 ngmres->fnorms[0] = fnorm; 439 nu = fnorm*fnorm; 440 Q(0,0) = nu; 441 /* Fdot[0] = F */ 442 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 443 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 444 } else { 445 /* select the current size of the subspace */ 446 if (l < ngmres->msize) l++; 447 k_restart++; 448 /* place the current entry in the list of previous entries */ 449 ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr); 450 ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr); 451 ngmres->fnorms[ivec] = fnorm; 452 if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of r^A */ 453 for (i = 0; i < l; i++) { 454 ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr); 455 Q(i, ivec) = qentry; 456 Q(ivec, i) = qentry; 457 } 458 } 459 460 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 461 snes->iter = k; 462 snes->norm = fnorm; 463 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 464 SNESLogConvHistory(snes, snes->norm, snes->iter); 465 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 466 467 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 468 if (snes->reason) PetscFunctionReturn(0); 469 } 470 snes->reason = SNES_DIVERGED_MAX_IT; 471 PetscFunctionReturn(0); 472 } 473 474 /*MC 475 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 476 477 Level: beginner 478 479 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 480 SIAM Journal on Scientific Computing, 21(5), 2000. 481 482 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 483 J. Assoc. Comput. Mach., 12:547–560, 1965." 484 485 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 486 M*/ 487 488 EXTERN_C_BEGIN 489 #undef __FUNCT__ 490 #define __FUNCT__ "SNESCreate_NGMRES" 491 PetscErrorCode SNESCreate_NGMRES(SNES snes) 492 { 493 SNES_NGMRES *ngmres; 494 PetscErrorCode ierr; 495 496 PetscFunctionBegin; 497 snes->ops->destroy = SNESDestroy_NGMRES; 498 snes->ops->setup = SNESSetUp_NGMRES; 499 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 500 snes->ops->view = SNESView_NGMRES; 501 snes->ops->solve = SNESSolve_NGMRES; 502 snes->ops->reset = SNESReset_NGMRES; 503 504 snes->usespc = PETSC_TRUE; 505 snes->usesksp = PETSC_FALSE; 506 507 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 508 snes->data = (void*) ngmres; 509 ngmres->msize = 10; 510 511 ngmres->restart_it = 2; 512 ngmres->gammaA = 2.0; 513 ngmres->gammaC = 2.0; 514 ngmres->deltaB = 0.9; 515 ngmres->epsilonB = 0.1; 516 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 517 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 EXTERN_C_END 521