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