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