Actual source code: diagbrdn.c

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

      3: /*------------------------------------------------------------*/

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

     12:   VecCheckSameSize(F, 2, dX, 3);
     13:   VecCheckMatCompatible(B, dX, 3, F, 2);
     14:   VecPointwiseMult(dX, ldb->invD, F);
     15:   return(0);
     16: }

     18: /*------------------------------------------------------------*/

     20: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
     21: {
     22:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     23:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
     24:   PetscErrorCode    ierr;

     27:   VecCheckSameSize(X, 2, Z, 3);
     28:   VecCheckMatCompatible(B, X, 2, Z, 3);
     29:   VecPointwiseDivide(Z, X, ldb->invD);
     30:   return(0);
     31: }

     33: /*------------------------------------------------------------*/

     35: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
     36: {
     37:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
     38:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
     39:   PetscErrorCode    ierr;
     40:   PetscInt          old_k, i, start;
     41:   PetscScalar       yty, ststmp, curvature, ytDy, stDs, ytDs;
     42:   PetscReal         curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;

     45:   if (!lmvm->m) return(0);
     46:   if (lmvm->prev_set) {
     47:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
     48:     VecAYPX(lmvm->Xprev, -1.0, X);
     49:     VecAYPX(lmvm->Fprev, -1.0, F);
     50:     /* Compute tolerance for accepting the update */
     51:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
     52:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
     53:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
     54:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
     55:     if (PetscRealPart(ststmp) < lmvm->eps){
     56:       curvtol = 0.0;
     57:     } else {
     58:       curvtol = lmvm->eps * PetscRealPart(ststmp);
     59:     }
     60:     /* Test the curvature for the update */
     61:     if (PetscRealPart(curvature) > curvtol) {
     62:       /* Update is good so we accept it */
     63:       old_k = lmvm->k;
     64:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
     65:       /* If we hit the memory limit, shift the yty and yts arrays */
     66:       if (old_k == lmvm->k) {
     67:         for (i = 0; i <= lmvm->k-1; ++i) {
     68:           ldb->yty[i] = ldb->yty[i+1];
     69:           ldb->yts[i] = ldb->yts[i+1];
     70:           ldb->sts[i] = ldb->sts[i+1];
     71:         }
     72:       }
     73:       /* Accept dot products into the history */
     74:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
     75:       ldb->yty[lmvm->k] = PetscRealPart(yty);
     76:       ldb->yts[lmvm->k] = PetscRealPart(curvature);
     77:       ldb->sts[lmvm->k] = PetscRealPart(ststmp);
     78:       if (ldb->forward) {
     79:         /* We are doing diagonal scaling of the forward Hessian B */
     80:         /*  BFGS = DFP = inv(D); */
     81:         VecCopy(ldb->invD, ldb->invDnew);
     82:         VecReciprocal(ldb->invDnew);

     84:         /*  V = y*y */
     85:         VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);

     87:         /*  W = inv(D)*s */
     88:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
     89:         VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);

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

     94:         if (1.0 != ldb->theta) {
     95:           /*  BFGS portion of the update */
     96:           /*  U = (inv(D)*s)*(inv(D)*s) */
     97:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

     99:           /*  Assemble */
    100:           VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);
    101:         }
    102:         if (0.0 != ldb->theta) {
    103:           /*  DFP portion of the update */
    104:           /*  U = inv(D)*s*y */
    105:           VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);

    107:           /*  Assemble */
    108:           VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);
    109:           VecAXPY(ldb->DFP, -2.0, ldb->U);
    110:         }

    112:         if (0.0 == ldb->theta) {
    113:           VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
    114:         } else if (1.0 == ldb->theta) {
    115:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);
    116:         } else {
    117:           /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
    118:           VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
    119:         }

    121:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
    122:         /*  Obtain inverse and ensure positive definite */
    123:         VecReciprocal(ldb->invDnew);
    124:         VecAbs(ldb->invDnew);

    126:       } else {
    127:         /* Inverse Hessian update instead. */
    128:         VecCopy(ldb->invD, ldb->invDnew);

    130:         /*  V = s*s */
    131:         VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);

    133:         /*  W = D*y */
    134:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
    135:         VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);

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

    140:         if (1.0 != ldb->theta) {
    141:           /*  BFGS portion of the update */
    142:           /*  U = s*Dy */
    143:           VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);

    145:           /*  Assemble */
    146:           VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);
    147:           VecAXPY(ldb->BFGS, -2.0, ldb->U);
    148:         }
    149:         if (0.0 != ldb->theta) {
    150:           /*  DFP portion of the update */

    152:           /*  U = (inv(D)*y)*(inv(D)*y) */
    153:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

    155:           /*  Assemble */
    156:           VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);
    157:         }

    159:         if (0.0 == ldb->theta) {
    160:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);
    161:         } else if (1.0 == ldb->theta) {
    162:           VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
    163:         } else {
    164:           /*  Broyden update U=(1-theta)*P + theta*Q */
    165:           VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
    166:         }
    167:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
    168:         /*  Ensure positive definite */
    169:         VecAbs(ldb->invDnew);
    170:       }
    171:       if (ldb->sigma_hist > 0){
    172:         /*  Start with re-scaling on the newly computed diagonal */
    173:         if (0.5 == ldb->beta) {
    174:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
    175:             VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
    176:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

    178:             VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
    179:             VecDotBegin(ldb->W, lmvm->S[0], &stDs);
    180:             VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
    181:             VecDotEnd(ldb->W, lmvm->S[0], &stDs);

    183:             ss_sum = PetscRealPart(stDs);
    184:             yy_sum = PetscRealPart(ytDy);
    185:             ys_sum = ldb->yts[0];
    186:           } else {
    187:             VecCopy(ldb->invDnew, ldb->U);
    188:             VecReciprocal(ldb->U);

    190:             /*  Compute summations for scalar scaling */
    191:             yy_sum = 0;       /*  No safeguard required */
    192:             ys_sum = 0;       /*  No safeguard required */
    193:             ss_sum = 0;       /*  No safeguard required */
    194:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
    195:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    196:               VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
    197:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

    199:               VecDotBegin(ldb->W, lmvm->S[i], &stDs);
    200:               VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
    201:               VecDotEnd(ldb->W, lmvm->S[i], &stDs);
    202:               VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);

    204:               ss_sum += PetscRealPart(stDs);
    205:               ys_sum += ldb->yts[i];
    206:               yy_sum += PetscRealPart(ytDy);
    207:             }
    208:           }
    209:         } else if (0.0 == ldb->beta) {
    210:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
    211:             /*  Compute summations for scalar scaling */
    212:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

    214:             VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
    215:             VecDotBegin(ldb->W, ldb->W, &stDs);
    216:             VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
    217:             VecDotEnd(ldb->W, ldb->W, &stDs);

    219:             ys_sum = PetscRealPart(ytDs);
    220:             ss_sum = PetscRealPart(stDs);
    221:             yy_sum = ldb->yty[0];
    222:           } else {
    223:             VecCopy(ldb->invDnew, ldb->U);
    224:             VecReciprocal(ldb->U);

    226:             /*  Compute summations for scalar scaling */
    227:             yy_sum = 0;       /*  No safeguard required */
    228:             ys_sum = 0;       /*  No safeguard required */
    229:             ss_sum = 0;       /*  No safeguard required */
    230:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
    231:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    232:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

    234:               VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
    235:               VecDotBegin(ldb->W, ldb->W, &stDs);
    236:               VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
    237:               VecDotEnd(ldb->W, ldb->W, &stDs);

    239:               ss_sum += PetscRealPart(stDs);
    240:               ys_sum += PetscRealPart(ytDs);
    241:               yy_sum += ldb->yty[i];
    242:             }
    243:           }
    244:         } else if (1.0 == ldb->beta) {
    245:           /*  Compute summations for scalar scaling */
    246:           yy_sum = 0; /*  No safeguard required */
    247:           ys_sum = 0; /*  No safeguard required */
    248:           ss_sum = 0; /*  No safeguard required */
    249:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
    250:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    251:             VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);

    253:             VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
    254:             VecDotBegin(ldb->V, ldb->V, &ytDy);
    255:             VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
    256:             VecDotEnd(ldb->V, ldb->V, &ytDy);

    258:             yy_sum += PetscRealPart(ytDy);
    259:             ys_sum += PetscRealPart(ytDs);
    260:             ss_sum += ldb->sts[i];
    261:           }
    262:         } else {
    263:           VecCopy(ldb->invDnew, ldb->U);
    264:           VecPow(ldb->U, ldb->beta-1);

    266:           /*  Compute summations for scalar scaling */
    267:           yy_sum = 0; /*  No safeguard required */
    268:           ys_sum = 0; /*  No safeguard required */
    269:           ss_sum = 0; /*  No safeguard required */
    270:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
    271:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
    272:             VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
    273:             VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);

    275:             VecDotBegin(ldb->V, ldb->W, &ytDs);
    276:             VecDotBegin(ldb->V, ldb->V, &ytDy);
    277:             VecDotBegin(ldb->W, ldb->W, &stDs);
    278:             VecDotEnd(ldb->V, ldb->W, &ytDs);
    279:             VecDotEnd(ldb->V, ldb->V, &ytDy);
    280:             VecDotEnd(ldb->W, ldb->W, &stDs);

    282:             yy_sum += PetscRealPart(ytDy);
    283:             ys_sum += PetscRealPart(ytDs);
    284:             ss_sum += PetscRealPart(stDs);
    285:           }
    286:         }

    288:         if (0.0 == ldb->alpha) {
    289:           /*  Safeguard ys_sum  */
    290:           ys_sum = PetscMax(ldb->tol, ys_sum);

    292:           sigma = ss_sum / ys_sum;
    293:         } else if (1.0 == ldb->alpha) {
    294:           /* yy_sum is never 0; if it were, we'd be at the minimum */
    295:           sigma = ys_sum / yy_sum;
    296:         } else {
    297:           denom = 2.0*ldb->alpha*yy_sum;

    299:           /*  Safeguard denom */
    300:           denom = PetscMax(ldb->tol, denom);

    302:           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;
    303:         }
    304:       } else {
    305:         sigma = 1.0;
    306:       }
    307:       /*  If Q has small values, then Q^(r_beta - 1)
    308:        can have very large values.  Hence, ys_sum
    309:        and ss_sum can be infinity.  In this case,
    310:        sigma can either be not-a-number or infinity. */

    312:       if (PetscIsInfOrNanScalar(sigma)) {
    313:         /*  sigma is not-a-number; skip rescaling */
    314:       } else if (0.0 == sigma) {
    315:         /*  sigma is zero; this is a bad case; skip rescaling */
    316:       } else {
    317:         /*  sigma is positive */
    318:         VecScale(ldb->invDnew, sigma);
    319:       }

    321:       /* Combine the old diagonal and the new diagonal using a convex limiter */
    322:       if (1.0 == ldb->rho) {
    323:         VecCopy(ldb->invDnew, ldb->invD);
    324:       } else if (ldb->rho) {
    325:         VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);
    326:       }
    327:     } else {
    328:       MatLMVMReset(B, PETSC_FALSE);
    329:     }
    330:     /* End DiagBrdn update */

    332:   }
    333:   /* Save the solution and function to be used in the next update */
    334:   VecCopy(X, lmvm->Xprev);
    335:   VecCopy(F, lmvm->Fprev);
    336:   lmvm->prev_set = PETSC_TRUE;
    337:   return(0);
    338: }

    340: /*------------------------------------------------------------*/

    342: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
    343: {
    344:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
    345:   Mat_DiagBrdn      *bctx = (Mat_DiagBrdn*)bdata->ctx;
    346:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
    347:   Mat_DiagBrdn      *mctx = (Mat_DiagBrdn*)mdata->ctx;
    348:   PetscErrorCode    ierr;
    349:   PetscInt          i;

    352:   mctx->theta = bctx->theta;
    353:   mctx->alpha = bctx->alpha;
    354:   mctx->beta = bctx->beta;
    355:   mctx->rho = bctx->rho;
    356:   mctx->delta = bctx->delta;
    357:   mctx->delta_min = bctx->delta_min;
    358:   mctx->delta_max = bctx->delta_max;
    359:   mctx->tol = bctx->tol;
    360:   mctx->sigma = bctx->sigma;
    361:   mctx->sigma_hist = bctx->sigma_hist;
    362:   mctx->forward = bctx->forward;
    363:   VecCopy(bctx->invD, mctx->invD);
    364:   for (i=0; i<=bdata->k; ++i) {
    365:     mctx->yty[i] = bctx->yty[i];
    366:     mctx->yts[i] = bctx->yts[i];
    367:     mctx->sts[i] = bctx->sts[i];
    368:   }
    369:   return(0);
    370: }

    372: /*------------------------------------------------------------*/

    374: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
    375: {
    376:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    377:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    378:   PetscErrorCode    ierr;
    379:   PetscBool         isascii;

    382:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
    383:   if (isascii) {
    384:     PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);
    385:     PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
    386:     PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);
    387:   }
    388:   MatView_LMVM(B, pv);
    389:   return(0);
    390: }

    392: /*------------------------------------------------------------*/

    394: static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
    395: {
    396:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    397:   Mat_DiagBrdn       *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    398:   PetscErrorCode    ierr;

    401:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
    402:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
    403:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);
    404:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);
    405:   PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);
    406:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);
    407:   PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);
    408:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);
    409:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);
    410:   PetscOptionsTail();
    411:   if ((ldb->theta < 0.0) || (ldb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
    412:   if ((ldb->alpha < 0.0) || (ldb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
    413:   if ((ldb->rho < 0.0) || (ldb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
    414:   if (ldb->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
    415:   return(0);
    416: }

    418: /*------------------------------------------------------------*/

    420: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
    421: {
    422:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    423:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    424:   PetscErrorCode    ierr;

    427:   VecSet(ldb->invD, ldb->delta);
    428:   if (destructive && ldb->allocated) {
    429:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
    430:     VecDestroy(&ldb->invDnew);
    431:     VecDestroy(&ldb->invD);
    432:     VecDestroy(&ldb->BFGS);
    433:     VecDestroy(&ldb->DFP);
    434:     VecDestroy(&ldb->U);
    435:     VecDestroy(&ldb->V);
    436:     VecDestroy(&ldb->W);
    437:     ldb->allocated = PETSC_FALSE;
    438:   }
    439:   MatReset_LMVM(B, destructive);
    440:   return(0);
    441: }

    443: /*------------------------------------------------------------*/

    445: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
    446: {
    447:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    448:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    449:   PetscErrorCode    ierr;

    452:   MatAllocate_LMVM(B, X, F);
    453:   if (!ldb->allocated) {
    454:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
    455:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
    456:     VecDuplicate(lmvm->Xprev, &ldb->invD);
    457:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
    458:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
    459:     VecDuplicate(lmvm->Xprev, &ldb->U);
    460:     VecDuplicate(lmvm->Xprev, &ldb->V);
    461:     VecDuplicate(lmvm->Xprev, &ldb->W);
    462:     ldb->allocated = PETSC_TRUE;
    463:   }
    464:   return(0);
    465: }

    467: /*------------------------------------------------------------*/

    469: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
    470: {
    471:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    472:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    473:   PetscErrorCode    ierr;

    476:   if (ldb->allocated) {
    477:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
    478:     VecDestroy(&ldb->invDnew);
    479:     VecDestroy(&ldb->invD);
    480:     VecDestroy(&ldb->BFGS);
    481:     VecDestroy(&ldb->DFP);
    482:     VecDestroy(&ldb->U);
    483:     VecDestroy(&ldb->V);
    484:     VecDestroy(&ldb->W);
    485:     ldb->allocated = PETSC_FALSE;
    486:   }
    487:   PetscFree(lmvm->ctx);
    488:   MatDestroy_LMVM(B);
    489:   return(0);
    490: }

    492: /*------------------------------------------------------------*/

    494: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
    495: {
    496:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
    497:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
    498:   PetscErrorCode    ierr;

    501:   MatSetUp_LMVM(B);
    502:   if (!ldb->allocated) {
    503:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
    504:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
    505:     VecDuplicate(lmvm->Xprev, &ldb->invD);
    506:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
    507:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
    508:     VecDuplicate(lmvm->Xprev, &ldb->U);
    509:     VecDuplicate(lmvm->Xprev, &ldb->V);
    510:     VecDuplicate(lmvm->Xprev, &ldb->W);
    511:     ldb->allocated = PETSC_TRUE;
    512:   }
    513:   return(0);
    514: }

    516: /*------------------------------------------------------------*/

    518: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
    519: {
    520:   Mat_LMVM          *lmvm;
    521:   Mat_DiagBrdn      *ldb;
    522:   PetscErrorCode    ierr;

    525:   MatCreate_LMVM(B);
    526:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);
    527:   B->ops->setup = MatSetUp_DiagBrdn;
    528:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
    529:   B->ops->destroy = MatDestroy_DiagBrdn;
    530:   B->ops->solve = MatSolve_DiagBrdn;
    531:   B->ops->view = MatView_DiagBrdn;

    533:   lmvm = (Mat_LMVM*)B->data;
    534:   lmvm->square = PETSC_TRUE;
    535:   lmvm->m = 1;
    536:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
    537:   lmvm->ops->reset = MatReset_DiagBrdn;
    538:   lmvm->ops->mult = MatMult_DiagBrdn;
    539:   lmvm->ops->update = MatUpdate_DiagBrdn;
    540:   lmvm->ops->copy = MatCopy_DiagBrdn;

    542:   PetscNewLog(B, &ldb);
    543:   lmvm->ctx = (void*)ldb;
    544:   ldb->theta = 0.0;
    545:   ldb->alpha = 1.0;
    546:   ldb->rho = 1.0;
    547:   ldb->forward = PETSC_TRUE;
    548:   ldb->beta = 0.5;
    549:   ldb->sigma = 1.0;
    550:   ldb->delta = 1.0;
    551:   ldb->delta_min = 1e-7;
    552:   ldb->delta_max = 100.0;
    553:   ldb->tol = 1e-8;
    554:   ldb->sigma_hist = 1;
    555:   ldb->allocated = PETSC_FALSE;
    556:   return(0);
    557: }

    559: /*------------------------------------------------------------*/

    561: /*@
    562:    MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
    563:    for approximating Hessians. It consists of a convex combination of DFP and BFGS
    564:    diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
    565:    To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
    566:    We also ensure positive definiteness by taking the VecAbs() of the final vector.

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

    574:    In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults
    575:    and MatSolves, the matrix must first be created using MatCreate() and MatSetType(),
    576:    followed by MatLMVMAllocate(). Then it will be available for updating
    577:    (via MatLMVMUpdate) in one's favored solver implementation.
    578:    This allows for MPI compatibility.

    580:    Collective

    582:    Input Parameters:
    583: +  comm - MPI communicator, set to PETSC_COMM_SELF
    584: .  n - number of local rows for storage vectors
    585: -  N - global size of the storage vectors

    587:    Output Parameter:
    588: .  B - the matrix

    590:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
    591:    paradigm instead of this routine directly.

    593:    Options Database Keys:
    594: +   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    595: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
    596: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    597: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
    598: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
    599: .   -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
    600: -   -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal

    602:    Level: intermediate

    604: .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
    605:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
    606: @*/
    607: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    608: {
    609:   PetscErrorCode    ierr;

    612:   MatCreate(comm, B);
    613:   MatSetSizes(*B, n, n, N, N);
    614:   MatSetType(*B, MATLMVMDIAGBROYDEN);
    615:   MatSetUp(*B);
    616:   return(0);
    617: }