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_CRITICAL: 144 ierr = SNESLineSearchSet(snes,SNESLineSearchCriticalSecant,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 method. 503 504 Level: beginner 505 506 Options Database: 507 + -snes_ngmres_additive - Use a variant that line-searches between the candidate solution and the combined solution. 508 . -snes_ngmres_anderson - Use Anderson mixing NGMRES variant which combines candidate solutions instead of actual solutions. 509 . -snes_ngmres_m - Number of stored previous solutions and residuals. 510 . -snes_ngmres_restart_it - Number of iterations the restart conditions hold before restart. 511 . -snes_ngmres_gammaA - Residual tolerance for solution selection between the candidate and combination. 512 . -snes_ngmres_gammaC - Residual tolerance for restart. 513 . -snes_ngmres_epsilonB - Difference tolerance between subsequent solutions triggering restart. 514 . -snes_ngmres_deltaB - Difference tolerance between residuals triggering restart. 515 . -snes_ngmres_monitor - Prints relevant information about the ngmres iteration. 516 - -snes_ls <basic,basicnonorms,quadratic,critical> - Line search type. 517 518 Notes: 519 520 The N-GMRES method combines m previous solutions into a minimum-residual solution by solving a small linearized 521 optimization problem at each iteration. 522 523 References: 524 525 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 526 SIAM Journal on Scientific Computing, 21(5), 2000. 527 528 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 529 J. Assoc. Comput. Mach., 12:547–560, 1965." 530 531 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 532 M*/ 533 534 EXTERN_C_BEGIN 535 #undef __FUNCT__ 536 #define __FUNCT__ "SNESCreate_NGMRES" 537 PetscErrorCode SNESCreate_NGMRES(SNES snes) 538 { 539 SNES_NGMRES *ngmres; 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 snes->ops->destroy = SNESDestroy_NGMRES; 544 snes->ops->setup = SNESSetUp_NGMRES; 545 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 546 snes->ops->view = SNESView_NGMRES; 547 snes->ops->solve = SNESSolve_NGMRES; 548 snes->ops->reset = SNESReset_NGMRES; 549 550 snes->usespc = PETSC_TRUE; 551 snes->usesksp = PETSC_FALSE; 552 553 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 554 snes->data = (void*) ngmres; 555 ngmres->msize = 30; 556 557 snes->max_funcs = 30000; 558 snes->max_its = 10000; 559 560 ngmres->additive = PETSC_FALSE; 561 ngmres->anderson = PETSC_FALSE; 562 563 ngmres->restart_it = 2; 564 ngmres->gammaA = 2.0; 565 ngmres->gammaC = 2.0; 566 ngmres->deltaB = 0.9; 567 ngmres->epsilonB = 0.1; 568 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr); 569 ierr = SNESLineSearchSetType(snes, SNES_LS_BASIC);CHKERRQ(ierr); 570 PetscFunctionReturn(0); 571 } 572 EXTERN_C_END 573