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 /* scale the initial update */ 197 if (qn->scalingtype == SNES_QN_JACOBIANSCALE) { 198 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 199 } 200 201 for(i = 0, i_r = 0; i < snes->max_its; i++, i_r++) { 202 ierr = LBGFSApplyJinv_Private(snes, i_r, D, Y);CHKERRQ(ierr); 203 /* line search for lambda */ 204 ynorm = 1; gnorm = fnorm; 205 ierr = VecCopy(D, Dold);CHKERRQ(ierr); 206 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 207 ierr = SNESLineSearchPreCheckApply(snes,X,Y,&changed);CHKERRQ(ierr); 208 ierr = SNESLineSearchApply(snes,X,F,Y,fnorm,xnorm,W,G,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 209 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 210 if (snes->domainerror) { 211 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 212 PetscFunctionReturn(0); 213 } 214 if (!lssucceed) { 215 if (++snes->numFailures >= snes->maxFailures) { 216 snes->reason = SNES_DIVERGED_LINE_SEARCH; 217 break; 218 } 219 } 220 if (qn->scalingtype == SNES_QN_LSSCALE) { 221 qn->scaling = ynorm; 222 } 223 /* Update function and solution vectors */ 224 fnorm = gnorm; 225 ierr = VecCopy(G,F);CHKERRQ(ierr); 226 ierr = VecCopy(W,X);CHKERRQ(ierr); 227 228 /* convergence monitoring */ 229 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); 230 231 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 232 ierr = SNESSetFunctionNorm(snes, fnorm);CHKERRQ(ierr); 233 234 SNESLogConvHistory(snes,snes->norm,snes->iter); 235 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 236 /* set parameter for default relative tolerance convergence test */ 237 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 238 if (snes->reason) PetscFunctionReturn(0); 239 240 241 if (snes->pc) { 242 if (qn->compositiontype == SNES_QN_SEQUENTIAL) { 243 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 244 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 245 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 246 snes->reason = SNES_DIVERGED_INNER; 247 PetscFunctionReturn(0); 248 } 249 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 250 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 251 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 252 ierr = VecCopy(F, D);CHKERRQ(ierr); 253 } else { 254 ierr = VecCopy(X, D);CHKERRQ(ierr); 255 ierr = SNESSolve(snes->pc, B, D);CHKERRQ(ierr); 256 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 257 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 258 snes->reason = SNES_DIVERGED_INNER; 259 PetscFunctionReturn(0); 260 } 261 ierr = VecAYPX(D,-1.0,X);CHKERRQ(ierr); 262 } 263 } else { 264 ierr = VecCopy(F, D);CHKERRQ(ierr); 265 } 266 267 /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */ 268 ierr = VecDot(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 269 ierr = VecDot(Dold, D, &DolddotD);CHKERRQ(ierr); 270 ierr = VecDot(D, D, &DdotD);CHKERRQ(ierr); 271 ierr = VecDot(Y, D, &YdotD);CHKERRQ(ierr); 272 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold)) || (i_r > qn->n_restart - 1 && qn->n_restart > 0)) { 273 if (qn->monitor) { 274 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 275 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); 276 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 277 } 278 i_r = -1; 279 /* general purpose update */ 280 if (snes->ops->update) { 281 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 282 } 283 if (qn->scalingtype == SNES_QN_JACOBIANSCALE) { 284 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 285 } 286 } else { 287 /* set the differences */ 288 k = i_r % m; 289 ierr = VecCopy(D, dF[k]);CHKERRQ(ierr); 290 ierr = VecAXPY(dF[k], -1.0, Dold);CHKERRQ(ierr); 291 ierr = VecCopy(X, dX[k]);CHKERRQ(ierr); 292 ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr); 293 ierr = VecDot(dX[k], dF[k], &rhosc);CHKERRQ(ierr); 294 295 /* set scaling to be shanno scaling */ 296 rho[k] = 1. / rhosc; 297 if (qn->scalingtype == SNES_QN_SHANNOSCALE) { 298 ierr = VecDot(dF[k], dF[k], &a);CHKERRQ(ierr); 299 qn->scaling = PetscRealPart(rhosc) / PetscRealPart(a); 300 } 301 /* general purpose update */ 302 if (snes->ops->update) { 303 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 304 } 305 } 306 } 307 if (i == snes->max_its) { 308 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 309 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 310 } 311 PetscFunctionReturn(0); 312 } 313 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "SNESSetUp_QN" 317 static PetscErrorCode SNESSetUp_QN(SNES snes) 318 { 319 SNES_QN *qn = (SNES_QN*)snes->data; 320 PetscErrorCode ierr; 321 PetscFunctionBegin; 322 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);CHKERRQ(ierr); 323 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);CHKERRQ(ierr); 324 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, qn->m, PetscScalar, &qn->beta, qn->m, PetscScalar, &qn->rho);CHKERRQ(ierr); 325 ierr = SNESDefaultGetWork(snes,5);CHKERRQ(ierr); 326 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "SNESReset_QN" 332 static PetscErrorCode SNESReset_QN(SNES snes) 333 { 334 PetscErrorCode ierr; 335 SNES_QN *qn; 336 PetscFunctionBegin; 337 if (snes->data) { 338 qn = (SNES_QN*)snes->data; 339 if (qn->dX) { 340 ierr = VecDestroyVecs(qn->m, &qn->dX);CHKERRQ(ierr); 341 } 342 if (qn->dF) { 343 ierr = VecDestroyVecs(qn->m, &qn->dF);CHKERRQ(ierr); 344 } 345 ierr = PetscFree3(qn->alpha, qn->beta, qn->rho);CHKERRQ(ierr); 346 } 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "SNESDestroy_QN" 352 static PetscErrorCode SNESDestroy_QN(SNES snes) 353 { 354 PetscErrorCode ierr; 355 PetscFunctionBegin; 356 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 357 ierr = PetscFree(snes->data);CHKERRQ(ierr); 358 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "SNESSetFromOptions_QN" 364 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 365 { 366 367 PetscErrorCode ierr; 368 SNES_QN *qn; 369 const char *compositions[] = {"sequential", "composed"}; 370 const char *scalings[] = {"shanno", "ls", "jacobian"}; 371 PetscInt indx = 0; 372 PetscBool flg; 373 PetscBool monflg = PETSC_FALSE; 374 PetscFunctionBegin; 375 376 qn = (SNES_QN*)snes->data; 377 378 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 379 ierr = PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-BFGS methods", "SNESQN", qn->m, &qn->m, PETSC_NULL);CHKERRQ(ierr); 380 ierr = PetscOptionsInt("-snes_qn_restart", "Maximum number of iterations between restarts", "SNESQN", qn->n_restart, &qn->n_restart, PETSC_NULL);CHKERRQ(ierr); 381 ierr = PetscOptionsReal("-snes_qn_powell_gamma", "Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, PETSC_NULL);CHKERRQ(ierr); 382 ierr = PetscOptionsReal("-snes_qn_powell_downhill", "Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, PETSC_NULL);CHKERRQ(ierr); 383 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, PETSC_NULL);CHKERRQ(ierr); 384 ierr = PetscOptionsEList("-snes_qn_composition", "Composition type", "SNESQN",compositions,2,"sequential",&indx,&flg);CHKERRQ(ierr); 385 if (flg) { 386 switch (indx) { 387 case 0: qn->compositiontype = SNES_QN_SEQUENTIAL; 388 break; 389 case 1: qn->compositiontype = SNES_QN_COMPOSED; 390 break; 391 } 392 } 393 ierr = PetscOptionsEList("-snes_qn_scaling", "Scaling type", "SNESQN",scalings,3,"shanno",&indx,&flg);CHKERRQ(ierr); 394 if (flg) { 395 switch (indx) { 396 case 0: qn->scalingtype = SNES_QN_SHANNOSCALE; 397 break; 398 case 1: qn->scalingtype = SNES_QN_LSSCALE; 399 break; 400 case 2: qn->scalingtype = SNES_QN_JACOBIANSCALE; 401 snes->usesksp = PETSC_TRUE; 402 break; 403 } 404 } 405 406 ierr = PetscOptionsTail();CHKERRQ(ierr); 407 if (monflg) { 408 qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 409 } 410 PetscFunctionReturn(0); 411 } 412 413 EXTERN_C_BEGIN 414 #undef __FUNCT__ 415 #define __FUNCT__ "SNESLineSearchSetType_QN" 416 PetscErrorCode SNESLineSearchSetType_QN(SNES snes, SNESLineSearchType type) 417 { 418 PetscErrorCode ierr; 419 PetscFunctionBegin; 420 421 switch (type) { 422 case SNES_LS_BASIC: 423 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 424 break; 425 case SNES_LS_BASIC_NONORMS: 426 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 427 break; 428 case SNES_LS_QUADRATIC: 429 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 430 break; 431 case SNES_LS_CRITICAL: 432 ierr = SNESLineSearchSet(snes,SNESLineSearchCriticalSecant,PETSC_NULL);CHKERRQ(ierr); 433 break; 434 default: 435 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type"); 436 break; 437 } 438 snes->ls_type = type; 439 PetscFunctionReturn(0); 440 } 441 EXTERN_C_END 442 443 444 /* -------------------------------------------------------------------------- */ 445 /*MC 446 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 447 448 Options Database: 449 450 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 451 . -snes_qn_powell_angle - Angle condition for restart. 452 . -snes_qn_powell_descent - Descent condition for restart. 453 . -snes_qn_composition <sequential, composed>- Type of composition. 454 . -snes_ls <basic, basicnonorms, quadratic, critical> - Type of line search. 455 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 456 457 Notes: This implements the L-BFGS algorithm for the solution of F(x) = b using previous change in F(x) and x to 458 form the approximate inverse Jacobian using a series of multiplicative rank-one updates. This will eventually be 459 generalized to implement several limited-memory Broyden methods. 460 461 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 462 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 463 iteration as the current iteration's values when constructing the approximate jacobian. The second, composed, 464 perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 465 466 References: 467 468 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 469 International Journal of Computer Mathematics, vol. 86, 2009. 470 471 472 Level: beginner 473 474 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 475 476 M*/ 477 EXTERN_C_BEGIN 478 #undef __FUNCT__ 479 #define __FUNCT__ "SNESCreate_QN" 480 PetscErrorCode SNESCreate_QN(SNES snes) 481 { 482 483 PetscErrorCode ierr; 484 SNES_QN *qn; 485 486 PetscFunctionBegin; 487 snes->ops->setup = SNESSetUp_QN; 488 snes->ops->solve = SNESSolve_QN; 489 snes->ops->destroy = SNESDestroy_QN; 490 snes->ops->setfromoptions = SNESSetFromOptions_QN; 491 snes->ops->view = 0; 492 snes->ops->reset = SNESReset_QN; 493 494 snes->usespc = PETSC_TRUE; 495 snes->usesksp = PETSC_FALSE; 496 497 snes->max_funcs = 30000; 498 snes->max_its = 10000; 499 500 ierr = PetscNewLog(snes,SNES_QN,&qn);CHKERRQ(ierr); 501 snes->data = (void *) qn; 502 qn->m = 10; 503 qn->scaling = 1.0; 504 qn->dX = PETSC_NULL; 505 qn->dF = PETSC_NULL; 506 qn->monitor = PETSC_NULL; 507 qn->powell_gamma = 0.9; 508 qn->powell_downhill = 0.2; 509 qn->compositiontype = SNES_QN_SEQUENTIAL; 510 qn->scalingtype = SNES_QN_SHANNOSCALE; 511 qn->n_restart = -1; 512 513 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_QN",SNESLineSearchSetType_QN);CHKERRQ(ierr); 514 ierr = SNESLineSearchSetType(snes, SNES_LS_CRITICAL);CHKERRQ(ierr); 515 516 PetscFunctionReturn(0); 517 } 518 EXTERN_C_END 519