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 136 137 /* present solution, residual, and preconditioned residual */ 138 Vec x, r, b, d; 139 Vec x_A, r_A; 140 141 /* previous iterations to construct the subspace */ 142 Vec *rdot = ngmres->rdot; 143 Vec *xdot = ngmres->xdot; 144 145 /* coefficients and RHS to the minimization problem */ 146 PetscScalar *beta = ngmres->beta; 147 PetscScalar *xi = ngmres->xi; 148 PetscReal r_norm, r_A_norm; 149 PetscReal nu; 150 PetscScalar alph_total = 0.; 151 PetscScalar qentry; 152 PetscInt i, j, k, k_restart, l, ivec; 153 154 /* solution selection data */ 155 PetscBool selectA, selectRestart; 156 PetscReal d_norm, d_min_norm, d_cur_norm; 157 PetscReal r_min_norm; 158 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 /* variable initialization */ 163 snes->reason = SNES_CONVERGED_ITERATING; 164 x = snes->vec_sol; 165 r = snes->vec_func; 166 b = snes->vec_rhs; 167 x_A = snes->vec_sol_update; 168 r_A = snes->work[0]; 169 d = snes->work[1]; 170 r = snes->work[2]; 171 172 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 173 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 174 snes->iter = 0; 175 snes->norm = 0.; 176 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 177 178 /* initialization */ 179 180 /* r = F(x) */ 181 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 182 if (snes->domainerror) { 183 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 184 PetscFunctionReturn(0); 185 } 186 187 /* nu = (r, r) */ 188 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 189 r_min_norm = r_norm; 190 nu = r_norm*r_norm; 191 if (PetscIsInfOrNanReal(r_norm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 192 193 /* q_{00} = nu */ 194 Q(0,0) = nu; 195 ngmres->r_norms[0] = r_norm; 196 /* rdot[0] = r */ 197 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 198 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 199 200 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 201 snes->norm = r_norm; 202 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 203 SNESLogConvHistory(snes, r_norm, 0); 204 ierr = SNESMonitor(snes, 0, r_norm);CHKERRQ(ierr); 205 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 206 if (snes->reason) PetscFunctionReturn(0); 207 208 k_restart = 1; 209 l = 1; 210 for (k=1; k<snes->max_its; k++) { 211 212 /* select which vector of the stored subspace will be updated */ 213 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 214 215 /* Computation of x^M */ 216 ierr = SNESSolve(pc, b, x);CHKERRQ(ierr); 217 /* r = F(x) */ 218 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 219 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 220 /* nu = (r, r) */ 221 ngmres->r_norms[ivec] = r_norm; 222 nu = r_norm*r_norm; 223 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^M */ 224 225 /* construct the right hand side and xi factors */ 226 for (i = 0; i < l; i++) { 227 VecDot(rdot[i], r, &xi[i]); 228 beta[i] = nu - xi[i]; 229 } 230 231 /* construct h */ 232 for (j = 0; j < l; j++) { 233 for (i = 0; i < l; i++) { 234 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 235 } 236 } 237 #ifdef PETSC_MISSING_LAPACK_GELSS 238 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 239 #else 240 ngmres->m = PetscBLASIntCast(l); 241 ngmres->n = PetscBLASIntCast(l); 242 ngmres->info = PetscBLASIntCast(0); 243 ngmres->rcond = -1.; 244 #ifdef PETSC_USE_COMPLEX 245 LAPACKgelss_(&ngmres->m, 246 &ngmres->n, 247 &ngmres->nrhs, 248 ngmres->h, 249 &ngmres->lda, 250 ngmres->beta, 251 &ngmres->ldb, 252 ngmres->s, 253 &ngmres->rcond, 254 &ngmres->rank, 255 ngmres->work, 256 &ngmres->lwork, 257 ngmres->rwork, 258 &ngmres->info); 259 #else 260 LAPACKgelss_(&ngmres->m, 261 &ngmres->n, 262 &ngmres->nrhs, 263 ngmres->h, 264 &ngmres->lda, 265 ngmres->beta, 266 &ngmres->ldb, 267 ngmres->s, 268 &ngmres->rcond, 269 &ngmres->rank, 270 ngmres->work, 271 &ngmres->lwork, 272 &ngmres->info); 273 #endif 274 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 275 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 276 #endif 277 278 alph_total = 0.; 279 for (i = 0; i < l; i++) { 280 alph_total += beta[i]; 281 } 282 ierr = VecCopy(x, x_A);CHKERRQ(ierr); 283 ierr = VecScale(x_A, 1. - alph_total);CHKERRQ(ierr); 284 285 for(i=0;i<l;i++){ 286 ierr= VecAXPY(x_A, beta[i], xdot[i]);CHKERRQ(ierr); 287 } 288 ierr = SNESComputeFunction(snes, x_A, r_A);CHKERRQ(ierr); 289 ierr = VecNorm(r_A, NORM_2, &r_A_norm);CHKERRQ(ierr); 290 291 selectA = PETSC_TRUE; 292 /* Conditions for choosing the accelerated answer */ 293 294 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 295 if (r_A_norm >= ngmres->gammaA*r_min_norm) { 296 selectA = PETSC_FALSE; 297 } 298 299 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 300 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 301 ierr=VecAXPY(d,-1,x);CHKERRQ(ierr); 302 ierr=VecNorm(d,NORM_2,&d_norm);CHKERRQ(ierr); 303 d_min_norm = -1.0; 304 for(i=0;i<l;i++) { 305 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 306 ierr=VecAXPY(d,-1,xdot[i]);CHKERRQ(ierr); 307 ierr=VecNorm(d,NORM_2,&d_cur_norm);CHKERRQ(ierr); 308 if((d_cur_norm < d_min_norm) || (d_min_norm < 0.0)) d_min_norm = d_cur_norm; 309 } 310 if (ngmres->epsilonB*d_norm<d_min_norm || sqrt(r_norm)<ngmres->deltaB*sqrt(r_min_norm)) { 311 } else { 312 selectA=PETSC_FALSE; 313 } 314 315 316 if (selectA) { 317 if (ngmres->monitor) { 318 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 319 } 320 /* copy it over */ 321 r_norm = r_A_norm; 322 nu = r_norm*r_norm; 323 ierr = VecCopy(r_A, r);CHKERRQ(ierr); 324 ierr = VecCopy(x_A, x);CHKERRQ(ierr); 325 } else { 326 if (ngmres->monitor) { 327 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 328 } 329 } 330 331 selectRestart = PETSC_FALSE; 332 333 /* maximum iteration criterion */ 334 if (k_restart > ngmres->k_rmax) { 335 selectRestart = PETSC_TRUE; 336 } 337 338 /* difference stagnation restart */ 339 if ((ngmres->epsilonB*d_norm > d_min_norm) && (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 340 if (ngmres->monitor) { 341 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm);CHKERRQ(ierr); 342 } 343 selectRestart = PETSC_TRUE; 344 } 345 346 /* residual stagnation restart */ 347 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 348 if (ngmres->monitor) { 349 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm));CHKERRQ(ierr); 350 } 351 selectRestart = PETSC_TRUE; 352 } 353 354 if (selectRestart) { 355 if (ngmres->monitor){ 356 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 357 } 358 k_restart = 1; 359 l = 1; 360 /* q_{00} = nu */ 361 ngmres->r_norms[0] = r_norm; 362 nu = r_norm*r_norm; 363 Q(0,0) = nu; 364 /* rdot[0] = r */ 365 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 366 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 367 } else { 368 /* select the current size of the subspace */ 369 if (l < ngmres->msize) { 370 l++; 371 } 372 k_restart++; 373 /* place the current entry in the list of previous entries */ 374 ierr = VecCopy(r, rdot[ivec]);CHKERRQ(ierr); 375 ierr = VecCopy(x, xdot[ivec]);CHKERRQ(ierr); 376 ngmres->r_norms[ivec] = r_norm; 377 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^A */ 378 for (i = 0; i < l; i++) { 379 VecDot(r, rdot[i], &qentry); 380 Q(i, ivec) = qentry; 381 Q(ivec, i) = qentry; 382 } 383 } 384 385 SNESLogConvHistory(snes, r_norm, k); 386 ierr = SNESMonitor(snes, k, r_norm);CHKERRQ(ierr); 387 388 snes->iter =k; 389 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 390 if (snes->reason) PetscFunctionReturn(0); 391 } 392 snes->reason = SNES_DIVERGED_MAX_IT; 393 PetscFunctionReturn(0); 394 } 395 396 /*MC 397 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 398 399 Level: beginner 400 401 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 402 SIAM Journal on Scientific Computing, 21(5), 2000. 403 404 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 405 J. Assoc. Comput. Mach., 12:547–560, 1965." 406 407 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 408 M*/ 409 410 EXTERN_C_BEGIN 411 #undef __FUNCT__ 412 #define __FUNCT__ "SNESCreate_NGMRES" 413 PetscErrorCode SNESCreate_NGMRES(SNES snes) 414 { 415 SNES_NGMRES *ngmres; 416 PetscErrorCode ierr; 417 418 PetscFunctionBegin; 419 snes->ops->destroy = SNESDestroy_NGMRES; 420 snes->ops->setup = SNESSetUp_NGMRES; 421 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 422 snes->ops->view = SNESView_NGMRES; 423 snes->ops->solve = SNESSolve_NGMRES; 424 snes->ops->reset = SNESReset_NGMRES; 425 426 snes->usesksp = PETSC_FALSE; 427 428 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 429 snes->data = (void*) ngmres; 430 ngmres->msize = 10; 431 432 ngmres->gammaA = 2.; 433 ngmres->gammaC = 2.; 434 ngmres->deltaB = 0.9; 435 ngmres->epsilonB = 0.1; 436 ngmres->k_rmax = 200; 437 438 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 439 PetscFunctionReturn(0); 440 } 441 EXTERN_C_END 442