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 for (i = 0; i < l; i++) { 273 ierr = VecDot(Fdot[i], FM, &xi[i]);CHKERRQ(ierr); 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 for(i=0;i<l;i++){ 339 ierr= VecAXPY(XA, beta[i], Xdot[i]);CHKERRQ(ierr); 340 } 341 ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr); 342 ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr); 343 344 /* differences for selection and restart */ 345 ierr=VecCopy(XA, D);CHKERRQ(ierr); 346 ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr); 347 ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr); 348 dminnorm = -1.0; 349 for(i=0;i<l;i++) { 350 ierr=VecCopy(XA, D);CHKERRQ(ierr); 351 ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr); 352 ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr); 353 if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm; 354 } 355 356 /* combination (additive) or selection (multiplicative) of the N-GMRES solution */ 357 if (ngmres->additive) { 358 /* X = X + \lambda(XA - X) */ 359 if (ngmres->monitor) { 360 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);CHKERRQ(ierr); 361 } 362 ierr = VecCopy(XA, Y);CHKERRQ(ierr); 363 ierr = VecAYPX(Y, -1.0, X);CHKERRQ(ierr); 364 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FA,XA,&ynorm,&fnorm,&lssucceed);CHKERRQ(ierr); 365 if (!lssucceed) { 366 if (++snes->numFailures >= snes->maxFailures) { 367 snes->reason = SNES_DIVERGED_LINE_SEARCH; 368 PetscFunctionReturn(0); 369 } 370 } 371 if (ngmres->monitor) { 372 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", fnorm);CHKERRQ(ierr); 373 } 374 ierr = VecCopy(FA, F);CHKERRQ(ierr); 375 ierr = VecCopy(XA, X);CHKERRQ(ierr); 376 } else { 377 selectA = PETSC_TRUE; 378 /* Conditions for choosing the accelerated answer */ 379 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 380 if (fAnorm >= ngmres->gammaA*fminnorm) { 381 selectA = PETSC_FALSE; 382 } 383 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 384 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 385 } else { 386 selectA=PETSC_FALSE; 387 } 388 if (selectA) { 389 if (ngmres->monitor) { 390 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 391 } 392 /* copy it over */ 393 fnorm = fAnorm; 394 nu = fnorm*fnorm; 395 ierr = VecCopy(FA, F);CHKERRQ(ierr); 396 ierr = VecCopy(XA, X);CHKERRQ(ierr); 397 } else { 398 if (ngmres->monitor) { 399 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 400 } 401 fnorm = fMnorm; 402 nu = fnorm*fnorm; 403 ierr = VecCopy(FM, F);CHKERRQ(ierr); 404 ierr = VecCopy(XM, X);CHKERRQ(ierr); 405 } 406 } 407 408 selectRestart = PETSC_FALSE; 409 410 /* difference stagnation restart */ 411 if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) { 412 if (ngmres->monitor) { 413 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr); 414 } 415 selectRestart = PETSC_TRUE; 416 } 417 /* residual stagnation restart */ 418 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 419 if (ngmres->monitor) { 420 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 421 } 422 selectRestart = PETSC_TRUE; 423 } 424 425 /* if the restart conditions persist for more than restart_it iterations, restart. */ 426 if (selectRestart) { 427 restart_count++; 428 } else { 429 restart_count = 0; 430 } 431 432 /* restart after restart conditions have persisted for a fixed number of iterations */ 433 if (restart_count >= ngmres->restart_it) { 434 if (ngmres->monitor){ 435 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 436 } 437 restart_count = 0; 438 k_restart = 1; 439 l = 1; 440 /* q_{00} = nu */ 441 ngmres->fnorms[0] = fnorm; 442 nu = fnorm*fnorm; 443 Q(0,0) = nu; 444 /* Fdot[0] = F */ 445 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 446 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 447 } else { 448 /* select the current size of the subspace */ 449 if (l < ngmres->msize) l++; 450 k_restart++; 451 /* place the current entry in the list of previous entries */ 452 ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr); 453 ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr); 454 ngmres->fnorms[ivec] = fnorm; 455 if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of r^A */ 456 for (i = 0; i < l; i++) { 457 ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr); 458 Q(i, ivec) = qentry; 459 Q(ivec, i) = qentry; 460 } 461 } 462 463 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 464 snes->iter = k; 465 snes->norm = fnorm; 466 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 467 SNESLogConvHistory(snes, snes->norm, snes->iter); 468 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 469 470 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 471 if (snes->reason) PetscFunctionReturn(0); 472 } 473 snes->reason = SNES_DIVERGED_MAX_IT; 474 PetscFunctionReturn(0); 475 } 476 477 /*MC 478 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 479 480 Level: beginner 481 482 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 483 SIAM Journal on Scientific Computing, 21(5), 2000. 484 485 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 486 J. Assoc. Comput. Mach., 12:547–560, 1965." 487 488 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 489 M*/ 490 491 EXTERN_C_BEGIN 492 #undef __FUNCT__ 493 #define __FUNCT__ "SNESCreate_NGMRES" 494 PetscErrorCode SNESCreate_NGMRES(SNES snes) 495 { 496 SNES_NGMRES *ngmres; 497 PetscErrorCode ierr; 498 499 PetscFunctionBegin; 500 snes->ops->destroy = SNESDestroy_NGMRES; 501 snes->ops->setup = SNESSetUp_NGMRES; 502 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 503 snes->ops->view = SNESView_NGMRES; 504 snes->ops->solve = SNESSolve_NGMRES; 505 snes->ops->reset = SNESReset_NGMRES; 506 507 snes->usespc = PETSC_TRUE; 508 snes->usesksp = PETSC_FALSE; 509 510 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 511 snes->data = (void*) ngmres; 512 ngmres->msize = 10; 513 514 ngmres->restart_it = 2; 515 ngmres->gammaA = 2.0; 516 ngmres->gammaC = 2.0; 517 ngmres->deltaB = 0.9; 518 ngmres->epsilonB = 0.1; 519 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 520 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 EXTERN_C_END 524