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 = PetscOptionsBool("-snes_ngmres_anderson", "Use Anderson mixing storage", "SNES", ngmres->anderson, &ngmres->anderson, PETSC_NULL);CHKERRQ(ierr); 89 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 90 ierr = PetscOptionsInt("-snes_ngmres_restart_it","Tolerance iterations before restart","SNES", ngmres->restart_it, &ngmres->restart_it, PETSC_NULL);CHKERRQ(ierr); 91 ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 92 if (debug) { 93 ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 94 } 95 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscOptionsReal("-snes_ngmres_gammaC", " Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 98 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 99 ierr = PetscOptionsTail();CHKERRQ(ierr); 100 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "SNESView_NGMRES" 106 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 107 { 108 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 109 PetscBool iascii; 110 const char *cstr; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 115 if (iascii) { 116 cstr = SNESLineSearchTypeName(snes->ls_type); 117 ierr = PetscViewerASCIIPrintf(viewer, " line search variant: %s\n",cstr);CHKERRQ(ierr); 118 ierr = PetscViewerASCIIPrintf(viewer, " Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer, " Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(viewer, " Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr); 121 } 122 PetscFunctionReturn(0); 123 } 124 125 EXTERN_C_BEGIN 126 #undef __FUNCT__ 127 #define __FUNCT__ "SNESLineSearchSetType_NGMRES" 128 PetscErrorCode SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type) 129 { 130 PetscErrorCode ierr; 131 PetscFunctionBegin; 132 133 switch (type) { 134 case SNES_LS_BASIC: 135 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 136 break; 137 case SNES_LS_BASIC_NONORMS: 138 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 139 break; 140 case SNES_LS_QUADRATIC: 141 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 142 break; 143 default: 144 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type."); 145 break; 146 } 147 snes->ls_type = type; 148 PetscFunctionReturn(0); 149 } 150 EXTERN_C_END 151 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "SNESSolve_NGMRES" 155 156 PetscErrorCode SNESSolve_NGMRES(SNES snes) 157 { 158 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 159 /* present solution, residual, and preconditioned residual */ 160 Vec X, F, B, D, Y; 161 162 /* candidate linear combination answers */ 163 Vec XA, FA, XM, FM, FPC; 164 165 /* previous iterations to construct the subspace */ 166 Vec *Fdot = ngmres->Fdot; 167 Vec *Xdot = ngmres->Xdot; 168 169 /* coefficients and RHS to the minimization problem */ 170 PetscScalar *beta = ngmres->beta; 171 PetscScalar *xi = ngmres->xi; 172 PetscReal fnorm, fMnorm, fAnorm, ynorm, xnorm = 0.0; 173 PetscReal nu; 174 PetscScalar alph_total = 0.; 175 PetscScalar qentry; 176 PetscInt i, j, k, k_restart, l, ivec, restart_count = 0; 177 178 /* solution selection data */ 179 PetscBool selectA, selectRestart; 180 PetscReal dnorm, dminnorm, dcurnorm; 181 PetscReal fminnorm; 182 183 SNESConvergedReason reason; 184 PetscBool lssucceed; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 /* variable initialization */ 189 snes->reason = SNES_CONVERGED_ITERATING; 190 X = snes->vec_sol; 191 F = snes->vec_func; 192 B = snes->vec_rhs; 193 XA = snes->vec_sol_update; 194 FA = snes->work[0]; 195 D = snes->work[1]; 196 197 /* work for the line search */ 198 Y = snes->work[2]; 199 XM = snes->work[3]; 200 FM = snes->work[4]; 201 202 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 203 snes->iter = 0; 204 snes->norm = 0.; 205 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 206 207 /* initialization */ 208 209 /* r = F(x) */ 210 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 211 if (snes->domainerror) { 212 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 213 PetscFunctionReturn(0); 214 } 215 216 /* nu = (r, r) */ 217 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 218 fminnorm = fnorm; 219 nu = fnorm*fnorm; 220 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 221 222 /* q_{00} = nu */ 223 Q(0,0) = nu; 224 ngmres->fnorms[0] = fnorm; 225 /* Fdot[0] = F */ 226 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 227 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 228 229 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 230 snes->norm = fnorm; 231 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 232 SNESLogConvHistory(snes, fnorm, 0); 233 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 234 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 235 if (snes->reason) PetscFunctionReturn(0); 236 237 k_restart = 1; 238 l = 1; 239 for (k=1; k < snes->max_its+1; k++) { 240 /* select which vector of the stored subspace will be updated */ 241 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 242 243 /* Computation of x^M */ 244 if (snes->pc) { 245 ierr = VecCopy(X, XM);CHKERRQ(ierr); 246 ierr = SNESSolve(snes->pc, B, XM);CHKERRQ(ierr); 247 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 248 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 249 snes->reason = SNES_DIVERGED_INNER; 250 PetscFunctionReturn(0); 251 } 252 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 253 ierr = VecCopy(FPC, FM);CHKERRQ(ierr); 254 ierr = SNESGetFunctionNorm(snes->pc, &fMnorm);CHKERRQ(ierr); 255 } else { 256 /* no preconditioner -- just take gradient descent with line search */ 257 ierr = VecCopy(F, Y);CHKERRQ(ierr); 258 ierr = SNESLineSearchPreCheckApply(snes,X,Y,PETSC_NULL); 259 ierr = SNESLineSearchApply(snes,X,F,Y,fnorm,xnorm,XM,FM,&ynorm,&fMnorm,&lssucceed);CHKERRQ(ierr); 260 if (!lssucceed) { 261 if (++snes->numFailures >= snes->maxFailures) { 262 snes->reason = SNES_DIVERGED_LINE_SEARCH; 263 PetscFunctionReturn(0); 264 } 265 } 266 } 267 268 /* r = F(x) */ 269 nu = fMnorm*fMnorm; 270 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */ 271 272 /* construct the right hand side and xi factors */ 273 ierr = VecMDot(FM, l, Fdot, xi);CHKERRQ(ierr); 274 275 for (i = 0; i < l; i++) { 276 beta[i] = nu - xi[i]; 277 } 278 279 /* construct h */ 280 for (j = 0; j < l; j++) { 281 for (i = 0; i < l; i++) { 282 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 283 } 284 } 285 286 if(l == 1) { 287 /* simply set alpha[0] = beta[0] / H[0, 0] */ 288 beta[0] = beta[0] / H(0, 0); 289 } else { 290 #ifdef PETSC_MISSING_LAPACK_GELSS 291 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 292 #else 293 ngmres->m = PetscBLASIntCast(l); 294 ngmres->n = PetscBLASIntCast(l); 295 ngmres->info = PetscBLASIntCast(0); 296 ngmres->rcond = -1.; 297 #ifdef PETSC_USE_COMPLEX 298 LAPACKgelss_(&ngmres->m, 299 &ngmres->n, 300 &ngmres->nrhs, 301 ngmres->h, 302 &ngmres->lda, 303 ngmres->beta, 304 &ngmres->ldb, 305 ngmres->s, 306 &ngmres->rcond, 307 &ngmres->rank, 308 ngmres->work, 309 &ngmres->lwork, 310 ngmres->rwork, 311 &ngmres->info); 312 #else 313 LAPACKgelss_(&ngmres->m, 314 &ngmres->n, 315 &ngmres->nrhs, 316 ngmres->h, 317 &ngmres->lda, 318 ngmres->beta, 319 &ngmres->ldb, 320 ngmres->s, 321 &ngmres->rcond, 322 &ngmres->rank, 323 ngmres->work, 324 &ngmres->lwork, 325 &ngmres->info); 326 #endif 327 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 328 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 329 #endif 330 } 331 332 alph_total = 0.; 333 for (i = 0; i < l; i++) { 334 alph_total += beta[i]; 335 } 336 337 ierr = VecCopy(XM, XA);CHKERRQ(ierr); 338 ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr); 339 340 ierr = VecMAXPY(XA, l, beta, Xdot);CHKERRQ(ierr); 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 if (ngmres->anderson) { 442 ngmres->fnorms[0] = fMnorm; 443 nu = fMnorm*fMnorm; 444 Q(0,0) = nu; 445 /* Fdot[0] = F */ 446 ierr = VecCopy(XM, Xdot[0]);CHKERRQ(ierr); 447 ierr = VecCopy(FM, Fdot[0]);CHKERRQ(ierr); 448 } else { 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 } 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 if (ngmres->anderson) { 462 ierr = VecCopy(FM, Fdot[ivec]);CHKERRQ(ierr); 463 ierr = VecCopy(XM, Xdot[ivec]);CHKERRQ(ierr); 464 ngmres->fnorms[ivec] = fMnorm; 465 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of FM */ 466 for (i = 0; i < l; i++) { 467 ierr = VecDot(FM, Fdot[i], &qentry);CHKERRQ(ierr); 468 Q(i, ivec) = qentry; 469 Q(ivec, i) = qentry; 470 } 471 } else { 472 ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr); 473 ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr); 474 ngmres->fnorms[ivec] = fnorm; 475 if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of FA */ 476 for (i = 0; i < l; i++) { 477 ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr); 478 Q(i, ivec) = qentry; 479 Q(ivec, i) = qentry; 480 } 481 } 482 } 483 484 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 485 snes->iter = k; 486 snes->norm = fnorm; 487 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 488 SNESLogConvHistory(snes, snes->norm, snes->iter); 489 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 490 491 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 492 if (snes->reason) PetscFunctionReturn(0); 493 } 494 snes->reason = SNES_DIVERGED_MAX_IT; 495 PetscFunctionReturn(0); 496 } 497 498 /*MC 499 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 500 501 Level: beginner 502 503 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 504 SIAM Journal on Scientific Computing, 21(5), 2000. 505 506 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 507 J. Assoc. Comput. Mach., 12:547–560, 1965." 508 509 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 510 M*/ 511 512 EXTERN_C_BEGIN 513 #undef __FUNCT__ 514 #define __FUNCT__ "SNESCreate_NGMRES" 515 PetscErrorCode SNESCreate_NGMRES(SNES snes) 516 { 517 SNES_NGMRES *ngmres; 518 PetscErrorCode ierr; 519 520 PetscFunctionBegin; 521 snes->ops->destroy = SNESDestroy_NGMRES; 522 snes->ops->setup = SNESSetUp_NGMRES; 523 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 524 snes->ops->view = SNESView_NGMRES; 525 snes->ops->solve = SNESSolve_NGMRES; 526 snes->ops->reset = SNESReset_NGMRES; 527 528 snes->usespc = PETSC_TRUE; 529 snes->usesksp = PETSC_FALSE; 530 531 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 532 snes->data = (void*) ngmres; 533 ngmres->msize = 30; 534 535 snes->max_funcs = 30000; 536 snes->max_its = 10000; 537 538 ngmres->additive = PETSC_FALSE; 539 ngmres->anderson = PETSC_FALSE; 540 541 ngmres->restart_it = 2; 542 ngmres->gammaA = 2.0; 543 ngmres->gammaC = 2.0; 544 ngmres->deltaB = 0.9; 545 ngmres->epsilonB = 0.1; 546 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 547 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 548 PetscFunctionReturn(0); 549 } 550 EXTERN_C_END 551