Actual source code: qn.c

      1: #include <petsc/private/snesimpl.h>
      2: #include <petscdm.h>

      4: #define H(i,j)  qn->dXdFmat[i*qn->m + j]

      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};

     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;

     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;

     37:   /* basically just a regular newton's method except for the application of the Jacobian */

     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);

     42:   PetscCitationsRegister(SNESCitation,&SNEScite);
     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];

     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];

     53:   snes->reason = SNES_CONVERGED_ITERATING;

     55:   PetscObjectSAWsTakeAccess((PetscObject)snes);
     56:   snes->iter = 0;
     57:   snes->norm = 0.;
     58:   PetscObjectSAWsGrantAccess((PetscObject)snes);

     60:   if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
     61:     SNESApplyNPC(snes,X,NULL,F);
     62:     SNESGetConvergedReason(snes->npc,&reason);
     63:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
     64:       snes->reason = SNES_DIVERGED_INNER;
     65:       return(0);
     66:     }
     67:     VecNorm(F,NORM_2,&fnorm);
     68:   } else {
     69:     if (!snes->vec_func_init_set) {
     70:       SNESComputeFunction(snes,X,F);
     71:     } else snes->vec_func_init_set = PETSC_FALSE;

     73:     VecNorm(F,NORM_2,&fnorm);
     74:     SNESCheckFunctionNorm(snes,fnorm);
     75:   }
     76:   if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
     77:     SNESApplyNPC(snes,X,F,D);
     78:     SNESGetConvergedReason(snes->npc,&reason);
     79:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
     80:       snes->reason = SNES_DIVERGED_INNER;
     81:       return(0);
     82:     }
     83:   } else {
     84:     VecCopy(F,D);
     85:   }

     87:   PetscObjectSAWsTakeAccess((PetscObject)snes);
     88:   snes->norm = fnorm;
     89:   PetscObjectSAWsGrantAccess((PetscObject)snes);
     90:   SNESLogConvergenceHistory(snes,fnorm,0);
     91:   SNESMonitor(snes,0,fnorm);

     93:   /* test convergence */
     94:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
     95:   if (snes->reason) return(0);

     97:   if (snes->npc && snes->npcside== PC_RIGHT) {
     98:     PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
     99:     SNESSolve(snes->npc,snes->vec_rhs,X);
    100:     PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
    101:     SNESGetConvergedReason(snes->npc,&reason);
    102:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    103:       snes->reason = SNES_DIVERGED_INNER;
    104:       return(0);
    105:     }
    106:     SNESGetNPCFunction(snes,F,&fnorm);
    107:     VecCopy(F,D);
    108:   }

    110:   /* general purpose update */
    111:   if (snes->ops->update) {
    112:     (*snes->ops->update)(snes, snes->iter);
    113:   }

    115:   /* scale the initial update */
    116:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
    117:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
    118:     SNESCheckJacobianDomainerror(snes);
    119:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
    120:     MatLMVMSetJ0KSP(qn->B, snes->ksp);
    121:   }

    123:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
    124:     /* update QN approx and calculate step */
    125:     MatLMVMUpdate(qn->B, X, D);
    126:     MatSolve(qn->B, D, Y);

    128:     /* line search for lambda */
    129:     ynorm = 1; gnorm = fnorm;
    130:     VecCopy(D, Dold);
    131:     VecCopy(X, Xold);
    132:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
    133:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
    134:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
    135:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
    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:     }

    145:     /* convergence monitoring */
    146:     PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);

    148:     if (snes->npc && snes->npcside== PC_RIGHT) {
    149:       PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
    150:       SNESSolve(snes->npc,snes->vec_rhs,X);
    151:       PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
    152:       SNESGetConvergedReason(snes->npc,&reason);
    153:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
    154:         snes->reason = SNES_DIVERGED_INNER;
    155:         return(0);
    156:       }
    157:       SNESGetNPCFunction(snes,F,&fnorm);
    158:     }

    160:     SNESSetIterationNumber(snes, i+1);
    161:     snes->norm = fnorm;
    162:     snes->xnorm = xnorm;
    163:     snes->ynorm = ynorm;

    165:     SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
    166:     SNESMonitor(snes,snes->iter,snes->norm);

    168:     /* set parameter for default relative tolerance convergence test */
    169:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
    170:     if (snes->reason) return(0);
    171:     if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
    172:       SNESApplyNPC(snes,X,F,D);
    173:       SNESGetConvergedReason(snes->npc,&reason);
    174:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
    175:         snes->reason = SNES_DIVERGED_INNER;
    176:         return(0);
    177:       }
    178:     } else {
    179:       VecCopy(F, D);
    180:     }

    182:     /* general purpose update */
    183:     if (snes->ops->update) {
    184:       (*snes->ops->update)(snes, snes->iter);
    185:     }

    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:         MatMult(snes->jacobian_pre,Dold,W);
    193:       } else {
    194:         VecCopy(Dold,W);
    195:       }
    196:       VecDotBegin(W, Dold, &DolddotDold);
    197:       VecDotBegin(W, D, &DolddotD);
    198:       VecDotEnd(W, Dold, &DolddotDold);
    199:       VecDotEnd(W, D, &DolddotD);
    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:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
    210:         if (powell) {
    211:           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);
    212:         } else {
    213:           PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %D\n", i_r);
    214:         }
    215:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
    216:       }
    217:       i_r = -1;
    218:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
    219:         SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
    220:         SNESCheckJacobianDomainerror(snes);
    221:       }
    222:       MatLMVMReset(qn->B, PETSC_FALSE);
    223:     }
    224:   }
    225:   if (i == snes->max_its) {
    226:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
    227:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
    228:   }
    229:   return(0);
    230: }

    232: static PetscErrorCode SNESSetUp_QN(SNES snes)
    233: {
    234:   SNES_QN        *qn = (SNES_QN*)snes->data;
    236:   DM             dm;
    237:   PetscInt       n, N;

    241:   if (!snes->vec_sol) {
    242:     SNESGetDM(snes,&dm);
    243:     DMCreateGlobalVector(dm,&snes->vec_sol);
    244:   }
    245:   SNESSetWorkVecs(snes,4);

    247:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
    248:     SNESSetUpMatrices(snes);
    249:   }
    250:   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}

    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:   }

    268:   /* Set up the LMVM matrix */
    269:   switch (qn->type) {
    270:     case SNES_QN_BROYDEN:
    271:       MatSetType(qn->B, MATLMVMBROYDEN);
    272:       qn->scale_type = SNES_QN_SCALE_NONE;
    273:       break;
    274:     case SNES_QN_BADBROYDEN:
    275:       MatSetType(qn->B, MATLMVMBADBROYDEN);
    276:       qn->scale_type = SNES_QN_SCALE_NONE;
    277:       break;
    278:     default:
    279:       MatSetType(qn->B, MATLMVMBFGS);
    280:       switch (qn->scale_type) {
    281:         case SNES_QN_SCALE_NONE:
    282:           MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE);
    283:           break;
    284:         case SNES_QN_SCALE_SCALAR:
    285:           MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR);
    286:           break;
    287:         case SNES_QN_SCALE_JACOBIAN:
    288:           MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_USER);
    289:           break;
    290:         case SNES_QN_SCALE_DIAGONAL:
    291:         case SNES_QN_SCALE_DEFAULT:
    292:         default:
    293:           break;
    294:       }
    295:       break;
    296:   }
    297:   VecGetLocalSize(snes->vec_sol, &n);
    298:   VecGetSize(snes->vec_sol, &N);
    299:   MatSetSizes(qn->B, n, n, N, N);
    300:   MatSetUp(qn->B);
    301:   MatLMVMReset(qn->B, PETSC_TRUE);
    302:   MatLMVMSetHistorySize(qn->B, qn->m);
    303:   MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func);
    304:   return(0);
    305: }

    307: static PetscErrorCode SNESReset_QN(SNES snes)
    308: {
    310:   SNES_QN        *qn;

    313:   if (snes->data) {
    314:     qn = (SNES_QN*)snes->data;
    315:     MatDestroy(&qn->B);
    316:   }
    317:   return(0);
    318: }

    320: static PetscErrorCode SNESDestroy_QN(SNES snes)
    321: {

    325:   SNESReset_QN(snes);
    326:   PetscFree(snes->data);
    327:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
    328:   return(0);
    329: }

    331: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
    332: {

    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;

    343:   PetscOptionsHead(PetscOptionsObject,"SNES QN options");
    344:   PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);
    345:   PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance",          "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
    346:   PetscOptionsBool("-snes_qn_monitor",         "Monitor for the QN methods",      "SNESQN", qn->monflg, &qn->monflg, NULL);
    347:   PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
    348:   if (flg) SNESQNSetScaleType(snes,stype);

    350:   PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);
    351:   if (flg) SNESQNSetRestartType(snes,rtype);

    353:   PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);
    354:   if (flg) {SNESQNSetType(snes,qtype);}
    355:   MatSetFromOptions(qn->B);
    356:   PetscOptionsTail();
    357:   if (!snes->linesearch) {
    358:     SNESGetLineSearch(snes, &linesearch);
    359:     if (!((PetscObject)linesearch)->type_name) {
    360:       if (qn->type == SNES_QN_LBFGS) {
    361:         SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
    362:       } else if (qn->type == SNES_QN_BROYDEN) {
    363:         SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
    364:       } else {
    365:         SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
    366:       }
    367:     }
    368:   }
    369:   if (qn->monflg) {
    370:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor);
    371:   }
    372:   return(0);
    373: }

    375: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
    376: {
    377:   SNES_QN        *qn    = (SNES_QN*)snes->data;
    378:   PetscBool      iascii;

    382:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
    383:   if (iascii) {
    384:     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]);
    385:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %D\n", qn->m);
    386:   }
    387:   return(0);
    388: }

    390: /*@
    391:     SNESQNSetRestartType - Sets the restart type for SNESQN.

    393:     Logically Collective on SNES

    395:     Input Parameters:
    396: +   snes - the iterative context
    397: -   rtype - restart type

    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

    403:     Level: intermediate

    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

    410: @*/
    411: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
    412: {

    417:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
    418:   return(0);
    419: }

    421: /*@
    422:     SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN.

    424:     Logically Collective on SNES

    426:     Input Parameters:
    427: +   snes - the iterative context
    428: -   stype - scale type

    430:     Options Database:
    431: .   -snes_qn_scale_type <diagonal,none,scalar,jacobian>

    433:     Level: intermediate

    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.

    442: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESSetJacobian()
    443: @*/

    445: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
    446: {

    451:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
    452:   return(0);
    453: }

    455: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
    456: {
    457:   SNES_QN *qn = (SNES_QN*)snes->data;

    460:   qn->scale_type = stype;
    461:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
    462:   return(0);
    463: }

    465: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
    466: {
    467:   SNES_QN *qn = (SNES_QN*)snes->data;

    470:   qn->restart_type = rtype;
    471:   return(0);
    472: }

    474: /*@
    475:     SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN.

    477:     Logically Collective on SNES

    479:     Input Parameters:
    480: +   snes - the iterative context
    481: -   qtype - variant type

    483:     Options Database:
    484: .   -snes_qn_type <lbfgs,broyden,badbroyden>

    486:     Level: beginner

    488:     SNESQNTypes:
    489: +   SNES_QN_LBFGS - LBFGS variant
    490: .   SNES_QN_BROYDEN - Broyden variant
    491: -   SNES_QN_BADBROYDEN - Bad Broyden variant

    493: @*/

    495: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
    496: {

    501:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
    502:   return(0);
    503: }

    505: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
    506: {
    507:   SNES_QN *qn = (SNES_QN*)snes->data;

    510:   qn->type = qtype;
    511:   return(0);
    512: }

    514: /* -------------------------------------------------------------------------- */
    515: /*MC
    516:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

    518:       Options Database:

    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.

    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.

    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.

    539:       Uses left nonlinear preconditioning by default.

    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

    555:       Level: beginner

    557: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR

    559: M*/
    560: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
    561: {
    563:   SNES_QN        *qn;
    564:   const char     *optionsprefix;

    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;

    574:   snes->npcside= PC_LEFT;

    576:   snes->usesnpc = PETSC_TRUE;
    577:   snes->usesksp = PETSC_FALSE;

    579:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    581:   if (!snes->tolerancesset) {
    582:     snes->max_funcs = 30000;
    583:     snes->max_its   = 10000;
    584:   }

    586:   PetscNewLog(snes,&qn);
    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;

    597:   MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);
    598:   SNESGetOptionsPrefix(snes, &optionsprefix);
    599:   MatSetOptionsPrefix(qn->B, optionsprefix);

    601:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
    602:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
    603:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
    604:   return(0);
    605: }