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 ierr = PetscMemzero(ngmres->h, hsize*sizeof(PetscScalar));CHKERRQ(ierr); 67 ierr = PetscMemzero(ngmres->q, hsize*sizeof(PetscScalar));CHKERRQ(ierr); 68 ierr = PetscMemzero(ngmres->xi, msize*sizeof(PetscScalar));CHKERRQ(ierr); 69 ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr); 70 ngmres->lwork = 12*msize; 71 #if PETSC_USE_COMPLEX 72 ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork); 73 #endif 74 ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work); 75 76 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->xdot);CHKERRQ(ierr); 77 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->rdot);CHKERRQ(ierr); 78 ierr = SNESDefaultGetWork(snes, 3);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 82 #undef __FUNCT__ 83 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 84 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 85 { 86 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 87 PetscErrorCode ierr; 88 PetscBool debug; 89 90 PetscFunctionBegin; 91 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 92 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 93 ierr = PetscOptionsInt("-snes_ngmres_restart", "Maximum iterations before restart.", "SNES", ngmres->k_rmax, &ngmres->k_rmax, PETSC_NULL);CHKERRQ(ierr); 94 ierr = PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr); 95 if (debug) { 96 ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 97 } 98 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 99 ierr = PetscOptionsReal("-snes_ngmres_gammaC", " Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 100 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 101 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 102 ierr = PetscOptionsTail();CHKERRQ(ierr); 103 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "SNESView_NGMRES" 109 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 110 { 111 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 112 PetscBool iascii; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 117 if (iascii) { 118 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer, " Maximum iterations before restart %d\n", ngmres->k_rmax);CHKERRQ(ierr); 120 } 121 PetscFunctionReturn(0); 122 } 123 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "SNESSolve_NGMRES" 127 128 PetscErrorCode SNESSolve_NGMRES(SNES snes) 129 { 130 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 131 /* present solution, residual, and preconditioned residual */ 132 Vec x, r, b, d; 133 Vec x_A, r_A; 134 135 /* previous iterations to construct the subspace */ 136 Vec *rdot = ngmres->rdot; 137 Vec *xdot = ngmres->xdot; 138 139 /* coefficients and RHS to the minimization problem */ 140 PetscScalar *beta = ngmres->beta; 141 PetscScalar *xi = ngmres->xi; 142 PetscReal r_norm, r_A_norm; 143 PetscReal nu; 144 PetscScalar alph_total = 0.; 145 PetscScalar qentry; 146 PetscInt i, j, k, k_restart, l, ivec; 147 148 /* solution selection data */ 149 PetscBool selectA, selectRestart; 150 PetscReal d_norm, d_min_norm, d_cur_norm; 151 PetscReal r_min_norm; 152 153 SNESConvergedReason reason; 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 /* variable initialization */ 158 snes->reason = SNES_CONVERGED_ITERATING; 159 x = snes->vec_sol; 160 r = snes->vec_func; 161 b = snes->vec_rhs; 162 x_A = snes->vec_sol_update; 163 r_A = snes->work[0]; 164 d = snes->work[1]; 165 166 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 167 snes->iter = 0; 168 snes->norm = 0.; 169 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 170 171 /* initialization */ 172 173 /* r = F(x) */ 174 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 175 if (snes->domainerror) { 176 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 177 PetscFunctionReturn(0); 178 } 179 180 /* nu = (r, r) */ 181 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 182 r_min_norm = r_norm; 183 nu = r_norm*r_norm; 184 if (PetscIsInfOrNanReal(r_norm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation"); 185 186 /* q_{00} = nu */ 187 Q(0,0) = nu; 188 ngmres->r_norms[0] = r_norm; 189 /* rdot[0] = r */ 190 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 191 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 192 193 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 194 snes->norm = r_norm; 195 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 196 SNESLogConvHistory(snes, r_norm, 0); 197 ierr = SNESMonitor(snes, 0, r_norm);CHKERRQ(ierr); 198 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 199 if (snes->reason) PetscFunctionReturn(0); 200 201 k_restart = 1; 202 l = 1; 203 for (k=1; k<snes->max_its; k++) { 204 /* select which vector of the stored subspace will be updated */ 205 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 206 207 /* Computation of x^M */ 208 if (!snes->pc) { 209 /* no preconditioner -- just take gradient descent */ 210 ierr = VecAXPY(x, -1.0, r);CHKERRQ(ierr); 211 } else { 212 ierr = SNESSolve(snes->pc, b, x);CHKERRQ(ierr); 213 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 214 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 215 snes->reason = SNES_DIVERGED_INNER; 216 PetscFunctionReturn(0); 217 } 218 } 219 220 /* r = F(x) */ 221 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 222 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 223 /* nu = (r, r) */ 224 ngmres->r_norms[ivec] = r_norm; 225 nu = r_norm*r_norm; 226 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^M */ 227 228 /* construct the right hand side and xi factors */ 229 for (i = 0; i < l; i++) { 230 ierr = VecDot(rdot[i], r, &xi[i]);CHKERRQ(ierr); 231 beta[i] = nu - xi[i]; 232 } 233 234 /* construct h */ 235 for (j = 0; j < l; j++) { 236 for (i = 0; i < l; i++) { 237 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 238 } 239 } 240 #ifdef PETSC_MISSING_LAPACK_GELSS 241 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 242 #else 243 ngmres->m = PetscBLASIntCast(l); 244 ngmres->n = PetscBLASIntCast(l); 245 ngmres->info = PetscBLASIntCast(0); 246 ngmres->rcond = -1.; 247 #ifdef PETSC_USE_COMPLEX 248 LAPACKgelss_(&ngmres->m, 249 &ngmres->n, 250 &ngmres->nrhs, 251 ngmres->h, 252 &ngmres->lda, 253 ngmres->beta, 254 &ngmres->ldb, 255 ngmres->s, 256 &ngmres->rcond, 257 &ngmres->rank, 258 ngmres->work, 259 &ngmres->lwork, 260 ngmres->rwork, 261 &ngmres->info); 262 #else 263 LAPACKgelss_(&ngmres->m, 264 &ngmres->n, 265 &ngmres->nrhs, 266 ngmres->h, 267 &ngmres->lda, 268 ngmres->beta, 269 &ngmres->ldb, 270 ngmres->s, 271 &ngmres->rcond, 272 &ngmres->rank, 273 ngmres->work, 274 &ngmres->lwork, 275 &ngmres->info); 276 #endif 277 if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS"); 278 if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge"); 279 #endif 280 281 alph_total = 0.; 282 for (i = 0; i < l; i++) { 283 alph_total += beta[i]; 284 } 285 ierr = VecCopy(x, x_A);CHKERRQ(ierr); 286 ierr = VecScale(x_A, 1. - alph_total);CHKERRQ(ierr); 287 288 for(i=0;i<l;i++){ 289 ierr= VecAXPY(x_A, beta[i], xdot[i]);CHKERRQ(ierr); 290 } 291 ierr = SNESComputeFunction(snes, x_A, r_A);CHKERRQ(ierr); 292 ierr = VecNorm(r_A, NORM_2, &r_A_norm);CHKERRQ(ierr); 293 294 selectA = PETSC_TRUE; 295 /* Conditions for choosing the accelerated answer */ 296 297 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 298 if (r_A_norm >= ngmres->gammaA*r_min_norm) { 299 selectA = PETSC_FALSE; 300 } 301 302 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 303 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 304 ierr=VecAXPY(d,-1,x);CHKERRQ(ierr); 305 ierr=VecNorm(d,NORM_2,&d_norm);CHKERRQ(ierr); 306 d_min_norm = -1.0; 307 for(i=0;i<l;i++) { 308 ierr=VecCopy(x_A,d);CHKERRQ(ierr); 309 ierr=VecAXPY(d,-1,xdot[i]);CHKERRQ(ierr); 310 ierr=VecNorm(d,NORM_2,&d_cur_norm);CHKERRQ(ierr); 311 if((d_cur_norm < d_min_norm) || (d_min_norm < 0.0)) d_min_norm = d_cur_norm; 312 } 313 if (ngmres->epsilonB*d_norm<d_min_norm || sqrt(r_norm)<ngmres->deltaB*sqrt(r_min_norm)) { 314 } else { 315 selectA=PETSC_FALSE; 316 } 317 318 319 if (selectA) { 320 if (ngmres->monitor) { 321 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 322 } 323 /* copy it over */ 324 r_norm = r_A_norm; 325 nu = r_norm*r_norm; 326 ierr = VecCopy(r_A, r);CHKERRQ(ierr); 327 ierr = VecCopy(x_A, x);CHKERRQ(ierr); 328 } else { 329 if (ngmres->monitor) { 330 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm);CHKERRQ(ierr); 331 } 332 } 333 334 selectRestart = PETSC_FALSE; 335 336 /* maximum iteration criterion */ 337 if (k_restart > ngmres->k_rmax) { 338 selectRestart = PETSC_TRUE; 339 } 340 341 /* difference stagnation restart */ 342 if((ngmres->epsilonB*d_norm > d_min_norm) && (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 343 if (ngmres->monitor) { 344 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm);CHKERRQ(ierr); 345 } 346 selectRestart = PETSC_TRUE; 347 } 348 /* residual stagnation restart */ 349 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 350 if (ngmres->monitor) { 351 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm));CHKERRQ(ierr); 352 } 353 selectRestart = PETSC_TRUE; 354 } 355 356 if (selectRestart) { 357 if (ngmres->monitor){ 358 ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr); 359 } 360 k_restart = 1; 361 l = 1; 362 /* q_{00} = nu */ 363 ngmres->r_norms[0] = r_norm; 364 nu = r_norm*r_norm; 365 Q(0,0) = nu; 366 /* rdot[0] = r */ 367 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 368 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 369 } else { 370 /* select the current size of the subspace */ 371 if (l < ngmres->msize) l++; 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 ierr = VecDot(r, rdot[i], &qentry);CHKERRQ(ierr); 380 Q(i, ivec) = qentry; 381 Q(ivec, i) = qentry; 382 } 383 } 384 385 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 386 snes->iter = k; 387 snes->norm = r_norm; 388 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 389 SNESLogConvHistory(snes, snes->norm, snes->iter); 390 ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr); 391 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 PetscFunctionReturn(0); 443 } 444 EXTERN_C_END 445