1 /* Defines the basic SNES object */ 2 #include <../src/snes/impls/ngmres/snesngmres.h> 3 #include <petscblaslapack.h> 4 5 6 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "SNESReset_NGMRES" 10 PetscErrorCode SNESReset_NGMRES(SNES snes) 11 { 12 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = VecDestroyVecs(ngmres->msize, &ngmres->rdot);CHKERRQ(ierr); 17 ierr = VecDestroyVecs(ngmres->msize, &ngmres->xdot);CHKERRQ(ierr); 18 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "SNESDestroy_NGMRES" 24 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 25 { 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 30 if (snes->data) { 31 SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data; 32 ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->r_norms, ngmres->q);CHKERRQ(ierr); 33 ierr = PetscFree(ngmres->s);CHKERRQ(ierr); 34 #if PETSC_USE_COMPLEX 35 ierr = PetscFree(ngmres->rwork); 36 #endif 37 ierr = PetscFree(ngmres->work); 38 } 39 ierr = PetscFree(snes->data); 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "SNESSetUp_NGMRES" 45 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 46 { 47 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 48 PetscInt msize,hsize; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 msize = ngmres->msize; /* restart size */ 53 hsize = msize * msize; 54 55 56 /* explicit least squares minimization solve */ 57 ierr = PetscMalloc5(hsize,PetscScalar,&ngmres->h, 58 msize,PetscScalar,&ngmres->beta, 59 msize,PetscScalar,&ngmres->xi, 60 msize,PetscReal,&ngmres->r_norms, 61 hsize,PetscScalar,&ngmres->q);CHKERRQ(ierr); 62 ngmres->nrhs = 1; 63 ngmres->lda = msize; 64 ngmres->ldb = msize; 65 ierr = PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);CHKERRQ(ierr); 66 67 ierr = PetscMemzero(ngmres->h,hsize*sizeof(PetscScalar));CHKERRQ(ierr); 68 ierr = PetscMemzero(ngmres->q,hsize*sizeof(PetscScalar));CHKERRQ(ierr); 69 ierr = PetscMemzero(ngmres->xi,msize*sizeof(PetscScalar));CHKERRQ(ierr); 70 ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr); 71 72 ngmres->lwork = 12*msize; 73 #if PETSC_USE_COMPLEX 74 ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork); 75 #endif 76 ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work); 77 78 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->xdot);CHKERRQ(ierr); 79 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->rdot);CHKERRQ(ierr); 80 ierr = SNESDefaultGetWork(snes, 3);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 86 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 87 { 88 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 89 PetscErrorCode ierr; 90 PetscBool debug; 91 92 PetscFunctionBegin; 93 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 94 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 95 ierr = PetscOptionsInt("-snes_ngmres_restart", "Maximum iterations before restart.", "SNES", ngmres->k_rmax, &ngmres->k_rmax, PETSC_NULL);CHKERRQ(ierr); 96 ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 97 if (debug) { 98 ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 99 } 100 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 101 ierr = PetscOptionsReal("-snes_ngmres_gammaC", "Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 102 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 103 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 104 ierr = PetscOptionsTail();CHKERRQ(ierr); 105 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 106 PetscFunctionReturn(0); 107 } 108 109 #undef __FUNCT__ 110 #define __FUNCT__ "SNESView_NGMRES" 111 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 112 { 113 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 114 PetscBool iascii; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 119 if (iascii) { 120 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPrintf(viewer, " Maximum iterations before restart %d\n", ngmres->k_rmax);CHKERRQ(ierr); 122 } 123 PetscFunctionReturn(0); 124 } 125 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "SNESSolve_NGMRES" 129 130 PetscErrorCode SNESSolve_NGMRES(SNES snes) 131 { 132 SNES pc; 133 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 134 135 /* present solution, residual, and preconditioned residual */ 136 Vec x, r, b, d; 137 Vec x_A, r_A; 138 139 /* previous iterations to construct the subspace */ 140 Vec *rdot = ngmres->rdot; 141 Vec *xdot = ngmres->xdot; 142 143 /* coefficients and RHS to the minimization problem */ 144 PetscScalar *beta = ngmres->beta; 145 PetscScalar *xi = ngmres->xi; 146 PetscReal r_norm, r_A_norm; 147 PetscReal nu; 148 PetscScalar alph_total = 0.; 149 PetscScalar qentry; 150 PetscInt i, j, k, k_restart, l, ivec; 151 152 /* solution selection data */ 153 PetscBool selectA, selectRestart; 154 PetscReal d_norm, d_min_norm, d_cur_norm; 155 PetscReal r_min_norm; 156 157 SNESConvergedReason reason; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 /* variable initialization */ 162 snes->reason = SNES_CONVERGED_ITERATING; 163 x = snes->vec_sol; 164 r = snes->vec_func; 165 b = snes->vec_rhs; 166 x_A = snes->vec_sol_update; 167 r_A = snes->work[0]; 168 d = snes->work[1]; 169 r = snes->work[2]; 170 171 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 172 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 173 snes->iter = 0; 174 snes->norm = 0.; 175 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 176 177 /* initialization */ 178 179 /* r = F(x) */ 180 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 181 if (snes->domainerror) { 182 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 183 PetscFunctionReturn(0); 184 } 185 186 /* nu = (r, r) */ 187 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 188 r_min_norm = r_norm; 189 nu = r_norm*r_norm; 190 if (PetscIsInfOrNanReal(r_norm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 191 192 /* q_{00} = nu */ 193 Q(0,0) = nu; 194 ngmres->r_norms[0] = r_norm; 195 /* rdot[0] = r */ 196 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 197 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 198 199 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 200 snes->norm = r_norm; 201 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 202 SNESLogConvHistory(snes, r_norm, 0); 203 ierr = SNESMonitor(snes, 0, r_norm);CHKERRQ(ierr); 204 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 205 if (snes->reason) PetscFunctionReturn(0); 206 207 k_restart = 1; 208 l = 1; 209 for (k=1; k<snes->max_its; k++) { 210 211 /* select which vector of the stored subspace will be updated */ 212 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 213 214 /* Computation of x^M */ 215 ierr = SNESSolve(pc, b, x);CHKERRQ(ierr); 216 ierr = SNESGetConvergedReason(pc,&reason);CHKERRQ(ierr); 217 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 218 snes->reason = SNES_DIVERGED_INNER; 219 PetscFunctionReturn(0); 220 } 221 222 /* r = F(x) */ 223 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 224 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 225 /* nu = (r, r) */ 226 ngmres->r_norms[ivec] = r_norm; 227 nu = r_norm*r_norm; 228 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^M */ 229 230 /* construct the right hand side and xi factors */ 231 for (i = 0; i < l; i++) { 232 ierr = VecDot(rdot[i], r, &xi[i]);CHKERRQ(ierr); 233 beta[i] = nu - xi[i]; 234 } 235 236 /* construct h */ 237 for (j = 0; j < l; j++) { 238 for (i = 0; i < l; i++) { 239 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 240 } 241 } 242 #ifdef PETSC_MISSING_LAPACK_GELSS 243 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 244 #else 245 ngmres->m = PetscBLASIntCast(l); 246 ngmres->n = PetscBLASIntCast(l); 247 ngmres->info = PetscBLASIntCast(0); 248 ngmres->rcond = -1.; 249 #ifdef PETSC_USE_COMPLEX 250 LAPACKgelss_(&ngmres->m, 251 &ngmres->n, 252 &ngmres->nrhs, 253 ngmres->h, 254 &ngmres->lda, 255 ngmres->beta, 256 &ngmres->ldb, 257 ngmres->s, 258 &ngmres->rcond, 259 &ngmres->rank, 260 ngmres->work, 261 &ngmres->lwork, 262 ngmres->rwork, 263 &ngmres->info); 264 #else 265 LAPACKgelss_(&ngmres->m, 266 &ngmres->n, 267 &ngmres->nrhs, 268 ngmres->h, 269 &ngmres->lda, 270 ngmres->beta, 271 &ngmres->ldb, 272 ngmres->s, 273 &ngmres->rcond, 274 &ngmres->rank, 275 ngmres->work, 276 &ngmres->lwork, 277 &ngmres->info); 278 #endif 279 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 280 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 281 #endif 282 283 alph_total = 0.; 284 for (i = 0; i < l; i++) { 285 alph_total += beta[i]; 286 } 287 ierr = VecCopy(x, x_A);CHKERRQ(ierr); 288 ierr = VecScale(x_A, 1. - alph_total);CHKERRQ(ierr); 289 290 for(i=0;i<l;i++){ 291 ierr= VecAXPY(x_A, beta[i], xdot[i]);CHKERRQ(ierr); 292 } 293 ierr = SNESComputeFunction(snes, x_A, r_A);CHKERRQ(ierr); 294 ierr = VecNorm(r_A, NORM_2, &r_A_norm);CHKERRQ(ierr); 295 296 selectA = PETSC_TRUE; 297 /* Conditions for choosing the accelerated answer */ 298 299 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 300 if (r_A_norm >= ngmres->gammaA*r_min_norm) { 301 selectA = PETSC_FALSE; 302 } 303 304 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 305 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 306 ierr=VecAXPY(d,-1,x);CHKERRQ(ierr); 307 ierr=VecNorm(d,NORM_2,&d_norm);CHKERRQ(ierr); 308 d_min_norm = -1.0; 309 for(i=0;i<l;i++) { 310 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 311 ierr=VecAXPY(d,-1,xdot[i]);CHKERRQ(ierr); 312 ierr=VecNorm(d,NORM_2,&d_cur_norm);CHKERRQ(ierr); 313 if((d_cur_norm < d_min_norm) || (d_min_norm < 0.0)) d_min_norm = d_cur_norm; 314 } 315 if (ngmres->epsilonB*d_norm<d_min_norm || sqrt(r_norm)<ngmres->deltaB*sqrt(r_min_norm)) { 316 } else { 317 selectA=PETSC_FALSE; 318 } 319 320 321 if (selectA) { 322 if (ngmres->monitor) { 323 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 324 } 325 /* copy it over */ 326 r_norm = r_A_norm; 327 nu = r_norm*r_norm; 328 ierr = VecCopy(r_A, r);CHKERRQ(ierr); 329 ierr = VecCopy(x_A, x);CHKERRQ(ierr); 330 } else { 331 if (ngmres->monitor) { 332 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 333 } 334 } 335 336 selectRestart = PETSC_FALSE; 337 338 /* maximum iteration criterion */ 339 if (k_restart > ngmres->k_rmax) { 340 selectRestart = PETSC_TRUE; 341 } 342 343 /* difference stagnation restart */ 344 if ((ngmres->epsilonB*d_norm > d_min_norm) && (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 345 if (ngmres->monitor) { 346 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm);CHKERRQ(ierr); 347 } 348 selectRestart = PETSC_TRUE; 349 } 350 351 /* residual stagnation restart */ 352 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 353 if (ngmres->monitor) { 354 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm));CHKERRQ(ierr); 355 } 356 selectRestart = PETSC_TRUE; 357 } 358 359 if (selectRestart) { 360 if (ngmres->monitor){ 361 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 362 } 363 k_restart = 1; 364 l = 1; 365 /* q_{00} = nu */ 366 ngmres->r_norms[0] = r_norm; 367 nu = r_norm*r_norm; 368 Q(0,0) = nu; 369 /* rdot[0] = r */ 370 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 371 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 372 } else { 373 /* select the current size of the subspace */ 374 if (l < ngmres->msize) l++; 375 k_restart++; 376 /* place the current entry in the list of previous entries */ 377 ierr = VecCopy(r, rdot[ivec]);CHKERRQ(ierr); 378 ierr = VecCopy(x, xdot[ivec]);CHKERRQ(ierr); 379 ngmres->r_norms[ivec] = r_norm; 380 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^A */ 381 for (i = 0; i < l; i++) { 382 ierr = VecDot(r, rdot[i], &qentry);CHKERRQ(ierr); 383 Q(i, ivec) = qentry; 384 Q(ivec, i) = qentry; 385 } 386 } 387 388 SNESLogConvHistory(snes, r_norm, k); 389 ierr = SNESMonitor(snes, k, r_norm);CHKERRQ(ierr); 390 391 snes->iter =k; 392 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 393 if (snes->reason) PetscFunctionReturn(0); 394 } 395 snes->reason = SNES_DIVERGED_MAX_IT; 396 PetscFunctionReturn(0); 397 } 398 399 /*MC 400 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 401 402 Level: beginner 403 404 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 405 SIAM Journal on Scientific Computing, 21(5), 2000. 406 407 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 408 J. Assoc. Comput. Mach., 12:547–560, 1965." 409 410 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 411 M*/ 412 413 EXTERN_C_BEGIN 414 #undef __FUNCT__ 415 #define __FUNCT__ "SNESCreate_NGMRES" 416 PetscErrorCode SNESCreate_NGMRES(SNES snes) 417 { 418 SNES_NGMRES *ngmres; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 snes->ops->destroy = SNESDestroy_NGMRES; 423 snes->ops->setup = SNESSetUp_NGMRES; 424 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 425 snes->ops->view = SNESView_NGMRES; 426 snes->ops->solve = SNESSolve_NGMRES; 427 snes->ops->reset = SNESReset_NGMRES; 428 429 snes->usespc = PETSC_TRUE; 430 snes->usesksp = PETSC_FALSE; 431 432 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 433 snes->data = (void*) ngmres; 434 ngmres->msize = 10; 435 436 ngmres->gammaA = 2.; 437 ngmres->gammaC = 2.; 438 ngmres->deltaB = 0.9; 439 ngmres->epsilonB = 0.1; 440 ngmres->k_rmax = 200; 441 442 if (!snes->pc) { 443 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 444 ierr = SNESSetType(snes->pc,SNESRICHARDSON);CHKERRQ(ierr); 445 } 446 PetscFunctionReturn(0); 447 } 448 EXTERN_C_END 449