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