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