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