1 #include <private/snesimpl.h> 2 3 typedef enum {SNES_QN_SEQUENTIAL, SNES_QN_COMPOSED} SNESQNCompositionType; 4 5 typedef struct { 6 Vec *dX; /* The change in X */ 7 Vec *dF; /* The change in F */ 8 PetscInt m; /* the number of kept previous steps */ 9 PetscScalar *alpha; 10 PetscScalar *beta; 11 PetscScalar *rho; 12 PetscViewer monitor; 13 PetscReal powell_gamma; /* Powell angle restart condition */ 14 PetscReal powell_downhill; /* Powell descent restart condition */ 15 PetscReal scaling; /* scaling of H0 */ 16 17 SNESQNCompositionType compositiontype; /* determine if the composition is done sequentially or as a composition */ 18 19 } SNES_QN; 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "LBGFSApplyJinv_Private" 23 PetscErrorCode LBGFSApplyJinv_Private(SNES snes, PetscInt it, Vec D, Vec Y) { 24 25 PetscErrorCode ierr; 26 27 SNES_QN *qn = (SNES_QN*)snes->data; 28 29 Vec *dX = qn->dX; 30 Vec *dF = qn->dF; 31 32 PetscScalar *alpha = qn->alpha; 33 PetscScalar *beta = qn->beta; 34 PetscScalar *rho = qn->rho; 35 36 PetscInt k, i; 37 PetscInt m = qn->m; 38 PetscScalar t; 39 PetscInt l = m; 40 41 PetscFunctionBegin; 42 43 ierr = VecCopy(D, Y);CHKERRQ(ierr); 44 45 if (it < m) l = it; 46 47 /* outward recursion starting at iteration k's update and working back */ 48 for (i = 0; i < l; i++) { 49 k = (it - i - 1) % l; 50 ierr = VecDot(dX[k], Y, &t);CHKERRQ(ierr); 51 alpha[k] = t*rho[k]; 52 if (qn->monitor) { 53 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 54 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha: %14.12e\n", it, k, PetscRealPart(alpha[k]));CHKERRQ(ierr); 55 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 56 } 57 ierr = VecAXPY(Y, -alpha[k], dF[k]);CHKERRQ(ierr); 58 } 59 60 ierr = VecScale(Y, qn->scaling);CHKERRQ(ierr); 61 62 /* inward recursion starting at the first update and working forward */ 63 for (i = 0; i < l; i++) { 64 k = (it + i - l) % l; 65 ierr = VecDot(dF[k], Y, &t);CHKERRQ(ierr); 66 beta[k] = rho[k]*t; 67 ierr = VecAXPY(Y, (alpha[k] - beta[k]), dX[k]); 68 if (qn->monitor) { 69 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 70 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha - beta: %14.12e\n", it, k, PetscRealPart(alpha[k] - beta[k]));CHKERRQ(ierr); 71 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 72 } 73 } 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "SNESSolve_QN" 79 static PetscErrorCode SNESSolve_QN(SNES snes) 80 { 81 82 PetscErrorCode ierr; 83 SNES_QN *qn = (SNES_QN*) snes->data; 84 85 Vec X, Xold; 86 Vec F, G, B; 87 Vec W, Y, FPC, D, Dold; 88 SNESConvergedReason reason; 89 PetscInt i, i_r, k; 90 91 PetscReal fnorm, xnorm = 0, ynorm, gnorm; 92 PetscInt m = qn->m; 93 PetscBool lssucceed, changed; 94 95 PetscScalar rhosc; 96 97 Vec *dX = qn->dX; 98 Vec *dF = qn->dF; 99 PetscScalar *rho = qn->rho; 100 PetscScalar DolddotD, DolddotDold, DdotD, YdotD, a; 101 102 /* basically just a regular newton's method except for the application of the jacobian */ 103 PetscFunctionBegin; 104 105 X = snes->vec_sol; /* solution vector */ 106 F = snes->vec_func; /* residual vector */ 107 Y = snes->vec_sol_update; /* search direction generated by J^-1D*/ 108 B = snes->vec_rhs; 109 G = snes->work[0]; 110 W = snes->work[1]; 111 Xold = snes->work[2]; 112 113 /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */ 114 D = snes->work[3]; 115 Dold = snes->work[4]; 116 117 snes->reason = SNES_CONVERGED_ITERATING; 118 119 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 120 snes->iter = 0; 121 snes->norm = 0.; 122 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 123 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 124 if (snes->domainerror) { 125 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 126 PetscFunctionReturn(0); 127 } 128 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 129 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 130 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 131 snes->norm = fnorm; 132 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 133 SNESLogConvHistory(snes,fnorm,0); 134 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 135 136 /* set parameter for default relative tolerance convergence test */ 137 snes->ttol = fnorm*snes->rtol; 138 /* test convergence */ 139 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 140 if (snes->reason) PetscFunctionReturn(0); 141 142 /* composed solve -- either sequential or composed */ 143 if (snes->pc) { 144 if (qn->compositiontype == SNES_QN_SEQUENTIAL) { 145 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 146 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 147 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 148 snes->reason = SNES_DIVERGED_INNER; 149 PetscFunctionReturn(0); 150 } 151 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 152 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 153 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 154 ierr = VecCopy(F, Y);CHKERRQ(ierr); 155 } else { 156 ierr = VecCopy(X, Y);CHKERRQ(ierr); 157 ierr = SNESSolve(snes->pc, B, Y);CHKERRQ(ierr); 158 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 159 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 160 snes->reason = SNES_DIVERGED_INNER; 161 PetscFunctionReturn(0); 162 } 163 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); 164 } 165 } else { 166 ierr = VecCopy(F, Y);CHKERRQ(ierr); 167 } 168 ierr = VecCopy(Y, D);CHKERRQ(ierr); 169 170 for(i = 0, i_r = 0; i < snes->max_its; i++, i_r++) { 171 /* line search for lambda */ 172 ynorm = 1; gnorm = fnorm; 173 ierr = VecCopy(D, Dold);CHKERRQ(ierr); 174 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 175 ierr = SNESLineSearchPreCheckApply(snes,X,Y,&changed);CHKERRQ(ierr); 176 ierr = SNESLineSearchApply(snes,X,F,Y,fnorm,xnorm,W,G,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 177 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 178 if (snes->domainerror) { 179 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 180 PetscFunctionReturn(0); 181 } 182 if (!lssucceed) { 183 if (++snes->numFailures >= snes->maxFailures) { 184 snes->reason = SNES_DIVERGED_LINE_SEARCH; 185 break; 186 } 187 } 188 189 /* Update function and solution vectors */ 190 fnorm = gnorm; 191 ierr = VecCopy(G,F);CHKERRQ(ierr); 192 ierr = VecCopy(W,X);CHKERRQ(ierr); 193 194 /* convergence monitoring */ 195 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);CHKERRQ(ierr); 196 197 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 198 ierr = SNESSetFunctionNorm(snes, fnorm);CHKERRQ(ierr); 199 200 SNESLogConvHistory(snes,snes->norm,snes->iter); 201 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 202 /* set parameter for default relative tolerance convergence test */ 203 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 204 if (snes->reason) PetscFunctionReturn(0); 205 206 207 if (snes->pc) { 208 if (qn->compositiontype == SNES_QN_SEQUENTIAL) { 209 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 210 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 211 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 212 snes->reason = SNES_DIVERGED_INNER; 213 PetscFunctionReturn(0); 214 } 215 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 216 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 217 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 218 ierr = VecCopy(F, D);CHKERRQ(ierr); 219 } else { 220 ierr = VecCopy(X, D);CHKERRQ(ierr); 221 ierr = SNESSolve(snes->pc, B, D);CHKERRQ(ierr); 222 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 223 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 224 snes->reason = SNES_DIVERGED_INNER; 225 PetscFunctionReturn(0); 226 } 227 ierr = VecAYPX(D,-1.0,X);CHKERRQ(ierr); 228 } 229 } else { 230 ierr = VecCopy(F, D);CHKERRQ(ierr); 231 } 232 233 /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */ 234 ierr = VecDot(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 235 ierr = VecDot(Dold, D, &DolddotD);CHKERRQ(ierr); 236 ierr = VecDot(D, D, &DdotD);CHKERRQ(ierr); 237 ierr = VecDot(Y, D, &YdotD);CHKERRQ(ierr); 238 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) { 239 if (qn->monitor) { 240 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 241 ierr = PetscViewerASCIIPrintf(qn->monitor, "restart! |%14.12e| > %4.2f*|%14.12e|\n", k, PetscRealPart(DolddotD), qn->powell_gamma, PetscRealPart(DolddotDold));CHKERRQ(ierr); 242 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 243 } 244 i_r = -1; 245 /* general purpose update */ 246 if (snes->ops->update) { 247 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 248 } 249 ierr = VecCopy(D, Y);CHKERRQ(ierr); 250 } else { 251 /* set the differences */ 252 k = i_r % m; 253 ierr = VecCopy(D, dF[k]);CHKERRQ(ierr); 254 ierr = VecAXPY(dF[k], -1.0, Dold);CHKERRQ(ierr); 255 ierr = VecCopy(X, dX[k]);CHKERRQ(ierr); 256 ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr); 257 ierr = VecDot(dX[k], dF[k], &rhosc);CHKERRQ(ierr); 258 259 /* set scaling to be shanno scaling */ 260 rho[k] = 1. / rhosc; 261 ierr = VecDot(dF[k], dF[k], &a);CHKERRQ(ierr); 262 qn->scaling = PetscRealPart(rhosc) / PetscRealPart(a); 263 264 /* general purpose update */ 265 if (snes->ops->update) { 266 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 267 } 268 /* apply the current iteration of the approximate jacobian in order to get the next search direction*/ 269 ierr = LBGFSApplyJinv_Private(snes, i_r+1, D, Y);CHKERRQ(ierr); 270 } 271 } 272 if (i == snes->max_its) { 273 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 274 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 275 } 276 PetscFunctionReturn(0); 277 } 278 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "SNESSetUp_QN" 282 static PetscErrorCode SNESSetUp_QN(SNES snes) 283 { 284 SNES_QN *qn = (SNES_QN*)snes->data; 285 PetscErrorCode ierr; 286 PetscFunctionBegin; 287 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);CHKERRQ(ierr); 288 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);CHKERRQ(ierr); 289 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, qn->m, PetscScalar, &qn->beta, qn->m, PetscScalar, &qn->rho);CHKERRQ(ierr); 290 ierr = SNESDefaultGetWork(snes,5);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "SNESReset_QN" 296 static PetscErrorCode SNESReset_QN(SNES snes) 297 { 298 PetscErrorCode ierr; 299 SNES_QN *qn; 300 PetscFunctionBegin; 301 if (snes->data) { 302 qn = (SNES_QN*)snes->data; 303 if (qn->dX) { 304 ierr = VecDestroyVecs(qn->m, &qn->dX);CHKERRQ(ierr); 305 } 306 if (qn->dF) { 307 ierr = VecDestroyVecs(qn->m, &qn->dF);CHKERRQ(ierr); 308 } 309 ierr = PetscFree3(qn->alpha, qn->beta, qn->rho);CHKERRQ(ierr); 310 } 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "SNESDestroy_QN" 316 static PetscErrorCode SNESDestroy_QN(SNES snes) 317 { 318 PetscErrorCode ierr; 319 PetscFunctionBegin; 320 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 321 ierr = PetscFree(snes->data);CHKERRQ(ierr); 322 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "SNESSetFromOptions_QN" 328 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 329 { 330 331 PetscErrorCode ierr; 332 SNES_QN *qn; 333 const char *compositions[] = {"sequential", "composed"}; 334 PetscInt indx = 0; 335 PetscBool flg; 336 PetscBool monflg = PETSC_FALSE; 337 PetscFunctionBegin; 338 339 qn = (SNES_QN*)snes->data; 340 341 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 342 ierr = PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-BFGS methods", "SNESQN", qn->m, &qn->m, PETSC_NULL);CHKERRQ(ierr); 343 ierr = PetscOptionsReal("-snes_qn_powell_gamma", "Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, PETSC_NULL);CHKERRQ(ierr); 344 ierr = PetscOptionsReal("-snes_qn_powell_downhill", "Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, PETSC_NULL);CHKERRQ(ierr); 345 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, PETSC_NULL);CHKERRQ(ierr); 346 ierr = PetscOptionsEList("-snes_qn_composition", "Composition type", "SNESQN",compositions,2,"sequential",&indx,&flg);CHKERRQ(ierr); 347 if (flg) { 348 switch (indx) { 349 case 0: qn->compositiontype = SNES_QN_SEQUENTIAL; 350 break; 351 case 1: qn->compositiontype = SNES_QN_COMPOSED; 352 break; 353 } 354 } 355 ierr = PetscOptionsTail();CHKERRQ(ierr); 356 if (monflg) { 357 qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 EXTERN_C_BEGIN 363 #undef __FUNCT__ 364 #define __FUNCT__ "SNESLineSearchSetType_QN" 365 PetscErrorCode SNESLineSearchSetType_QN(SNES snes, SNESLineSearchType type) 366 { 367 PetscErrorCode ierr; 368 PetscFunctionBegin; 369 370 switch (type) { 371 case SNES_LS_BASIC: 372 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 373 break; 374 case SNES_LS_BASIC_NONORMS: 375 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 376 break; 377 case SNES_LS_QUADRATIC: 378 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 379 break; 380 case SNES_LS_CRITICAL: 381 ierr = SNESLineSearchSet(snes,SNESLineSearchCriticalSecant,PETSC_NULL);CHKERRQ(ierr); 382 break; 383 default: 384 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type"); 385 break; 386 } 387 snes->ls_type = type; 388 PetscFunctionReturn(0); 389 } 390 EXTERN_C_END 391 392 393 /* -------------------------------------------------------------------------- */ 394 /*MC 395 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 396 397 Options Database: 398 399 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 400 . -snes_qn_powell_angle - Angle condition for restart. 401 . -snes_qn_powell_descent - Descent condition for restart. 402 . -snes_qn_composition <sequential, composed>- Type of composition. 403 . -snes_ls <basic, basicnonorms, quadratic, critical> - Type of line search. 404 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 405 406 Notes: This implements the L-BFGS algorithm for the solution of F(x) = b using previous change in F(x) and x to 407 form the approximate inverse Jacobian using a series of multiplicative rank-one updates. This will eventually be 408 generalized to implement several limited-memory Broyden methods. 409 410 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 411 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 412 iteration as the current iteration's values when constructing the approximate jacobian. The second, composed, 413 perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 414 415 References: 416 417 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 418 International Journal of Computer Mathematics, vol. 86, 2009. 419 420 421 Level: beginner 422 423 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 424 425 M*/ 426 EXTERN_C_BEGIN 427 #undef __FUNCT__ 428 #define __FUNCT__ "SNESCreate_QN" 429 PetscErrorCode SNESCreate_QN(SNES snes) 430 { 431 432 PetscErrorCode ierr; 433 SNES_QN *qn; 434 435 PetscFunctionBegin; 436 snes->ops->setup = SNESSetUp_QN; 437 snes->ops->solve = SNESSolve_QN; 438 snes->ops->destroy = SNESDestroy_QN; 439 snes->ops->setfromoptions = SNESSetFromOptions_QN; 440 snes->ops->view = 0; 441 snes->ops->reset = SNESReset_QN; 442 443 snes->usespc = PETSC_TRUE; 444 snes->usesksp = PETSC_FALSE; 445 446 snes->max_funcs = 30000; 447 snes->max_its = 10000; 448 449 ierr = PetscNewLog(snes,SNES_QN,&qn);CHKERRQ(ierr); 450 snes->data = (void *) qn; 451 qn->m = 30; 452 qn->scaling = 1.0; 453 qn->dX = PETSC_NULL; 454 qn->dF = PETSC_NULL; 455 qn->monitor = PETSC_NULL; 456 qn->powell_gamma = 0.9; 457 qn->powell_downhill = 0.2; 458 qn->compositiontype = SNES_QN_SEQUENTIAL; 459 460 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_QN",SNESLineSearchSetType_QN);CHKERRQ(ierr); 461 ierr = SNESLineSearchSetType(snes, SNES_LS_CRITICAL);CHKERRQ(ierr); 462 463 PetscFunctionReturn(0); 464 } 465 EXTERN_C_END 466