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 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 16 PetscFunctionReturn(0); 17 } 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "SNESDestroy_NGMRES" 21 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 22 { 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 27 if (snes->data) { 28 SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data; 29 ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->fnorms, ngmres->q);CHKERRQ(ierr); 30 ierr = PetscFree(ngmres->s);CHKERRQ(ierr); 31 #if PETSC_USE_COMPLEX 32 ierr = PetscFree(ngmres->rwork); 33 #endif 34 ierr = PetscFree(ngmres->work); 35 } 36 ierr = PetscFree(snes->data); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "SNESSetUp_NGMRES" 42 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 43 { 44 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 45 PetscInt msize,hsize; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 msize = ngmres->msize; /* restart size */ 50 hsize = msize * msize; 51 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 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr); 74 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr); 75 ierr = SNESDefaultGetWork(snes, 5);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 81 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 82 { 83 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 84 PetscErrorCode ierr; 85 PetscBool debug; 86 87 PetscFunctionBegin; 88 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 89 ierr = PetscOptionsBool("-snes_ngmres_additive", "Use additive variant vs. choice", "SNES", ngmres->additive, &ngmres->additive, PETSC_NULL);CHKERRQ(ierr); 90 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 91 ierr = PetscOptionsInt("-snes_ngmres_restart_it","Tolerance iterations before restart","SNES", ngmres->restart_it, &ngmres->restart_it, PETSC_NULL);CHKERRQ(ierr); 92 ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 93 if (debug) { 94 ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 95 } 96 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsReal("-snes_ngmres_gammaC", " Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 98 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 99 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 100 ierr = PetscOptionsTail();CHKERRQ(ierr); 101 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 102 PetscFunctionReturn(0); 103 } 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "SNESView_NGMRES" 107 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 108 { 109 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 110 PetscBool iascii; 111 const char *cstr; 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 116 if (iascii) { 117 cstr = SNESLineSearchTypeName(snes->ls_type); 118 ierr = PetscViewerASCIIPrintf(viewer, " line search variant: %s\n",cstr);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer, " Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(viewer, " Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPrintf(viewer, " Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr); 122 } 123 PetscFunctionReturn(0); 124 } 125 126 EXTERN_C_BEGIN 127 #undef __FUNCT__ 128 #define __FUNCT__ "SNESLineSearchSetType_NGMRES" 129 PetscErrorCode SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type) 130 { 131 PetscErrorCode ierr; 132 PetscFunctionBegin; 133 134 switch (type) { 135 case SNES_LS_BASIC: 136 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 137 break; 138 case SNES_LS_BASIC_NONORMS: 139 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 140 break; 141 case SNES_LS_QUADRATIC: 142 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 143 break; 144 default: 145 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type."); 146 break; 147 } 148 snes->ls_type = type; 149 PetscFunctionReturn(0); 150 } 151 EXTERN_C_END 152 153 154 #undef __FUNCT__ 155 #define __FUNCT__ "SNESSolve_NGMRES" 156 157 PetscErrorCode SNESSolve_NGMRES(SNES snes) 158 { 159 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 160 /* present solution, residual, and preconditioned residual */ 161 Vec X, F, B, D, Y; 162 163 /* candidate linear combination answers */ 164 Vec XA, FA, XM, FM; 165 166 /* previous iterations to construct the subspace */ 167 Vec *Fdot = ngmres->Fdot; 168 Vec *Xdot = ngmres->Xdot; 169 170 /* coefficients and RHS to the minimization problem */ 171 PetscScalar *beta = ngmres->beta; 172 PetscScalar *xi = ngmres->xi; 173 PetscReal fnorm, fMnorm, fAnorm, ynorm, xnorm = 0.0; 174 PetscReal nu; 175 PetscScalar alph_total = 0.; 176 PetscScalar qentry; 177 PetscInt i, j, k, k_restart, l, ivec, restart_count = 0; 178 179 /* solution selection data */ 180 PetscBool selectA, selectRestart; 181 PetscReal dnorm, dminnorm, dcurnorm; 182 PetscReal fminnorm; 183 184 SNESConvergedReason reason; 185 PetscBool lssucceed; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 /* variable initialization */ 190 snes->reason = SNES_CONVERGED_ITERATING; 191 X = snes->vec_sol; 192 F = snes->vec_func; 193 B = snes->vec_rhs; 194 XA = snes->vec_sol_update; 195 FA = snes->work[0]; 196 D = snes->work[1]; 197 198 /* work for the line search */ 199 Y = snes->work[2]; 200 XM = snes->work[3]; 201 FM = snes->work[4]; 202 203 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 204 snes->iter = 0; 205 snes->norm = 0.; 206 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 207 208 /* initialization */ 209 210 /* r = F(x) */ 211 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 212 if (snes->domainerror) { 213 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 214 PetscFunctionReturn(0); 215 } 216 217 /* nu = (r, r) */ 218 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 219 fminnorm = fnorm; 220 nu = fnorm*fnorm; 221 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 222 223 /* q_{00} = nu */ 224 Q(0,0) = nu; 225 ngmres->fnorms[0] = fnorm; 226 /* Fdot[0] = F */ 227 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 228 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 229 230 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 231 snes->norm = fnorm; 232 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 233 SNESLogConvHistory(snes, fnorm, 0); 234 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 235 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 236 if (snes->reason) PetscFunctionReturn(0); 237 238 k_restart = 1; 239 l = 1; 240 for (k=1; k < snes->max_its+1; k++) { 241 /* select which vector of the stored subspace will be updated */ 242 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 243 244 /* Computation of x^M */ 245 if (snes->pc) { 246 ierr = VecCopy(X, XM);CHKERRQ(ierr); 247 ierr = SNESSolve(snes->pc, B, XM);CHKERRQ(ierr); 248 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 249 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 250 snes->reason = SNES_DIVERGED_INNER; 251 PetscFunctionReturn(0); 252 } 253 ierr = SNESComputeFunction(snes, XM, FM);CHKERRQ(ierr); 254 ierr = VecNorm(FM, NORM_2, &fMnorm);CHKERRQ(ierr); 255 } else { 256 /* no preconditioner -- just take gradient descent with line search */ 257 ierr = VecCopy(F, Y);CHKERRQ(ierr); 258 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FM,XM,&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 ngmres->restart_it = 2; 514 ngmres->gammaA = 2.0; 515 ngmres->gammaC = 2.0; 516 ngmres->deltaB = 0.9; 517 ngmres->epsilonB = 0.1; 518 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 519 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 EXTERN_C_END 523