1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdm.h> 3 4 #define H(i,j) qn->dXdFmat[i*qn->m + j] 5 6 const char *const SNESQNScaleTypes[] = {"DEFAULT","NONE","SCALAR","DIAGONAL","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",NULL}; 7 const char *const SNESQNRestartTypes[] = {"DEFAULT","NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",NULL}; 8 const char *const SNESQNTypes[] = {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_",NULL}; 9 10 typedef struct { 11 Mat B; /* Quasi-Newton approximation Matrix (MATLMVM) */ 12 PetscInt m; /* The number of kept previous steps */ 13 PetscReal *lambda; /* The line search history of the method */ 14 PetscBool monflg; 15 PetscViewer monitor; 16 PetscReal powell_gamma; /* Powell angle restart condition */ 17 PetscReal scaling; /* scaling of H0 */ 18 SNESQNType type; /* the type of quasi-newton method used */ 19 SNESQNScaleType scale_type; /* the type of scaling used */ 20 SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */ 21 } SNES_QN; 22 23 static PetscErrorCode SNESSolve_QN(SNES snes) 24 { 25 PetscErrorCode ierr; 26 SNES_QN *qn = (SNES_QN*) snes->data; 27 Vec X,Xold; 28 Vec F,W; 29 Vec Y,D,Dold; 30 PetscInt i, i_r; 31 PetscReal fnorm,xnorm,ynorm,gnorm; 32 SNESLineSearchReason lssucceed; 33 PetscBool badstep,powell,periodic; 34 PetscScalar DolddotD,DolddotDold; 35 SNESConvergedReason reason; 36 37 /* basically just a regular newton's method except for the application of the Jacobian */ 38 39 PetscFunctionBegin; 40 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 41 42 ierr = PetscCitationsRegister(SNESCitation,&SNEScite);CHKERRQ(ierr); 43 F = snes->vec_func; /* residual vector */ 44 Y = snes->vec_sol_update; /* search direction generated by J^-1D*/ 45 W = snes->work[3]; 46 X = snes->vec_sol; /* solution vector */ 47 Xold = snes->work[0]; 48 49 /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */ 50 D = snes->work[1]; 51 Dold = snes->work[2]; 52 53 snes->reason = SNES_CONVERGED_ITERATING; 54 55 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 56 snes->iter = 0; 57 snes->norm = 0.; 58 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 59 60 if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 61 ierr = SNESApplyNPC(snes,X,NULL,F);CHKERRQ(ierr); 62 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 63 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 64 snes->reason = SNES_DIVERGED_INNER; 65 PetscFunctionReturn(0); 66 } 67 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); 68 } else { 69 if (!snes->vec_func_init_set) { 70 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 71 } else snes->vec_func_init_set = PETSC_FALSE; 72 73 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); 74 SNESCheckFunctionNorm(snes,fnorm); 75 } 76 if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 77 ierr = SNESApplyNPC(snes,X,F,D);CHKERRQ(ierr); 78 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 79 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 80 snes->reason = SNES_DIVERGED_INNER; 81 PetscFunctionReturn(0); 82 } 83 } else { 84 ierr = VecCopy(F,D);CHKERRQ(ierr); 85 } 86 87 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 88 snes->norm = fnorm; 89 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 90 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 91 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 92 93 /* test convergence */ 94 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 95 if (snes->reason) PetscFunctionReturn(0); 96 97 if (snes->npc && snes->npcside== PC_RIGHT) { 98 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 99 ierr = SNESSolve(snes->npc,snes->vec_rhs,X);CHKERRQ(ierr); 100 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 101 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 102 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 103 snes->reason = SNES_DIVERGED_INNER; 104 PetscFunctionReturn(0); 105 } 106 ierr = SNESGetNPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 107 ierr = VecCopy(F,D);CHKERRQ(ierr); 108 } 109 110 /* general purpose update */ 111 if (snes->ops->update) { 112 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 113 } 114 115 /* scale the initial update */ 116 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 117 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 118 SNESCheckJacobianDomainerror(snes); 119 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 120 ierr = MatLMVMSetJ0KSP(qn->B, snes->ksp);CHKERRQ(ierr); 121 } 122 123 for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) { 124 /* update QN approx and calculate step */ 125 ierr = MatLMVMUpdate(qn->B, X, D);CHKERRQ(ierr); 126 ierr = MatSolve(qn->B, D, Y);CHKERRQ(ierr); 127 128 /* line search for lambda */ 129 ynorm = 1; gnorm = fnorm; 130 ierr = VecCopy(D, Dold);CHKERRQ(ierr); 131 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 132 ierr = SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 133 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 134 ierr = SNESLineSearchGetReason(snes->linesearch, &lssucceed);CHKERRQ(ierr); 135 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 136 badstep = PETSC_FALSE; 137 if (lssucceed) { 138 if (++snes->numFailures >= snes->maxFailures) { 139 snes->reason = SNES_DIVERGED_LINE_SEARCH; 140 break; 141 } 142 badstep = PETSC_TRUE; 143 } 144 145 /* convergence monitoring */ 146 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); 147 148 if (snes->npc && snes->npcside== PC_RIGHT) { 149 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 150 ierr = SNESSolve(snes->npc,snes->vec_rhs,X);CHKERRQ(ierr); 151 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 152 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 153 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 154 snes->reason = SNES_DIVERGED_INNER; 155 PetscFunctionReturn(0); 156 } 157 ierr = SNESGetNPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 158 } 159 160 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 161 snes->norm = fnorm; 162 snes->xnorm = xnorm; 163 snes->ynorm = ynorm; 164 165 ierr = SNESLogConvergenceHistory(snes,snes->norm,snes->iter);CHKERRQ(ierr); 166 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 167 168 /* set parameter for default relative tolerance convergence test */ 169 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 170 if (snes->reason) PetscFunctionReturn(0); 171 if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 172 ierr = SNESApplyNPC(snes,X,F,D);CHKERRQ(ierr); 173 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 174 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 175 snes->reason = SNES_DIVERGED_INNER; 176 PetscFunctionReturn(0); 177 } 178 } else { 179 ierr = VecCopy(F, D);CHKERRQ(ierr); 180 } 181 182 /* general purpose update */ 183 if (snes->ops->update) { 184 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 185 } 186 187 /* restart conditions */ 188 powell = PETSC_FALSE; 189 if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) { 190 /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */ 191 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 192 ierr = MatMult(snes->jacobian_pre,Dold,W);CHKERRQ(ierr); 193 } else { 194 ierr = VecCopy(Dold,W);CHKERRQ(ierr); 195 } 196 ierr = VecDotBegin(W, Dold, &DolddotDold);CHKERRQ(ierr); 197 ierr = VecDotBegin(W, D, &DolddotD);CHKERRQ(ierr); 198 ierr = VecDotEnd(W, Dold, &DolddotDold);CHKERRQ(ierr); 199 ierr = VecDotEnd(W, D, &DolddotD);CHKERRQ(ierr); 200 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE; 201 } 202 periodic = PETSC_FALSE; 203 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) { 204 if (i_r>qn->m-1) periodic = PETSC_TRUE; 205 } 206 /* restart if either powell or periodic restart is satisfied. */ 207 if (badstep || powell || periodic) { 208 if (qn->monflg) { 209 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 210 if (powell) { 211 ierr = PetscViewerASCIIPrintf(qn->monitor, "Powell restart! |%14.12e| > %6.4f*|%14.12e| i_r = %D\n", (double)PetscRealPart(DolddotD), (double)qn->powell_gamma, (double)PetscRealPart(DolddotDold),i_r);CHKERRQ(ierr); 212 } else { 213 ierr = PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %D\n", i_r);CHKERRQ(ierr); 214 } 215 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 216 } 217 i_r = -1; 218 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 219 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 220 SNESCheckJacobianDomainerror(snes); 221 } 222 ierr = MatLMVMReset(qn->B, PETSC_FALSE);CHKERRQ(ierr); 223 } 224 } 225 if (i == snes->max_its) { 226 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 227 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 228 } 229 PetscFunctionReturn(0); 230 } 231 232 static PetscErrorCode SNESSetUp_QN(SNES snes) 233 { 234 SNES_QN *qn = (SNES_QN*)snes->data; 235 PetscErrorCode ierr; 236 DM dm; 237 PetscInt n, N; 238 239 PetscFunctionBegin; 240 241 if (!snes->vec_sol) { 242 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 243 ierr = DMCreateGlobalVector(dm,&snes->vec_sol);CHKERRQ(ierr); 244 } 245 ierr = SNESSetWorkVecs(snes,4);CHKERRQ(ierr); 246 247 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 248 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 249 } 250 if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;} 251 252 /* set method defaults */ 253 if (qn->scale_type == SNES_QN_SCALE_DEFAULT) { 254 if (qn->type == SNES_QN_BADBROYDEN) { 255 qn->scale_type = SNES_QN_SCALE_NONE; 256 } else { 257 qn->scale_type = SNES_QN_SCALE_SCALAR; 258 } 259 } 260 if (qn->restart_type == SNES_QN_RESTART_DEFAULT) { 261 if (qn->type == SNES_QN_LBFGS) { 262 qn->restart_type = SNES_QN_RESTART_POWELL; 263 } else { 264 qn->restart_type = SNES_QN_RESTART_PERIODIC; 265 } 266 } 267 268 /* Set up the LMVM matrix */ 269 switch (qn->type) { 270 case SNES_QN_BROYDEN: 271 ierr = MatSetType(qn->B, MATLMVMBROYDEN);CHKERRQ(ierr); 272 qn->scale_type = SNES_QN_SCALE_NONE; 273 break; 274 case SNES_QN_BADBROYDEN: 275 ierr = MatSetType(qn->B, MATLMVMBADBROYDEN);CHKERRQ(ierr); 276 qn->scale_type = SNES_QN_SCALE_NONE; 277 break; 278 default: 279 ierr = MatSetType(qn->B, MATLMVMBFGS);CHKERRQ(ierr); 280 switch (qn->scale_type) { 281 case SNES_QN_SCALE_NONE: 282 ierr = MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE);CHKERRQ(ierr); 283 break; 284 case SNES_QN_SCALE_SCALAR: 285 ierr = MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR);CHKERRQ(ierr); 286 break; 287 case SNES_QN_SCALE_JACOBIAN: 288 ierr = MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_USER);CHKERRQ(ierr); 289 break; 290 case SNES_QN_SCALE_DIAGONAL: 291 case SNES_QN_SCALE_DEFAULT: 292 default: 293 break; 294 } 295 break; 296 } 297 ierr = VecGetLocalSize(snes->vec_sol, &n);CHKERRQ(ierr); 298 ierr = VecGetSize(snes->vec_sol, &N);CHKERRQ(ierr); 299 ierr = MatSetSizes(qn->B, n, n, N, N);CHKERRQ(ierr); 300 ierr = MatSetUp(qn->B);CHKERRQ(ierr); 301 ierr = MatLMVMReset(qn->B, PETSC_TRUE);CHKERRQ(ierr); 302 ierr = MatLMVMSetHistorySize(qn->B, qn->m);CHKERRQ(ierr); 303 ierr = MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func);CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 static PetscErrorCode SNESReset_QN(SNES snes) 308 { 309 PetscErrorCode ierr; 310 SNES_QN *qn; 311 312 PetscFunctionBegin; 313 if (snes->data) { 314 qn = (SNES_QN*)snes->data; 315 ierr = MatDestroy(&qn->B);CHKERRQ(ierr); 316 } 317 PetscFunctionReturn(0); 318 } 319 320 static PetscErrorCode SNESDestroy_QN(SNES snes) 321 { 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 326 ierr = PetscFree(snes->data);CHKERRQ(ierr); 327 ierr = PetscObjectComposeFunction((PetscObject)snes,"",NULL);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes) 332 { 333 334 PetscErrorCode ierr; 335 SNES_QN *qn = (SNES_QN*)snes->data; 336 PetscBool flg; 337 SNESLineSearch linesearch; 338 SNESQNRestartType rtype = qn->restart_type; 339 SNESQNScaleType stype = qn->scale_type; 340 SNESQNType qtype = qn->type; 341 342 PetscFunctionBegin; 343 ierr = PetscOptionsHead(PetscOptionsObject,"SNES QN options");CHKERRQ(ierr); 344 ierr = PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);CHKERRQ(ierr); 345 ierr = PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);CHKERRQ(ierr); 346 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", qn->monflg, &qn->monflg, NULL);CHKERRQ(ierr); 347 ierr = PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);CHKERRQ(ierr); 348 if (flg) {ierr = SNESQNSetScaleType(snes,stype);CHKERRQ(ierr);} 349 350 ierr = PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);CHKERRQ(ierr); 351 if (flg) {ierr = SNESQNSetRestartType(snes,rtype);CHKERRQ(ierr);} 352 353 ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);CHKERRQ(ierr); 354 if (flg) {ierr = SNESQNSetType(snes,qtype);CHKERRQ(ierr);} 355 ierr = MatSetFromOptions(qn->B);CHKERRQ(ierr); 356 ierr = PetscOptionsTail();CHKERRQ(ierr); 357 if (!snes->linesearch) { 358 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 359 if (!((PetscObject)linesearch)->type_name) { 360 if (qn->type == SNES_QN_LBFGS) { 361 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 362 } else if (qn->type == SNES_QN_BROYDEN) { 363 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);CHKERRQ(ierr); 364 } else { 365 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 366 } 367 } 368 } 369 if (qn->monflg) { 370 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor);CHKERRQ(ierr); 371 } 372 PetscFunctionReturn(0); 373 } 374 375 static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer) 376 { 377 SNES_QN *qn = (SNES_QN*)snes->data; 378 PetscBool iascii; 379 PetscErrorCode ierr; 380 381 PetscFunctionBegin; 382 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 383 if (iascii) { 384 ierr = PetscViewerASCIIPrintf(viewer," type is %s, restart type is %s, scale type is %s\n",SNESQNTypes[qn->type],SNESQNRestartTypes[qn->restart_type],SNESQNScaleTypes[qn->scale_type]);CHKERRQ(ierr); 385 ierr = PetscViewerASCIIPrintf(viewer," Stored subspace size: %D\n", qn->m);CHKERRQ(ierr); 386 } 387 PetscFunctionReturn(0); 388 } 389 390 /*@ 391 SNESQNSetRestartType - Sets the restart type for SNESQN. 392 393 Logically Collective on SNES 394 395 Input Parameters: 396 + snes - the iterative context 397 - rtype - restart type 398 399 Options Database: 400 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 401 - -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic 402 403 Level: intermediate 404 405 SNESQNRestartTypes: 406 + SNES_QN_RESTART_NONE - never restart 407 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 408 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 409 410 @*/ 411 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) 412 { 413 PetscErrorCode ierr; 414 415 PetscFunctionBegin; 416 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 417 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 /*@ 422 SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN. 423 424 Logically Collective on SNES 425 426 Input Parameters: 427 + snes - the iterative context 428 - stype - scale type 429 430 Options Database: 431 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> 432 433 Level: intermediate 434 435 SNESQNScaleTypes: 436 + SNES_QN_SCALE_NONE - don't scale the problem 437 . SNES_QN_SCALE_SCALAR - use Shanno scaling 438 . SNES_QN_SCALE_DIAGONAL - scale with a diagonalized BFGS formula (see Gilbert and Lemarechal 1989), available 439 - SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration 440 of QN and at ever restart. 441 442 .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESSetJacobian() 443 @*/ 444 445 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) 446 { 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 451 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) 456 { 457 SNES_QN *qn = (SNES_QN*)snes->data; 458 459 PetscFunctionBegin; 460 qn->scale_type = stype; 461 if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE; 462 PetscFunctionReturn(0); 463 } 464 465 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 466 { 467 SNES_QN *qn = (SNES_QN*)snes->data; 468 469 PetscFunctionBegin; 470 qn->restart_type = rtype; 471 PetscFunctionReturn(0); 472 } 473 474 /*@ 475 SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN. 476 477 Logically Collective on SNES 478 479 Input Parameters: 480 + snes - the iterative context 481 - qtype - variant type 482 483 Options Database: 484 . -snes_qn_type <lbfgs,broyden,badbroyden> 485 486 Level: beginner 487 488 SNESQNTypes: 489 + SNES_QN_LBFGS - LBFGS variant 490 . SNES_QN_BROYDEN - Broyden variant 491 - SNES_QN_BADBROYDEN - Bad Broyden variant 492 493 @*/ 494 495 PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype) 496 { 497 PetscErrorCode ierr; 498 499 PetscFunctionBegin; 500 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 501 ierr = PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));CHKERRQ(ierr); 502 PetscFunctionReturn(0); 503 } 504 505 PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype) 506 { 507 SNES_QN *qn = (SNES_QN*)snes->data; 508 509 PetscFunctionBegin; 510 qn->type = qtype; 511 PetscFunctionReturn(0); 512 } 513 514 /* -------------------------------------------------------------------------- */ 515 /*MC 516 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 517 518 Options Database: 519 520 + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods. 521 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 522 . -snes_qn_powell_gamma - Angle condition for restart. 523 . -snes_qn_powell_descent - Descent condition for restart. 524 . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type 525 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian 526 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 527 - -snes_qn_monitor - Monitors the quasi-newton Jacobian. 528 529 Notes: 530 This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 531 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 532 updates. 533 534 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 535 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 536 iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed, 537 perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 538 539 Uses left nonlinear preconditioning by default. 540 541 References: 542 + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995. 543 . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods, 544 Technical Report, Northwestern University, June 1992. 545 . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE 546 Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 547 . 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 548 SIAM Review, 57(4), 2015 549 . 5. - Griewank, Andreas. "Broyden updating, the good and the bad!." Doc. Math (2012): 301-315. 550 . 6. - Gilbert, Jean Charles, and Claude Lemar{\'e}chal. "Some numerical experiments with variable-storage quasi-Newton algorithms." 551 Mathematical programming 45.1-3 (1989): 407-435. 552 - 7. - Dener A., Munson T. "Accelerating Limited-Memory Quasi-Newton Convergence for Large-Scale Optimization" 553 Computational Science - ICCS 2019. ICCS 2019. Lecture Notes in Computer Science, vol 11538. Springer, Cham 554 555 Level: beginner 556 557 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR 558 559 M*/ 560 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes) 561 { 562 PetscErrorCode ierr; 563 SNES_QN *qn; 564 const char *optionsprefix; 565 566 PetscFunctionBegin; 567 snes->ops->setup = SNESSetUp_QN; 568 snes->ops->solve = SNESSolve_QN; 569 snes->ops->destroy = SNESDestroy_QN; 570 snes->ops->setfromoptions = SNESSetFromOptions_QN; 571 snes->ops->view = SNESView_QN; 572 snes->ops->reset = SNESReset_QN; 573 574 snes->npcside= PC_LEFT; 575 576 snes->usesnpc = PETSC_TRUE; 577 snes->usesksp = PETSC_FALSE; 578 579 snes->alwayscomputesfinalresidual = PETSC_TRUE; 580 581 if (!snes->tolerancesset) { 582 snes->max_funcs = 30000; 583 snes->max_its = 10000; 584 } 585 586 ierr = PetscNewLog(snes,&qn);CHKERRQ(ierr); 587 snes->data = (void*) qn; 588 qn->m = 10; 589 qn->scaling = 1.0; 590 qn->monitor = NULL; 591 qn->monflg = PETSC_FALSE; 592 qn->powell_gamma = 0.9999; 593 qn->scale_type = SNES_QN_SCALE_DEFAULT; 594 qn->restart_type = SNES_QN_RESTART_DEFAULT; 595 qn->type = SNES_QN_LBFGS; 596 597 ierr = MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);CHKERRQ(ierr); 598 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 599 ierr = MatSetOptionsPrefix(qn->B, optionsprefix);CHKERRQ(ierr); 600 601 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr); 602 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr); 603 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606