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