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