1 /* Defines the basic SNES object */ 2 #include <../src/snes/impls/ngmres/snesngmres.h> 3 #include <petscblaslapack.h> 4 5 6 /*MC 7 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 8 9 Level: beginner 10 11 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 12 SIAM Journal on Scientific Computing, 21(5), 2000. 13 14 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 15 M*/ 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "SNESReset_NGMRES" 19 PetscErrorCode SNESReset_NGMRES(SNES snes) 20 { 21 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = VecDestroyVecs(ngmres->msize, &ngmres->rdot);CHKERRQ(ierr); 26 ierr = VecDestroyVecs(ngmres->msize, &ngmres->xdot);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "SNESDestroy_NGMRES" 32 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 33 { 34 PetscErrorCode ierr; 35 36 PetscFunctionBegin; 37 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 38 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 39 if (snes->data) { 40 SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data; 41 ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->r_norms, ngmres->q);CHKERRQ(ierr); 42 ierr = PetscFree(ngmres->s);CHKERRQ(ierr); 43 #if PETSC_USE_COMPLEX 44 ierr = PetscFree(ngmres->rwork); 45 #endif 46 ierr = PetscFree(ngmres->work); 47 } 48 ierr = PetscFree(snes->data); 49 PetscFunctionReturn(0); 50 } 51 52 #undef __FUNCT__ 53 #define __FUNCT__ "SNESSetUp_NGMRES" 54 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 55 { 56 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 57 PetscInt msize,hsize; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 msize = ngmres->msize; /* restart size */ 62 hsize = msize * msize; 63 64 65 /* explicit least squares minimization solve */ 66 ierr = PetscMalloc5(hsize,PetscScalar,&ngmres->h, 67 msize,PetscScalar,&ngmres->beta, 68 msize,PetscScalar,&ngmres->xi, 69 msize,PetscReal,&ngmres->r_norms, 70 hsize,PetscScalar,&ngmres->q);CHKERRQ(ierr); 71 ngmres->nrhs = 1; 72 ngmres->lda = msize; 73 ngmres->ldb = msize; 74 ierr = PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);CHKERRQ(ierr); 75 76 ierr = PetscMemzero(ngmres->h,hsize*sizeof(PetscScalar));CHKERRQ(ierr); 77 ierr = PetscMemzero(ngmres->q,hsize*sizeof(PetscScalar));CHKERRQ(ierr); 78 ierr = PetscMemzero(ngmres->xi,msize*sizeof(PetscScalar));CHKERRQ(ierr); 79 ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr); 80 81 ngmres->lwork = 12*msize; 82 #if PETSC_USE_COMPLEX 83 ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork); 84 #endif 85 ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work); 86 87 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->xdot);CHKERRQ(ierr); 88 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->rdot);CHKERRQ(ierr); 89 ierr = SNESDefaultGetWork(snes, 2);CHKERRQ(ierr); 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 95 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 96 { 97 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 102 ierr = PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 103 ierr = PetscOptionsInt("-snes_ngmres_restart", "Maximum iterations before restart.", "SNES", ngmres->k_rmax, &ngmres->k_rmax, PETSC_NULL);CHKERRQ(ierr); 104 ierr = PetscOptionsBool("-snes_ngmres_debug", "Debugging output for NGMRES", "SNES", ngmres->debug, &ngmres->debug, PETSC_NULL);CHKERRQ(ierr); 105 ierr = PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr); 106 ierr = PetscOptionsReal("-snes_ngmres_gammaC", "Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr); 107 ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr); 108 ierr = PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr); 109 ierr = PetscOptionsTail();CHKERRQ(ierr); 110 if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA; 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "SNESView_NGMRES" 116 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 117 { 118 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 119 PetscBool iascii; 120 PetscErrorCode ierr; 121 122 PetscFunctionBegin; 123 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 124 if (iascii) { 125 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(viewer, " Maximum iterations before restart %d\n", ngmres->k_rmax);CHKERRQ(ierr); 127 } 128 PetscFunctionReturn(0); 129 } 130 131 132 #undef __FUNCT__ 133 #define __FUNCT__ "SNESSolve_NGMRES" 134 135 PetscErrorCode SNESSolve_NGMRES(SNES snes) 136 { 137 SNES pc; 138 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 139 140 141 142 /* present solution, residual, and preconditioned residual */ 143 Vec x, r, b, d; 144 Vec x_A, r_A; 145 146 /* previous iterations to construct the subspace */ 147 Vec *rdot = ngmres->rdot; 148 Vec *xdot = ngmres->xdot; 149 150 /* coefficients and RHS to the minimization problem */ 151 PetscScalar *beta = ngmres->beta; 152 PetscScalar *xi = ngmres->xi; 153 PetscReal r_norm, r_A_norm; 154 PetscReal nu; 155 PetscScalar alph_total = 0.; 156 PetscScalar qentry; 157 PetscInt i, j, k, k_restart, l, ivec; 158 159 /* solution selection data */ 160 PetscBool selectA, selectRestart; 161 PetscReal d_norm, d_min_norm, d_cur_norm; 162 PetscReal r_min_norm; 163 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 /* variable initialization */ 168 snes->reason = SNES_CONVERGED_ITERATING; 169 x = snes->vec_sol; 170 r = snes->vec_func; 171 b = snes->vec_rhs; 172 x_A = snes->vec_sol_update; 173 r_A = snes->work[0]; 174 d = snes->work[1]; 175 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 176 177 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 178 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 179 snes->iter = 0; 180 snes->norm = 0.; 181 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 182 183 /* initialization */ 184 185 /* r = F(x) */ 186 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 187 if (snes->domainerror) { 188 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 189 PetscFunctionReturn(0); 190 } 191 192 /* nu = (r, r) */ 193 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 194 r_min_norm = r_norm; 195 nu = r_norm*r_norm; 196 if (PetscIsInfOrNanReal(r_norm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 197 198 /* q_{00} = nu */ 199 Q(0,0) = nu; 200 ngmres->r_norms[0] = r_norm; 201 /* rdot[0] = r */ 202 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 203 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 204 205 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 206 snes->norm = r_norm; 207 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 208 SNESLogConvHistory(snes, r_norm, 0); 209 ierr = SNESMonitor(snes, 0, r_norm);CHKERRQ(ierr); 210 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 211 if (snes->reason) PetscFunctionReturn(0); 212 213 k_restart = 1; 214 l = 1; 215 for (k=1; k<snes->max_its; k++) { 216 217 /* select which vector of the stored subspace will be updated */ 218 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 219 220 221 /* Computation of x^M */ 222 ierr = SNESSolve(pc, b, x);CHKERRQ(ierr); 223 /* r = F(x) */ 224 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 225 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 226 /* nu = (r, r) */ 227 ngmres->r_norms[ivec] = r_norm; 228 nu = r_norm*r_norm; 229 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^M */ 230 231 /* construct the right hand side and xi factors */ 232 for (i = 0; i < l; i++) { 233 VecDot(rdot[i], r, &xi[i]); 234 beta[i] = nu - xi[i]; 235 } 236 237 /* construct h */ 238 for (j = 0; j < l; j++) { 239 for (i = 0; i < l; i++) { 240 H(i, j) = Q(i, j) - xi[i] - xi[j] + nu; 241 } 242 } 243 #ifdef PETSC_MISSING_LAPACK_GELSS 244 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine."); 245 #else 246 ngmres->m = PetscBLASIntCast(l); 247 ngmres->n = PetscBLASIntCast(l); 248 ngmres->info = PetscBLASIntCast(0); 249 ngmres->rcond = -1.; 250 ngmres->nrhs = 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->debug) 325 PetscPrintf(PETSC_COMM_WORLD, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm); 326 /* copy it over */ 327 r_norm = r_A_norm; 328 nu = r_norm*r_norm; 329 ierr = VecCopy(r_A, r);CHKERRQ(ierr); 330 ierr = VecCopy(x_A, x);CHKERRQ(ierr); 331 } else { 332 if(ngmres->debug) 333 PetscPrintf(PETSC_COMM_WORLD, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm); 334 } 335 336 selectRestart = PETSC_FALSE; 337 338 /* maximum iteration criterion */ 339 if (k_restart > ngmres->k_rmax) { 340 selectRestart = PETSC_TRUE; 341 } 342 343 /* difference stagnation restart */ 344 if ((ngmres->epsilonB*d_norm > d_min_norm) && 345 (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 346 if (ngmres->debug) 347 PetscPrintf(PETSC_COMM_WORLD, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm); 348 selectRestart = PETSC_TRUE; 349 } 350 351 /* residual stagnation restart */ 352 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 353 if (ngmres->debug) 354 PetscPrintf(PETSC_COMM_WORLD, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm)); 355 selectRestart = PETSC_TRUE; 356 } 357 358 if (selectRestart) { 359 if (ngmres->debug) 360 PetscPrintf(PETSC_COMM_WORLD, "Restarted at iteration %d\n", k_restart); 361 k_restart = 1; 362 l = 1; 363 /* q_{00} = nu */ 364 ngmres->r_norms[0] = r_norm; 365 nu = r_norm*r_norm; 366 Q(0,0) = nu; 367 /* rdot[0] = r */ 368 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 369 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 370 } else { 371 /* select the current size of the subspace */ 372 if (l < ngmres->msize) { 373 l++; 374 } 375 k_restart++; 376 /* place the current entry in the list of previous entries */ 377 ierr = VecCopy(r, rdot[ivec]);CHKERRQ(ierr); 378 ierr = VecCopy(x, xdot[ivec]);CHKERRQ(ierr); 379 ngmres->r_norms[ivec] = r_norm; 380 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^A */ 381 for (i = 0; i < l; i++) { 382 VecDot(r, rdot[i], &qentry); 383 Q(i, ivec) = qentry; 384 Q(ivec, i) = qentry; 385 } 386 } 387 388 SNESLogConvHistory(snes, r_norm, k); 389 ierr = SNESMonitor(snes, k, r_norm);CHKERRQ(ierr); 390 391 snes->iter =k; 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 400 401 EXTERN_C_BEGIN 402 #undef __FUNCT__ 403 #define __FUNCT__ "SNESCreate_NGMRES" 404 PetscErrorCode SNESCreate_NGMRES(SNES snes) 405 { 406 SNES_NGMRES *ngmres; 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 snes->ops->destroy = SNESDestroy_NGMRES; 411 snes->ops->setup = SNESSetUp_NGMRES; 412 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 413 snes->ops->view = SNESView_NGMRES; 414 snes->ops->solve = SNESSolve_NGMRES; 415 snes->ops->reset = SNESReset_NGMRES; 416 417 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 418 snes->data = (void*) ngmres; 419 ngmres->msize = 10; 420 ngmres->debug = PETSC_FALSE; 421 422 ngmres->gammaA = 2.; 423 ngmres->gammaC = 2.; 424 ngmres->deltaB = 0.9; 425 ngmres->epsilonB = 0.1; 426 ngmres->k_rmax = 200; 427 428 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 EXTERN_C_END 432