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