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 case SNES_LS_SECANT: 144 ierr = SNESLineSearchSet(snes,SNESLineSearchSecant,PETSC_NULL);CHKERRQ(ierr); 145 break; 146 default: 147 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type."); 148 break; 149 } 150 snes->ls_type = type; 151 PetscFunctionReturn(0); 152 } 153 EXTERN_C_END 154 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "SNESSolve_NGMRES" 158 159 PetscErrorCode SNESSolve_NGMRES(SNES snes) 160 { 161 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 162 /* present solution, residual, and preconditioned residual */ 163 Vec X, F, B, D, Y; 164 165 /* candidate linear combination answers */ 166 Vec XA, FA, XM, FM, FPC; 167 168 /* previous iterations to construct the subspace */ 169 Vec *Fdot = ngmres->Fdot; 170 Vec *Xdot = ngmres->Xdot; 171 172 /* coefficients and RHS to the minimization problem */ 173 PetscScalar *beta = ngmres->beta; 174 PetscScalar *xi = ngmres->xi; 175 PetscReal fnorm, fMnorm, fAnorm, ynorm, xnorm = 0.0; 176 PetscReal nu; 177 PetscScalar alph_total = 0.; 178 PetscScalar qentry; 179 PetscInt i, j, k, k_restart, l, ivec, restart_count = 0; 180 181 /* solution selection data */ 182 PetscBool selectA, selectRestart; 183 PetscReal dnorm, dminnorm, dcurnorm; 184 PetscReal fminnorm; 185 186 SNESConvergedReason reason; 187 PetscBool lssucceed; 188 PetscErrorCode ierr; 189 190 PetscFunctionBegin; 191 /* variable initialization */ 192 snes->reason = SNES_CONVERGED_ITERATING; 193 X = snes->vec_sol; 194 F = snes->vec_func; 195 B = snes->vec_rhs; 196 XA = snes->vec_sol_update; 197 FA = snes->work[0]; 198 D = snes->work[1]; 199 200 /* work for the line search */ 201 Y = snes->work[2]; 202 XM = snes->work[3]; 203 FM = snes->work[4]; 204 205 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 206 snes->iter = 0; 207 snes->norm = 0.; 208 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 209 210 /* initialization */ 211 212 /* r = F(x) */ 213 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 214 if (snes->domainerror) { 215 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 216 PetscFunctionReturn(0); 217 } 218 219 /* nu = (r, r) */ 220 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 221 fminnorm = fnorm; 222 nu = fnorm*fnorm; 223 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 224 225 /* q_{00} = nu */ 226 Q(0,0) = nu; 227 ngmres->fnorms[0] = fnorm; 228 /* Fdot[0] = F */ 229 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 230 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 231 232 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 233 snes->norm = fnorm; 234 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 235 SNESLogConvHistory(snes, fnorm, 0); 236 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 237 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 238 if (snes->reason) PetscFunctionReturn(0); 239 240 k_restart = 1; 241 l = 1; 242 for (k=1; k < snes->max_its+1; k++) { 243 /* select which vector of the stored subspace will be updated */ 244 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 245 246 /* Computation of x^M */ 247 if (snes->pc) { 248 ierr = VecCopy(X, XM);CHKERRQ(ierr); 249 ierr = SNESSolve(snes->pc, B, XM);CHKERRQ(ierr); 250 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 251 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 252 snes->reason = SNES_DIVERGED_INNER; 253 PetscFunctionReturn(0); 254 } 255 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 256 ierr = VecCopy(FPC, FM);CHKERRQ(ierr); 257 ierr = SNESGetFunctionNorm(snes->pc, &fMnorm);CHKERRQ(ierr); 258 } else { 259 /* no preconditioner -- just take gradient descent with line search */ 260 ierr = VecCopy(F, Y);CHKERRQ(ierr); 261 ierr = SNESLineSearchPreCheckApply(snes,X,Y,PETSC_NULL); 262 ierr = SNESLineSearchApply(snes,X,F,Y,fnorm,xnorm,XM,FM,&ynorm,&fMnorm,&lssucceed);CHKERRQ(ierr); 263 if (!lssucceed) { 264 if (++snes->numFailures >= snes->maxFailures) { 265 snes->reason = SNES_DIVERGED_LINE_SEARCH; 266 PetscFunctionReturn(0); 267 } 268 } 269 } 270 271 /* r = F(x) */ 272 nu = fMnorm*fMnorm; 273 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */ 274 275 /* construct the right hand side and xi factors */ 276 ierr = VecMDot(FM, l, Fdot, xi);CHKERRQ(ierr); 277 278 for (i = 0; i < l; i++) { 279 beta[i] = nu - xi[i]; 280 } 281 282 /* construct h */ 283 for (j = 0; j < l; j++) { 284 for (i = 0; i < l; i++) { 285 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 286 } 287 } 288 289 if(l == 1) { 290 /* simply set alpha[0] = beta[0] / H[0, 0] */ 291 beta[0] = beta[0] / H(0, 0); 292 } else { 293 #ifdef PETSC_MISSING_LAPACK_GELSS 294 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 295 #else 296 ngmres->m = PetscBLASIntCast(l); 297 ngmres->n = PetscBLASIntCast(l); 298 ngmres->info = PetscBLASIntCast(0); 299 ngmres->rcond = -1.; 300 #ifdef PETSC_USE_COMPLEX 301 LAPACKgelss_(&ngmres->m, 302 &ngmres->n, 303 &ngmres->nrhs, 304 ngmres->h, 305 &ngmres->lda, 306 ngmres->beta, 307 &ngmres->ldb, 308 ngmres->s, 309 &ngmres->rcond, 310 &ngmres->rank, 311 ngmres->work, 312 &ngmres->lwork, 313 ngmres->rwork, 314 &ngmres->info); 315 #else 316 LAPACKgelss_(&ngmres->m, 317 &ngmres->n, 318 &ngmres->nrhs, 319 ngmres->h, 320 &ngmres->lda, 321 ngmres->beta, 322 &ngmres->ldb, 323 ngmres->s, 324 &ngmres->rcond, 325 &ngmres->rank, 326 ngmres->work, 327 &ngmres->lwork, 328 &ngmres->info); 329 #endif 330 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 331 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 332 #endif 333 } 334 335 alph_total = 0.; 336 for (i = 0; i < l; i++) { 337 alph_total += beta[i]; 338 } 339 340 ierr = VecCopy(XM, XA);CHKERRQ(ierr); 341 ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr); 342 343 ierr = VecMAXPY(XA, l, beta, Xdot);CHKERRQ(ierr); 344 ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr); 345 ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr); 346 347 /* differences for selection and restart */ 348 ierr=VecCopy(XA, D);CHKERRQ(ierr); 349 ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr); 350 ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr); 351 dminnorm = -1.0; 352 for(i=0;i<l;i++) { 353 ierr=VecCopy(XA, D);CHKERRQ(ierr); 354 ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr); 355 ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr); 356 if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm; 357 } 358 359 /* combination (additive) or selection (multiplicative) of the N-GMRES solution */ 360 if (ngmres->additive) { 361 /* X = X + \lambda(XA - X) */ 362 if (ngmres->monitor) { 363 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);CHKERRQ(ierr); 364 } 365 ierr = VecCopy(XA, Y);CHKERRQ(ierr); 366 ierr = VecAYPX(Y, -1.0, X);CHKERRQ(ierr); 367 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,FA,XA,&ynorm,&fnorm,&lssucceed);CHKERRQ(ierr); 368 if (!lssucceed) { 369 if (++snes->numFailures >= snes->maxFailures) { 370 snes->reason = SNES_DIVERGED_LINE_SEARCH; 371 PetscFunctionReturn(0); 372 } 373 } 374 if (ngmres->monitor) { 375 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", fnorm);CHKERRQ(ierr); 376 } 377 ierr = VecCopy(FA, F);CHKERRQ(ierr); 378 ierr = VecCopy(XA, X);CHKERRQ(ierr); 379 } else { 380 selectA = PETSC_TRUE; 381 /* Conditions for choosing the accelerated answer */ 382 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 383 if (fAnorm >= ngmres->gammaA*fminnorm) { 384 selectA = PETSC_FALSE; 385 } 386 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 387 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 388 } else { 389 selectA=PETSC_FALSE; 390 } 391 if (selectA) { 392 if (ngmres->monitor) { 393 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 394 } 395 /* copy it over */ 396 fnorm = fAnorm; 397 nu = fnorm*fnorm; 398 ierr = VecCopy(FA, F);CHKERRQ(ierr); 399 ierr = VecCopy(XA, X);CHKERRQ(ierr); 400 } else { 401 if (ngmres->monitor) { 402 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr); 403 } 404 fnorm = fMnorm; 405 nu = fnorm*fnorm; 406 ierr = VecCopy(FM, F);CHKERRQ(ierr); 407 ierr = VecCopy(XM, X);CHKERRQ(ierr); 408 } 409 } 410 411 selectRestart = PETSC_FALSE; 412 413 /* difference stagnation restart */ 414 if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) { 415 if (ngmres->monitor) { 416 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr); 417 } 418 selectRestart = PETSC_TRUE; 419 } 420 /* residual stagnation restart */ 421 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 422 if (ngmres->monitor) { 423 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 424 } 425 selectRestart = PETSC_TRUE; 426 } 427 428 /* if the restart conditions persist for more than restart_it iterations, restart. */ 429 if (selectRestart) { 430 restart_count++; 431 } else { 432 restart_count = 0; 433 } 434 435 /* restart after restart conditions have persisted for a fixed number of iterations */ 436 if (restart_count >= ngmres->restart_it) { 437 if (ngmres->monitor){ 438 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 439 } 440 restart_count = 0; 441 k_restart = 1; 442 l = 1; 443 /* q_{00} = nu */ 444 if (ngmres->anderson) { 445 ngmres->fnorms[0] = fMnorm; 446 nu = fMnorm*fMnorm; 447 Q(0,0) = nu; 448 /* Fdot[0] = F */ 449 ierr = VecCopy(XM, Xdot[0]);CHKERRQ(ierr); 450 ierr = VecCopy(FM, Fdot[0]);CHKERRQ(ierr); 451 } else { 452 ngmres->fnorms[0] = fnorm; 453 nu = fnorm*fnorm; 454 Q(0,0) = nu; 455 /* Fdot[0] = F */ 456 ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr); 457 ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr); 458 } 459 } else { 460 /* select the current size of the subspace */ 461 if (l < ngmres->msize) l++; 462 k_restart++; 463 /* place the current entry in the list of previous entries */ 464 if (ngmres->anderson) { 465 ierr = VecCopy(FM, Fdot[ivec]);CHKERRQ(ierr); 466 ierr = VecCopy(XM, Xdot[ivec]);CHKERRQ(ierr); 467 ngmres->fnorms[ivec] = fMnorm; 468 if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of FM */ 469 for (i = 0; i < l; i++) { 470 ierr = VecDot(FM, Fdot[i], &qentry);CHKERRQ(ierr); 471 Q(i, ivec) = qentry; 472 Q(ivec, i) = qentry; 473 } 474 } else { 475 ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr); 476 ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr); 477 ngmres->fnorms[ivec] = fnorm; 478 if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of FA */ 479 for (i = 0; i < l; i++) { 480 ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr); 481 Q(i, ivec) = qentry; 482 Q(ivec, i) = qentry; 483 } 484 } 485 } 486 487 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 488 snes->iter = k; 489 snes->norm = fnorm; 490 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 491 SNESLogConvHistory(snes, snes->norm, snes->iter); 492 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 493 494 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 495 if (snes->reason) PetscFunctionReturn(0); 496 } 497 snes->reason = SNES_DIVERGED_MAX_IT; 498 PetscFunctionReturn(0); 499 } 500 501 /*MC 502 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 503 504 Level: beginner 505 506 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 507 SIAM Journal on Scientific Computing, 21(5), 2000. 508 509 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 510 J. Assoc. Comput. Mach., 12:547–560, 1965." 511 512 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 513 M*/ 514 515 EXTERN_C_BEGIN 516 #undef __FUNCT__ 517 #define __FUNCT__ "SNESCreate_NGMRES" 518 PetscErrorCode SNESCreate_NGMRES(SNES snes) 519 { 520 SNES_NGMRES *ngmres; 521 PetscErrorCode ierr; 522 523 PetscFunctionBegin; 524 snes->ops->destroy = SNESDestroy_NGMRES; 525 snes->ops->setup = SNESSetUp_NGMRES; 526 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 527 snes->ops->view = SNESView_NGMRES; 528 snes->ops->solve = SNESSolve_NGMRES; 529 snes->ops->reset = SNESReset_NGMRES; 530 531 snes->usespc = PETSC_TRUE; 532 snes->usesksp = PETSC_FALSE; 533 534 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 535 snes->data = (void*) ngmres; 536 ngmres->msize = 30; 537 538 snes->max_funcs = 30000; 539 snes->max_its = 10000; 540 541 ngmres->additive = PETSC_FALSE; 542 ngmres->anderson = PETSC_FALSE; 543 544 ngmres->restart_it = 2; 545 ngmres->gammaA = 2.0; 546 ngmres->gammaC = 2.0; 547 ngmres->deltaB = 0.9; 548 ngmres->epsilonB = 0.1; 549 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 550 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 551 PetscFunctionReturn(0); 552 } 553 EXTERN_C_END 554