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_",0}; 7 const char *const SNESQNRestartTypes[] = {"DEFAULT","NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0}; 8 const char *const SNESQNTypes[] = {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_",0}; 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, SNESetJacobian() 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 PetscFunctionReturn(0); 462 } 463 464 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 465 { 466 SNES_QN *qn = (SNES_QN*)snes->data; 467 468 PetscFunctionBegin; 469 qn->restart_type = rtype; 470 PetscFunctionReturn(0); 471 } 472 473 /*@ 474 SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN. 475 476 Logically Collective on SNES 477 478 Input Parameters: 479 + snes - the iterative context 480 - qtype - variant type 481 482 Options Database: 483 . -snes_qn_type <lbfgs,broyden,badbroyden> 484 485 Level: beginner 486 487 SNESQNTypes: 488 + SNES_QN_LBFGS - LBFGS variant 489 . SNES_QN_BROYDEN - Broyden variant 490 - SNES_QN_BADBROYDEN - Bad Broyden variant 491 492 @*/ 493 494 PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype) 495 { 496 PetscErrorCode ierr; 497 498 PetscFunctionBegin; 499 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 500 ierr = PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));CHKERRQ(ierr); 501 PetscFunctionReturn(0); 502 } 503 504 PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype) 505 { 506 SNES_QN *qn = (SNES_QN*)snes->data; 507 508 PetscFunctionBegin; 509 qn->type = qtype; 510 PetscFunctionReturn(0); 511 } 512 513 /* -------------------------------------------------------------------------- */ 514 /*MC 515 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 516 517 Options Database: 518 519 + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods. 520 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 521 . -snes_qn_powell_gamma - Angle condition for restart. 522 . -snes_qn_powell_descent - Descent condition for restart. 523 . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type 524 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian 525 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 526 - -snes_qn_monitor - Monitors the quasi-newton Jacobian. 527 528 Notes: 529 This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 530 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 531 updates. 532 533 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 534 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 535 iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed, 536 perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 537 538 Uses left nonlinear preconditioning by default. 539 540 References: 541 + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995. 542 . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods, 543 Technical Report, Northwestern University, June 1992. 544 . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE 545 Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 546 . 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 547 SIAM Review, 57(4), 2015 548 . 5. - Griewank, Andreas. "Broyden updating, the good and the bad!." Doc. Math (2012): 301-315. 549 . 6. - Gilbert, Jean Charles, and Claude Lemar{\'e}chal. "Some numerical experiments with variable-storage quasi-Newton algorithms." 550 Mathematical programming 45.1-3 (1989): 407-435. 551 - 7. - Dener A., Munson T. "Accelerating Limited-Memory Quasi-Newton Convergence for Large-Scale Optimization" 552 Computational Science - ICCS 2019. ICCS 2019. Lecture Notes in Computer Science, vol 11538. Springer, Cham 553 554 Level: beginner 555 556 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR 557 558 M*/ 559 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes) 560 { 561 PetscErrorCode ierr; 562 SNES_QN *qn; 563 const char *optionsprefix; 564 565 PetscFunctionBegin; 566 snes->ops->setup = SNESSetUp_QN; 567 snes->ops->solve = SNESSolve_QN; 568 snes->ops->destroy = SNESDestroy_QN; 569 snes->ops->setfromoptions = SNESSetFromOptions_QN; 570 snes->ops->view = SNESView_QN; 571 snes->ops->reset = SNESReset_QN; 572 573 snes->npcside= PC_LEFT; 574 575 snes->usesnpc = PETSC_TRUE; 576 snes->usesksp = PETSC_FALSE; 577 578 snes->alwayscomputesfinalresidual = PETSC_TRUE; 579 580 if (!snes->tolerancesset) { 581 snes->max_funcs = 30000; 582 snes->max_its = 10000; 583 } 584 585 ierr = PetscNewLog(snes,&qn);CHKERRQ(ierr); 586 snes->data = (void*) qn; 587 qn->m = 10; 588 qn->scaling = 1.0; 589 qn->monitor = NULL; 590 qn->monflg = PETSC_FALSE; 591 qn->powell_gamma = 0.9999; 592 qn->scale_type = SNES_QN_SCALE_DEFAULT; 593 qn->restart_type = SNES_QN_RESTART_DEFAULT; 594 qn->type = SNES_QN_LBFGS; 595 596 ierr = MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);CHKERRQ(ierr); 597 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 598 ierr = MatSetOptionsPrefix(qn->B, optionsprefix);CHKERRQ(ierr); 599 600 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr); 601 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr); 602 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605