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