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, 3);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, f, 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 f = snes->work[0]; 172 b = snes->vec_rhs; 173 x_A = snes->vec_sol_update; 174 r_A = snes->work[1]; 175 d = snes->work[2]; 176 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 177 178 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 179 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 180 snes->iter = 0; 181 snes->norm = 0.; 182 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 183 184 /* initialization */ 185 186 /* r = F(x) */ 187 ierr = SNESComputeFunction(snes, x, r);CHKERRQ(ierr); 188 if (snes->domainerror) { 189 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 190 PetscFunctionReturn(0); 191 } 192 193 /* nu = (r, r) */ 194 ierr = VecNorm(r, NORM_2, &r_norm);CHKERRQ(ierr); 195 r_min_norm = r_norm; 196 nu = r_norm*r_norm; 197 if (PetscIsInfOrNanReal(r_norm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 198 199 /* q_{00} = nu */ 200 Q(0,0) = nu; 201 ngmres->r_norms[0] = r_norm; 202 /* rdot[0] = r */ 203 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 204 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 205 206 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 207 snes->norm = r_norm; 208 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 209 SNESLogConvHistory(snes, r_norm, 0); 210 ierr = SNESMonitor(snes, 0, r_norm);CHKERRQ(ierr); 211 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 212 if (snes->reason) PetscFunctionReturn(0); 213 214 k_restart = 1; 215 l = 1; 216 for (k=1; k<snes->max_its; k++) { 217 218 /* select which vector of the stored subspace will be updated */ 219 ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */ 220 221 222 /* Computation of x^M */ 223 ierr = SNESSolve(pc, b, x);CHKERRQ(ierr); 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 VecDot(rdot[i], r, &xi[i]); 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 ngmres->nrhs = 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=10000000; 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=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->debug) 326 PetscPrintf(PETSC_COMM_WORLD, "picked r_A, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm); 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->debug) 334 PetscPrintf(PETSC_COMM_WORLD, "picked r_M, ||r_A||_2 = %e, ||r_M||_2 = %e\n", r_A_norm, r_norm); 335 } 336 337 selectRestart = PETSC_FALSE; 338 339 /* maximum iteration criterion */ 340 if (k_restart > ngmres->k_rmax) { 341 selectRestart = PETSC_TRUE; 342 } 343 344 /* difference stagnation restart */ 345 if ((ngmres->epsilonB*d_norm > d_min_norm) && 346 (sqrt(r_A_norm) > ngmres->deltaB*sqrt(r_min_norm))) { 347 if (ngmres->debug) 348 PetscPrintf(PETSC_COMM_WORLD, "difference restart: %e > %e\n", ngmres->epsilonB*d_norm, d_min_norm); 349 selectRestart = PETSC_TRUE; 350 } 351 352 /* residual stagnation restart */ 353 if (sqrt(r_A_norm) > ngmres->gammaC*sqrt(r_min_norm)) { 354 if (ngmres->debug) 355 PetscPrintf(PETSC_COMM_WORLD, "residual restart: %e > %e\n", sqrt(r_A_norm), ngmres->gammaC*sqrt(r_min_norm)); 356 selectRestart = PETSC_TRUE; 357 } 358 359 if (selectRestart) { 360 if (ngmres->debug) 361 PetscPrintf(PETSC_COMM_WORLD, "Restarted at iteration %d\n", k_restart); 362 k_restart = 1; 363 l = 1; 364 /* q_{00} = nu */ 365 ngmres->r_norms[0] = r_norm; 366 nu = r_norm*r_norm; 367 Q(0,0) = nu; 368 /* rdot[0] = r */ 369 ierr = VecCopy(x, xdot[0]);CHKERRQ(ierr); 370 ierr = VecCopy(r, rdot[0]);CHKERRQ(ierr); 371 } else { 372 /* select the current size of the subspace */ 373 if (l < ngmres->msize) { 374 l++; 375 } 376 k_restart++; 377 /* place the current entry in the list of previous entries */ 378 ierr = VecCopy(r, rdot[ivec]);CHKERRQ(ierr); 379 ierr = VecCopy(x, xdot[ivec]);CHKERRQ(ierr); 380 ngmres->r_norms[ivec] = r_norm; 381 if (r_min_norm > r_norm) r_min_norm = r_norm; /* the minimum norm is now of r^A */ 382 for (i = 0; i < l; i++) { 383 VecDot(r, rdot[i], &qentry); 384 Q(i, ivec) = qentry; 385 Q(ivec, i) = qentry; 386 } 387 } 388 389 SNESLogConvHistory(snes, r_norm, k); 390 ierr = SNESMonitor(snes, k, r_norm);CHKERRQ(ierr); 391 392 snes->iter =k; 393 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,r_norm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 394 if (snes->reason) PetscFunctionReturn(0); 395 } 396 snes->reason = SNES_DIVERGED_MAX_IT; 397 PetscFunctionReturn(0); 398 } 399 400 401 402 EXTERN_C_BEGIN 403 #undef __FUNCT__ 404 #define __FUNCT__ "SNESCreate_NGMRES" 405 PetscErrorCode SNESCreate_NGMRES(SNES snes) 406 { 407 SNES_NGMRES *ngmres; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 snes->ops->destroy = SNESDestroy_NGMRES; 412 snes->ops->setup = SNESSetUp_NGMRES; 413 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 414 snes->ops->view = SNESView_NGMRES; 415 snes->ops->solve = SNESSolve_NGMRES; 416 snes->ops->reset = SNESReset_NGMRES; 417 418 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 419 snes->data = (void*) ngmres; 420 ngmres->msize = 10; 421 ngmres->debug = PETSC_FALSE; 422 423 ngmres->gammaA = 2.; 424 ngmres->gammaC = 2.; 425 ngmres->deltaB = 0.9; 426 ngmres->epsilonB = 0.1; 427 ngmres->k_rmax = 200; 428 429 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 430 PetscFunctionReturn(0); 431 } 432 EXTERN_C_END 433