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 (qn->type == SNES_QN_LBFGS) { 360 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 361 } else if (qn->type == SNES_QN_BROYDEN) { 362 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);CHKERRQ(ierr); 363 } else { 364 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 365 } 366 } 367 if (qn->monflg) { 368 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor);CHKERRQ(ierr); 369 } 370 PetscFunctionReturn(0); 371 } 372 373 static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer) 374 { 375 SNES_QN *qn = (SNES_QN*)snes->data; 376 PetscBool iascii; 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 381 if (iascii) { 382 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); 383 ierr = PetscViewerASCIIPrintf(viewer," Stored subspace size: %D\n", qn->m);CHKERRQ(ierr); 384 } 385 PetscFunctionReturn(0); 386 } 387 388 /*@ 389 SNESQNSetRestartType - Sets the restart type for SNESQN. 390 391 Logically Collective on SNES 392 393 Input Parameters: 394 + snes - the iterative context 395 - rtype - restart type 396 397 Options Database: 398 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 399 - -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic 400 401 Level: intermediate 402 403 SNESQNRestartTypes: 404 + SNES_QN_RESTART_NONE - never restart 405 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 406 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 407 408 @*/ 409 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) 410 { 411 PetscErrorCode ierr; 412 413 PetscFunctionBegin; 414 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 415 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 416 PetscFunctionReturn(0); 417 } 418 419 /*@ 420 SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN. 421 422 Logically Collective on SNES 423 424 Input Parameters: 425 + snes - the iterative context 426 - stype - scale type 427 428 Options Database: 429 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> 430 431 Level: intermediate 432 433 SNESQNScaleTypes: 434 + SNES_QN_SCALE_NONE - don't scale the problem 435 . SNES_QN_SCALE_SCALAR - use shanno scaling 436 . SNES_QN_SCALE_DIAGONAL - scale with a diagonalized BFGS formula (see Gilbert and Lemarechal 1989), available 437 - SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration 438 of QN and at ever restart. 439 440 .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian() 441 @*/ 442 443 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) 444 { 445 PetscErrorCode ierr; 446 447 PetscFunctionBegin; 448 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 449 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 450 PetscFunctionReturn(0); 451 } 452 453 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) 454 { 455 SNES_QN *qn = (SNES_QN*)snes->data; 456 457 PetscFunctionBegin; 458 qn->scale_type = stype; 459 PetscFunctionReturn(0); 460 } 461 462 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 463 { 464 SNES_QN *qn = (SNES_QN*)snes->data; 465 466 PetscFunctionBegin; 467 qn->restart_type = rtype; 468 PetscFunctionReturn(0); 469 } 470 471 /*@ 472 SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN. 473 474 Logically Collective on SNES 475 476 Input Parameters: 477 + snes - the iterative context 478 - qtype - variant type 479 480 Options Database: 481 . -snes_qn_type <lbfgs,broyden,badbroyden> 482 483 Level: beginner 484 485 SNESQNTypes: 486 + SNES_QN_LBFGS - LBFGS variant 487 . SNES_QN_BROYDEN - Broyden variant 488 - SNES_QN_BADBROYDEN - Bad Broyden variant 489 490 @*/ 491 492 PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype) 493 { 494 PetscErrorCode ierr; 495 496 PetscFunctionBegin; 497 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 498 ierr = PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));CHKERRQ(ierr); 499 PetscFunctionReturn(0); 500 } 501 502 PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype) 503 { 504 SNES_QN *qn = (SNES_QN*)snes->data; 505 506 PetscFunctionBegin; 507 qn->type = qtype; 508 PetscFunctionReturn(0); 509 } 510 511 /* -------------------------------------------------------------------------- */ 512 /*MC 513 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 514 515 Options Database: 516 517 + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods. 518 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 519 . -snes_qn_powell_gamma - Angle condition for restart. 520 . -snes_qn_powell_descent - Descent condition for restart. 521 . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type 522 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian 523 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 524 - -snes_qn_monitor - Monitors the quasi-newton Jacobian. 525 526 Notes: 527 This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 528 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 529 updates. 530 531 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 532 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 533 iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed, 534 perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 535 536 Uses left nonlinear preconditioning by default. 537 538 References: 539 + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995. 540 . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods, 541 Technical Report, Northwestern University, June 1992. 542 . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE 543 Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 544 . 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 545 SIAM Review, 57(4), 2015 546 . 5. - Griewank, Andreas. "Broyden updating, the good and the bad!." Doc. Math (2012): 301-315. 547 . 6. - Gilbert, Jean Charles, and Claude Lemaréchal. "Some numerical experiments with variable-storage quasi-Newton algorithms." 548 Mathematical programming 45.1-3 (1989): 407-435. 549 - 7. - Dener A., Munson T. "Accelerating Limited-Memory Quasi-Newton Convergence for Large-Scale Optimization" 550 Computational Science – ICCS 2019. ICCS 2019. Lecture Notes in Computer Science, vol 11538. Springer, Cham 551 552 Level: beginner 553 554 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR 555 556 M*/ 557 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes) 558 { 559 PetscErrorCode ierr; 560 SNES_QN *qn; 561 const char *optionsprefix; 562 563 PetscFunctionBegin; 564 snes->ops->setup = SNESSetUp_QN; 565 snes->ops->solve = SNESSolve_QN; 566 snes->ops->destroy = SNESDestroy_QN; 567 snes->ops->setfromoptions = SNESSetFromOptions_QN; 568 snes->ops->view = SNESView_QN; 569 snes->ops->reset = SNESReset_QN; 570 571 snes->npcside= PC_LEFT; 572 573 snes->usesnpc = PETSC_TRUE; 574 snes->usesksp = PETSC_FALSE; 575 576 snes->alwayscomputesfinalresidual = PETSC_TRUE; 577 578 if (!snes->tolerancesset) { 579 snes->max_funcs = 30000; 580 snes->max_its = 10000; 581 } 582 583 ierr = PetscNewLog(snes,&qn);CHKERRQ(ierr); 584 snes->data = (void*) qn; 585 qn->m = 10; 586 qn->scaling = 1.0; 587 qn->monitor = NULL; 588 qn->monflg = PETSC_FALSE; 589 qn->powell_gamma = 0.9999; 590 qn->scale_type = SNES_QN_SCALE_DEFAULT; 591 qn->restart_type = SNES_QN_RESTART_DEFAULT; 592 qn->type = SNES_QN_LBFGS; 593 594 ierr = MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);CHKERRQ(ierr); 595 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 596 ierr = MatSetOptionsPrefix(qn->B, optionsprefix);CHKERRQ(ierr); 597 598 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr); 599 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr); 600 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603