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