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 ngmres->nrhs = 1; 245 #ifdef PETSC_USE_COMPLEX 246 LAPACKgelss_(&ngmres->m, 247 &ngmres->n, 248 &ngmres->nrhs, 249 ngmres->h, 250 &ngmres->lda, 251 ngmres->beta, 252 &ngmres->ldb, 253 ngmres->s, 254 &ngmres->rcond, 255 &ngmres->rank, 256 ngmres->work, 257 &ngmres->lwork, 258 ngmres->rwork, 259 &ngmres->info); 260 #else 261 LAPACKgelss_(&ngmres->m, 262 &ngmres->n, 263 &ngmres->nrhs, 264 ngmres->h, 265 &ngmres->lda, 266 ngmres->beta, 267 &ngmres->ldb, 268 ngmres->s, 269 &ngmres->rcond, 270 &ngmres->rank, 271 ngmres->work, 272 &ngmres->lwork, 273 &ngmres->info); 274 #endif 275 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 276 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 277 #endif 278 279 alph_total = 0.; 280 for (i = 0; i < l; i++) { 281 alph_total += beta[i]; 282 } 283 ierr = VecCopy(x, x_A);CHKERRQ(ierr); 284 ierr = VecScale(x_A, 1. - alph_total);CHKERRQ(ierr); 285 286 for(i=0;i<l;i++){ 287 ierr= VecAXPY(x_A, beta[i], xdot[i]);CHKERRQ(ierr); 288 } 289 ierr = SNESComputeFunction(snes, x_A, r_A);CHKERRQ(ierr); 290 ierr = VecNorm(r_A, NORM_2, &r_A_norm);CHKERRQ(ierr); 291 292 selectA = PETSC_TRUE; 293 /* Conditions for choosing the accelerated answer */ 294 295 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 296 if (r_A_norm >= ngmres->gammaA*r_min_norm) { 297 selectA = PETSC_FALSE; 298 } 299 300 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 301 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 302 ierr=VecAXPY(d,-1,x);CHKERRQ(ierr); 303 ierr=VecNorm(d,NORM_2,&d_norm);CHKERRQ(ierr); 304 d_min_norm = -1.0; 305 for(i=0;i<l;i++) { 306 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 307 ierr=VecAXPY(d,-1,xdot[i]);CHKERRQ(ierr); 308 ierr=VecNorm(d,NORM_2,&d_cur_norm);CHKERRQ(ierr); 309 if((d_cur_norm < d_min_norm) || (d_min_norm < 0.0)) d_min_norm = d_cur_norm; 310 } 311 if (ngmres->epsilonB*d_norm<d_min_norm || sqrt(r_norm)<ngmres->deltaB*sqrt(r_min_norm)) { 312 } else { 313 selectA=PETSC_FALSE; 314 } 315 316 317 if (selectA) { 318 if (ngmres->monitor) { 319 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 320 } 321 /* copy it over */ 322 r_norm = r_A_norm; 323 nu = r_norm*r_norm; 324 ierr = VecCopy(r_A, r);CHKERRQ(ierr); 325 ierr = VecCopy(x_A, x);CHKERRQ(ierr); 326 } else { 327 if (ngmres->monitor) { 328 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 329 } 330 } 331 332 selectRestart = PETSC_FALSE; 333 334 /* maximum iteration criterion */ 335 if (k_restart > ngmres->k_rmax) { 336 selectRestart = PETSC_TRUE; 337 } 338 339 /* difference stagnation restart */ 340 if ((ngmres->epsilonB*d_norm > d_min_norm) && (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 341 if (ngmres->monitor) { 342 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm);CHKERRQ(ierr); 343 } 344 selectRestart = PETSC_TRUE; 345 } 346 347 /* residual stagnation restart */ 348 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 349 if (ngmres->monitor) { 350 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm));CHKERRQ(ierr); 351 } 352 selectRestart = PETSC_TRUE; 353 } 354 355 if (selectRestart) { 356 if (ngmres->monitor){ 357 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 358 } 359 k_restart = 1; 360 l = 1; 361 /* q_{00} = nu */ 362 ngmres->r_norms[0] = r_norm; 363 nu = r_norm*r_norm; 364 Q(0,0) = nu; 365 /* rdot[0] = r */ 366 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 367 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 368 } else { 369 /* select the current size of the subspace */ 370 if (l < ngmres->msize) { 371 l++; 372 } 373 k_restart++; 374 /* place the current entry in the list of previous entries */ 375 ierr = VecCopy(r, rdot[ivec]);CHKERRQ(ierr); 376 ierr = VecCopy(x, xdot[ivec]);CHKERRQ(ierr); 377 ngmres->r_norms[ivec] = r_norm; 378 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^A */ 379 for (i = 0; i < l; i++) { 380 VecDot(r, rdot[i], &qentry); 381 Q(i, ivec) = qentry; 382 Q(ivec, i) = qentry; 383 } 384 } 385 386 SNESLogConvHistory(snes, r_norm, k); 387 ierr = SNESMonitor(snes, k, r_norm);CHKERRQ(ierr); 388 389 snes->iter =k; 390 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 391 if (snes->reason) PetscFunctionReturn(0); 392 } 393 snes->reason = SNES_DIVERGED_MAX_IT; 394 PetscFunctionReturn(0); 395 } 396 397 /*MC 398 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 399 400 Level: beginner 401 402 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 403 SIAM Journal on Scientific Computing, 21(5), 2000. 404 405 This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations. 406 J. Assoc. Comput. Mach., 12:547–560, 1965." 407 408 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 409 M*/ 410 411 EXTERN_C_BEGIN 412 #undef __FUNCT__ 413 #define __FUNCT__ "SNESCreate_NGMRES" 414 PetscErrorCode SNESCreate_NGMRES(SNES snes) 415 { 416 SNES_NGMRES *ngmres; 417 PetscErrorCode ierr; 418 419 PetscFunctionBegin; 420 snes->ops->destroy = SNESDestroy_NGMRES; 421 snes->ops->setup = SNESSetUp_NGMRES; 422 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 423 snes->ops->view = SNESView_NGMRES; 424 snes->ops->solve = SNESSolve_NGMRES; 425 snes->ops->reset = SNESReset_NGMRES; 426 427 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 428 snes->data = (void*) ngmres; 429 ngmres->msize = 10; 430 431 ngmres->gammaA = 2.; 432 ngmres->gammaC = 2.; 433 ngmres->deltaB = 0.9; 434 ngmres->epsilonB = 0.1; 435 ngmres->k_rmax = 200; 436 437 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 EXTERN_C_END 441