Actual source code: diagbrdn.c

      1: #include <petscdevice.h>
      2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

      4: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
      5: {
      6:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
      7:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

      9:   PetscFunctionBegin;
     10:   PetscCall(VecPointwiseMult(dX, ldb->invD, F));
     11:   PetscFunctionReturn(PETSC_SUCCESS);
     12: }

     14: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
     15: {
     16:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
     17:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

     19:   PetscFunctionBegin;
     20:   PetscCall(VecPointwiseDivide(Z, X, ldb->invD));
     21:   PetscFunctionReturn(PETSC_SUCCESS);
     22: }

     24: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
     25: {
     26:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
     27:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;
     28:   PetscInt      old_k, i, start;
     29:   PetscScalar   curvature, ytDy, sts, stDs, ytDs;
     30:   PetscReal     curvtol, sigma, yy_sum, ss_sum, ys_sum, denom, ytytmp;
     31:   PetscReal     stDsr, ytDyr;

     33:   PetscFunctionBegin;
     34:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
     35:   if (lmvm->prev_set) {
     36:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
     37:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
     38:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

     40:     /* Test if the updates can be accepted */
     41:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ytytmp));
     42:     if (ytytmp < lmvm->eps) curvtol = 0.0;
     43:     else curvtol = lmvm->eps * ytytmp;

     45:     /* Test the curvature for the update */
     46:     if (PetscRealPart(curvature) > curvtol) {
     47:       /* Update is good so we accept it */
     48:       old_k = lmvm->k;
     49:       PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
     50:       /* If we hit the memory limit, shift the yty and yts arrays */
     51:       if (old_k == lmvm->k) {
     52:         for (i = 0; i <= lmvm->k - 1; ++i) {
     53:           ldb->yty[i] = ldb->yty[i + 1];
     54:           ldb->yts[i] = ldb->yts[i + 1];
     55:           ldb->sts[i] = ldb->sts[i + 1];
     56:         }
     57:       }
     58:       /* Accept dot products into the history */
     59:       PetscCall(VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &sts));
     60:       ldb->yty[lmvm->k] = ytytmp;
     61:       ldb->yts[lmvm->k] = PetscRealPart(curvature);
     62:       ldb->sts[lmvm->k] = PetscRealPart(sts);
     63:       if (ldb->forward) {
     64:         /* We are doing diagonal scaling of the forward Hessian B */
     65:         /*  BFGS = DFP = inv(D); */
     66:         PetscCall(VecCopy(ldb->invD, ldb->invDnew));
     67:         PetscCall(VecReciprocal(ldb->invDnew));

     69:         /*  V = y*y */
     70:         PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]));

     72:         /*  W = inv(D)*s */
     73:         PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]));
     74:         PetscCall(VecDot(ldb->W, lmvm->S[lmvm->k], &stDs));

     76:         /*  Safeguard stDs */
     77:         stDs = PetscMax(PetscRealPart(stDs), ldb->tol);

     79:         if (1.0 != ldb->theta) {
     80:           /*  BFGS portion of the update */
     81:           /*  U = (inv(D)*s)*(inv(D)*s) */
     82:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));

     84:           /*  Assemble */
     85:           PetscCall(VecAXPBY(ldb->BFGS, -1.0 / stDs, 0.0, ldb->U));
     86:         }
     87:         if (0.0 != ldb->theta) {
     88:           /*  DFP portion of the update */
     89:           /*  U = inv(D)*s*y */
     90:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]));

     92:           /*  Assemble */
     93:           PetscCall(VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V));
     94:           PetscCall(VecAXPY(ldb->DFP, -2.0, ldb->U));
     95:         }

     97:         if (0.0 == ldb->theta) {
     98:           PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->BFGS));
     99:         } else if (1.0 == ldb->theta) {
    100:           PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->DFP));
    101:         } else {
    102:           /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
    103:           PetscCall(VecAXPBYPCZ(ldb->invDnew, 1.0 - ldb->theta, (ldb->theta) / ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP));
    104:         }

    106:         PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
    107:         /*  Obtain inverse and ensure positive definite */
    108:         PetscCall(VecReciprocal(ldb->invDnew));
    109:         PetscCall(VecAbs(ldb->invDnew));

    111:       } else {
    112:         /* Inverse Hessian update instead. */
    113:         PetscCall(VecCopy(ldb->invD, ldb->invDnew));

    115:         /*  V = s*s */
    116:         PetscCall(VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]));

    118:         /*  W = D*y */
    119:         PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]));
    120:         PetscCall(VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy));

    122:         /*  Safeguard ytDy */
    123:         ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);

    125:         if (1.0 != ldb->theta) {
    126:           /*  BFGS portion of the update */
    127:           /*  U = s*Dy */
    128:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]));

    130:           /*  Assemble */
    131:           PetscCall(VecAXPBY(ldb->BFGS, ytDy / ldb->yts[lmvm->k], 0.0, ldb->V));
    132:           PetscCall(VecAXPY(ldb->BFGS, -2.0, ldb->U));
    133:         }
    134:         if (0.0 != ldb->theta) {
    135:           /*  DFP portion of the update */

    137:           /*  U = (inv(D)*y)*(inv(D)*y) */
    138:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));

    140:           /*  Assemble */
    141:           PetscCall(VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U));
    142:         }

    144:         if (0.0 == ldb->theta) {
    145:           PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->BFGS));
    146:         } else if (1.0 == ldb->theta) {
    147:           PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->DFP));
    148:         } else {
    149:           /*  Broyden update U=(1-theta)*P + theta*Q */
    150:           PetscCall(VecAXPBYPCZ(ldb->invDnew, (1.0 - ldb->theta) / ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP));
    151:         }
    152:         PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
    153:         /*  Ensure positive definite */
    154:         PetscCall(VecAbs(ldb->invDnew));
    155:       }
    156:       if (ldb->sigma_hist > 0) {
    157:         /*  Start with re-scaling on the newly computed diagonal */
    158:         if (0.5 == ldb->beta) {
    159:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
    160:             PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew));
    161:             PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));

    163:             PetscCall(VecDot(ldb->V, lmvm->Y[0], &ytDy));
    164:             PetscCall(VecDot(ldb->W, lmvm->S[0], &stDs));

    166:             ss_sum = PetscRealPart(stDs);
    167:             yy_sum = PetscRealPart(ytDy);
    168:             ys_sum = ldb->yts[0];
    169:           } else {
    170:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
    171:             PetscCall(VecReciprocal(ldb->U));

    173:             /*  Compute summations for scalar scaling */
    174:             yy_sum = 0; /*  No safeguard required */
    175:             ys_sum = 0; /*  No safeguard required */
    176:             ss_sum = 0; /*  No safeguard required */
    177:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
    178:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    179:               PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U));
    180:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

    182:               PetscCall(VecDot(ldb->W, lmvm->S[i], &stDs));
    183:               PetscCall(VecDot(ldb->V, lmvm->Y[i], &ytDy));

    185:               ss_sum += PetscRealPart(stDs);
    186:               ys_sum += ldb->yts[i];
    187:               yy_sum += PetscRealPart(ytDy);
    188:             }
    189:           }
    190:         } else if (0.0 == ldb->beta) {
    191:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
    192:             /*  Compute summations for scalar scaling */
    193:             PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));

    195:             PetscCall(VecDotNorm2(lmvm->Y[0], ldb->W, &ytDs, &stDsr));

    197:             ys_sum = PetscRealPart(ytDs);
    198:             ss_sum = stDsr;
    199:             yy_sum = ldb->yty[0];
    200:           } else {
    201:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
    202:             PetscCall(VecReciprocal(ldb->U));

    204:             /*  Compute summations for scalar scaling */
    205:             yy_sum = 0; /*  No safeguard required */
    206:             ys_sum = 0; /*  No safeguard required */
    207:             ss_sum = 0; /*  No safeguard required */
    208:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
    209:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    210:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

    212:               PetscCall(VecDotNorm2(lmvm->Y[i], ldb->W, &ytDs, &stDsr));

    214:               ss_sum += stDsr;
    215:               ys_sum += PetscRealPart(ytDs);
    216:               yy_sum += ldb->yty[i];
    217:             }
    218:           }
    219:         } else if (1.0 == ldb->beta) {
    220:           /*  Compute summations for scalar scaling */
    221:           yy_sum = 0; /*  No safeguard required */
    222:           ys_sum = 0; /*  No safeguard required */
    223:           ss_sum = 0; /*  No safeguard required */
    224:           start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
    225:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    226:             PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew));

    228:             PetscCall(VecDotNorm2(lmvm->S[i], ldb->V, &ytDs, &ytDyr));

    230:             yy_sum += ytDyr;
    231:             ys_sum += PetscRealPart(ytDs);
    232:             ss_sum += ldb->sts[i];
    233:           }
    234:         } else {
    235:           PetscCall(VecCopy(ldb->invDnew, ldb->U));
    236:           PetscCall(VecPow(ldb->U, ldb->beta - 1));

    238:           /*  Compute summations for scalar scaling */
    239:           yy_sum = 0; /*  No safeguard required */
    240:           ys_sum = 0; /*  No safeguard required */
    241:           ss_sum = 0; /*  No safeguard required */
    242:           start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
    243:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    244:             PetscCall(VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]));
    245:             PetscCall(VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]));

    247:             PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs, &ytDyr));
    248:             PetscCall(VecDot(ldb->W, ldb->W, &stDs));

    250:             yy_sum += ytDyr;
    251:             ys_sum += PetscRealPart(ytDs);
    252:             ss_sum += PetscRealPart(stDs);
    253:           }
    254:         }

    256:         if (0.0 == ldb->alpha) {
    257:           /*  Safeguard ys_sum  */
    258:           ys_sum = PetscMax(ldb->tol, ys_sum);

    260:           sigma = ss_sum / ys_sum;
    261:         } else if (1.0 == ldb->alpha) {
    262:           /* yy_sum is never 0; if it were, we'd be at the minimum */
    263:           sigma = ys_sum / yy_sum;
    264:         } else {
    265:           denom = 2.0 * ldb->alpha * yy_sum;

    267:           /*  Safeguard denom */
    268:           denom = PetscMax(ldb->tol, denom);

    270:           sigma = ((2.0 * ldb->alpha - 1) * ys_sum + PetscSqrtReal((2.0 * ldb->alpha - 1) * (2.0 * ldb->alpha - 1) * ys_sum * ys_sum - 4.0 * ldb->alpha * (ldb->alpha - 1) * yy_sum * ss_sum)) / denom;
    271:         }
    272:       } else {
    273:         sigma = 1.0;
    274:       }
    275:       /*  If Q has small values, then Q^(r_beta - 1)
    276:        can have very large values.  Hence, ys_sum
    277:        and ss_sum can be infinity.  In this case,
    278:        sigma can either be not-a-number or infinity. */

    280:       if (PetscIsInfOrNanScalar(sigma)) {
    281:         /*  sigma is not-a-number; skip rescaling */
    282:       } else if (0.0 == sigma) {
    283:         /*  sigma is zero; this is a bad case; skip rescaling */
    284:       } else {
    285:         /*  sigma is positive */
    286:         PetscCall(VecScale(ldb->invDnew, sigma));
    287:       }

    289:       /* Combine the old diagonal and the new diagonal using a convex limiter */
    290:       if (1.0 == ldb->rho) {
    291:         PetscCall(VecCopy(ldb->invDnew, ldb->invD));
    292:       } else if (ldb->rho) PetscCall(VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
    293:     } else {
    294:       PetscCall(MatLMVMReset(B, PETSC_FALSE));
    295:     }
    296:     /* End DiagBrdn update */
    297:   }
    298:   /* Save the solution and function to be used in the next update */
    299:   PetscCall(VecCopy(X, lmvm->Xprev));
    300:   PetscCall(VecCopy(F, lmvm->Fprev));
    301:   lmvm->prev_set = PETSC_TRUE;
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
    306: {
    307:   Mat_LMVM     *bdata = (Mat_LMVM *)B->data;
    308:   Mat_DiagBrdn *bctx  = (Mat_DiagBrdn *)bdata->ctx;
    309:   Mat_LMVM     *mdata = (Mat_LMVM *)M->data;
    310:   Mat_DiagBrdn *mctx  = (Mat_DiagBrdn *)mdata->ctx;
    311:   PetscInt      i;

    313:   PetscFunctionBegin;
    314:   mctx->theta      = bctx->theta;
    315:   mctx->alpha      = bctx->alpha;
    316:   mctx->beta       = bctx->beta;
    317:   mctx->rho        = bctx->rho;
    318:   mctx->delta      = bctx->delta;
    319:   mctx->delta_min  = bctx->delta_min;
    320:   mctx->delta_max  = bctx->delta_max;
    321:   mctx->tol        = bctx->tol;
    322:   mctx->sigma      = bctx->sigma;
    323:   mctx->sigma_hist = bctx->sigma_hist;
    324:   mctx->forward    = bctx->forward;
    325:   PetscCall(VecCopy(bctx->invD, mctx->invD));
    326:   for (i = 0; i <= bdata->k; ++i) {
    327:     mctx->yty[i] = bctx->yty[i];
    328:     mctx->yts[i] = bctx->yts[i];
    329:     mctx->sts[i] = bctx->sts[i];
    330:   }
    331:   PetscFunctionReturn(PETSC_SUCCESS);
    332: }

    334: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
    335: {
    336:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    337:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;
    338:   PetscBool     isascii;

    340:   PetscFunctionBegin;
    341:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
    342:   if (isascii) {
    343:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
    344:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
    345:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta));
    346:   }
    347:   PetscCall(MatView_LMVM(B, pv));
    348:   PetscFunctionReturn(PETSC_SUCCESS);
    349: }

    351: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
    352: {
    353:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    354:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

    356:   PetscFunctionBegin;
    357:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
    358:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
    359:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL, 0.0, 1.0));
    360:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL, 0.0, 1.0));
    361:   PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
    362:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL, 0.0, 1.0));
    363:   PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
    364:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
    365:   PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL, 0));
    366:   PetscOptionsHeadEnd();
    367:   PetscFunctionReturn(PETSC_SUCCESS);
    368: }

    370: static PetscErrorCode MatSetUp_DiagBrdn(Mat);
    371: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
    372: {
    373:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    374:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

    376:   PetscFunctionBegin;
    377:   if (!ldb->allocated) PetscCall(MatSetUp_DiagBrdn(B));
    378:   PetscCall(VecSet(ldb->invD, ldb->delta));
    379:   if (destructive && ldb->allocated) {
    380:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
    381:     PetscCall(VecDestroy(&ldb->invDnew));
    382:     PetscCall(VecDestroy(&ldb->invD));
    383:     PetscCall(VecDestroy(&ldb->BFGS));
    384:     PetscCall(VecDestroy(&ldb->DFP));
    385:     PetscCall(VecDestroy(&ldb->U));
    386:     PetscCall(VecDestroy(&ldb->V));
    387:     PetscCall(VecDestroy(&ldb->W));
    388:     ldb->allocated = PETSC_FALSE;
    389:   }
    390:   PetscCall(MatReset_LMVM(B, destructive));
    391:   PetscFunctionReturn(PETSC_SUCCESS);
    392: }

    394: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
    395: {
    396:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    397:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

    399:   PetscFunctionBegin;
    400:   PetscCall(MatAllocate_LMVM(B, X, F));
    401:   if (!ldb->allocated) {
    402:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
    403:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
    404:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
    405:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
    406:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
    407:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
    408:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
    409:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
    410:     ldb->allocated = PETSC_TRUE;
    411:   }
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

    415: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
    416: {
    417:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    418:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

    420:   PetscFunctionBegin;
    421:   if (ldb->allocated) {
    422:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
    423:     PetscCall(VecDestroy(&ldb->invDnew));
    424:     PetscCall(VecDestroy(&ldb->invD));
    425:     PetscCall(VecDestroy(&ldb->BFGS));
    426:     PetscCall(VecDestroy(&ldb->DFP));
    427:     PetscCall(VecDestroy(&ldb->U));
    428:     PetscCall(VecDestroy(&ldb->V));
    429:     PetscCall(VecDestroy(&ldb->W));
    430:     ldb->allocated = PETSC_FALSE;
    431:   }
    432:   PetscCall(PetscFree(lmvm->ctx));
    433:   PetscCall(MatDestroy_LMVM(B));
    434:   PetscFunctionReturn(PETSC_SUCCESS);
    435: }

    437: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
    438: {
    439:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    440:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

    442:   PetscFunctionBegin;
    443:   PetscCall(MatSetUp_LMVM(B));
    444:   if (!ldb->allocated) {
    445:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
    446:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
    447:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
    448:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
    449:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
    450:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
    451:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
    452:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
    453:     PetscCall(VecSet(ldb->invD, ldb->delta));
    454:     ldb->allocated = PETSC_TRUE;
    455:   }
    456:   PetscFunctionReturn(PETSC_SUCCESS);
    457: }

    459: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
    460: {
    461:   Mat_LMVM     *lmvm;
    462:   Mat_DiagBrdn *ldb;

    464:   PetscFunctionBegin;
    465:   PetscCall(MatCreate_LMVM(B));
    466:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN));
    467:   B->ops->setup          = MatSetUp_DiagBrdn;
    468:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
    469:   B->ops->destroy        = MatDestroy_DiagBrdn;
    470:   B->ops->view           = MatView_DiagBrdn;

    472:   lmvm                = (Mat_LMVM *)B->data;
    473:   lmvm->square        = PETSC_TRUE;
    474:   lmvm->m             = 1;
    475:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
    476:   lmvm->ops->reset    = MatReset_DiagBrdn;
    477:   lmvm->ops->mult     = MatMult_DiagBrdn;
    478:   lmvm->ops->solve    = MatSolve_DiagBrdn;
    479:   lmvm->ops->update   = MatUpdate_DiagBrdn;
    480:   lmvm->ops->copy     = MatCopy_DiagBrdn;

    482:   PetscCall(PetscNew(&ldb));
    483:   lmvm->ctx       = (void *)ldb;
    484:   ldb->theta      = 0.0;
    485:   ldb->alpha      = 1.0;
    486:   ldb->rho        = 1.0;
    487:   ldb->forward    = PETSC_TRUE;
    488:   ldb->beta       = 0.5;
    489:   ldb->sigma      = 1.0;
    490:   ldb->delta      = 1.0;
    491:   ldb->delta_min  = 1e-7;
    492:   ldb->delta_max  = 100.0;
    493:   ldb->tol        = 1e-8;
    494:   ldb->sigma_hist = 1;
    495:   ldb->allocated  = PETSC_FALSE;
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }

    499: /*@
    500:   MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
    501:   for approximating Hessians.

    503:   Collective

    505:   Input Parameters:
    506: + comm - MPI communicator
    507: . n    - number of local rows for storage vectors
    508: - N    - global size of the storage vectors

    510:   Output Parameter:
    511: . B - the matrix

    513:   Options Database Keys:
    514: + -mat_lmvm_theta      - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    515: . -mat_lmvm_rho        - (developer) update limiter for the J0 scaling
    516: . -mat_lmvm_alpha      - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    517: . -mat_lmvm_beta       - (developer) exponential factor for the diagonal J0 scaling
    518: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
    519: . -mat_lmvm_tol        - (developer) tolerance for bounding the denominator of the rescaling away from 0.
    520: - -mat_lmvm_forward    - (developer) whether or not to use the forward or backward Broyden update to the diagonal

    522:   Level: intermediate

    524:   Notes:
    525:   It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
    526:   paradigm instead of this routine directly.

    528:   It consists of a convex combination of DFP and BFGS
    529:   diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
    530:   To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
    531:   We also ensure positive definiteness by taking the `VecAbs()` of the final vector.

    533:   There are two ways of approximating the diagonal: using the forward (B) update
    534:   schemes for BFGS and DFP and then taking the inverse, or directly working with
    535:   the inverse (H) update schemes for the BFGS and DFP updates, derived using the
    536:   Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
    537:   parameter below.

    539:   In order to use the DiagBrdn matrix with other vector types, i.e. doing matrix-vector products
    540:   and matrix solves, the matrix must first be created using `MatCreate()` and `MatSetType()`,
    541:   followed by `MatLMVMAllocate()`. Then it will be available for updating
    542:   (via `MatLMVMUpdate()`) in one's favored solver implementation.

    544: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
    545:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
    546: @*/
    547: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    548: {
    549:   PetscFunctionBegin;
    550:   PetscCall(KSPInitializePackage());
    551:   PetscCall(MatCreate(comm, B));
    552:   PetscCall(MatSetSizes(*B, n, n, N, N));
    553:   PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
    554:   PetscCall(MatSetUp(*B));
    555:   PetscFunctionReturn(PETSC_SUCCESS);
    556: }