Actual source code: symbrdn.c

      1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
      2: #include <../src/ksp/ksp/utils/lmvm/dense/denseqn.h>
      3: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
      4: #include <petsc/private/kspimpl.h>
      5: #include <petscdevice.h>

      7: const char *const MatLMVMSymBroydenScaleTypes[] = {"NONE", "SCALAR", "DIAGONAL", "USER", "MatLMVMSymBrdnScaleType", "MAT_LMVM_SYMBROYDEN_SCALING_", NULL};

      9: /*
     10:   The solution method below is the matrix-free implementation of
     11:   Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation
     12:   and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).

     14:   Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
     15:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
     16:   repeated calls of MatSolve without incurring redundant computation.

     18:   dX <- J0^{-1} * F

     20:   for i=0,1,2,...,k
     21:     # Q[i] = (B_i)^T{-1} Y[i]

     23:     rho = 1.0 / (Y[i]^T S[i])
     24:     alpha = rho * (S[i]^T F)
     25:     zeta = 1.0 / (Y[i]^T Q[i])
     26:     gamma = zeta * (Y[i]^T dX)

     28:     dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
     29:     W <- (rho * S[i]) - (zeta * Q[i])
     30:     dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
     31:   end
     32: */
     33: static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
     34: {
     35:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
     36:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
     37:   PetscInt     i, j;
     38:   PetscReal    numer;
     39:   PetscScalar  sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;

     41:   PetscFunctionBegin;
     42:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
     43:   if (lsb->phi == 0.0) {
     44:     PetscCall(MatSolve_LMVMBFGS(B, F, dX));
     45:     PetscFunctionReturn(PETSC_SUCCESS);
     46:   }
     47:   if (lsb->phi == 1.0) {
     48:     PetscCall(MatSolve_LMVMDFP(B, F, dX));
     49:     PetscFunctionReturn(PETSC_SUCCESS);
     50:   }

     52:   VecCheckSameSize(F, 2, dX, 3);
     53:   VecCheckMatCompatible(B, dX, 3, F, 2);

     55:   if (lsb->needP) {
     56:     /* Start the loop for (P[k] = (B_k) * S[k]) */
     57:     for (i = 0; i <= lmvm->k; ++i) {
     58:       PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
     59:       /* Compute the necessary dot products */
     60:       PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
     61:       for (j = 0; j < i; ++j) {
     62:         yjtsi = lsb->workscalar[j];
     63:         PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
     64:         /* Compute the pure BFGS component of the forward product */
     65:         PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
     66:         /* Tack on the convexly scaled extras to the forward product */
     67:         if (lsb->phi > 0.0) {
     68:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
     69:           PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
     70:           PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
     71:         }
     72:       }
     73:       PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
     74:       lsb->stp[i] = PetscRealPart(stp);
     75:     }
     76:     lsb->needP = PETSC_FALSE;
     77:   }
     78:   if (lsb->needQ) {
     79:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
     80:     for (i = 0; i <= lmvm->k; ++i) {
     81:       PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]));
     82:       /* Compute the necessary dot products */
     83:       PetscCall(VecMDot(lmvm->Y[i], i, lmvm->S, lsb->workscalar));
     84:       for (j = 0; j < i; ++j) {
     85:         sjtyi = lsb->workscalar[j];
     86:         PetscCall(VecDot(lmvm->Y[j], lsb->Q[i], &yjtqi));
     87:         /* Compute the pure DFP component of the inverse application*/
     88:         PetscCall(VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi) / lsb->ytq[j], PetscRealPart(sjtyi) / lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]));
     89:         /* Tack on the convexly scaled extras to the inverse application*/
     90:         if (lsb->psi[j] > 0.0) {
     91:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]));
     92:           PetscCall(VecDot(lsb->work, lmvm->Y[i], &wtyi));
     93:           PetscCall(VecAXPY(lsb->Q[i], lsb->psi[j] * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work));
     94:         }
     95:       }
     96:       PetscCall(VecDot(lmvm->Y[i], lsb->Q[i], &ytq));
     97:       lsb->ytq[i] = PetscRealPart(ytq);
     98:       if (lsb->phi == 1.0) {
     99:         lsb->psi[i] = 0.0;
    100:       } else if (lsb->phi == 0.0) {
    101:         lsb->psi[i] = 1.0;
    102:       } else {
    103:         numer       = (1.0 - lsb->phi) * lsb->yts[i] * lsb->yts[i];
    104:         lsb->psi[i] = numer / (numer + (lsb->phi * lsb->ytq[i] * lsb->stp[i]));
    105:       }
    106:     }
    107:     lsb->needQ = PETSC_FALSE;
    108:   }

    110:   /* Start the outer iterations for ((B^{-1}) * dX) */
    111:   PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
    112:   /* Get all the dot products we need */
    113:   PetscCall(VecMDot(F, lmvm->k + 1, lmvm->S, lsb->workscalar));
    114:   for (i = 0; i <= lmvm->k; ++i) {
    115:     stf = lsb->workscalar[i];
    116:     PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
    117:     /* Compute the pure DFP component */
    118:     PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]));
    119:     /* Tack on the convexly scaled extras */
    120:     PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]));
    121:     PetscCall(VecDot(lsb->work, F, &wtf));
    122:     PetscCall(VecAXPY(dX, lsb->psi[i] * lsb->ytq[i] * PetscRealPart(wtf), lsb->work));
    123:   }
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: /*
    128:   The forward-product below is the matrix-free implementation of
    129:   Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
    130:   Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).

    132:   P[i] = (B_i)*S[i] terms are computed ahead of time whenever
    133:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
    134:   repeated calls of MatMult inside KSP solvers without unnecessarily
    135:   recomputing P[i] terms in expensive nested-loops.

    137:   Z <- J0 * X

    139:   for i=0,1,2,...,k
    140:     # P[i] = (B_k) * S[i]

    142:     rho = 1.0 / (Y[i]^T S[i])
    143:     alpha = rho * (Y[i]^T F)
    144:     zeta = 1.0 / (S[i]^T P[i])
    145:     gamma = zeta * (S[i]^T dX)

    147:     dX <- dX - (gamma * P[i]) + (alpha * S[i])
    148:     W <- (rho * Y[i]) - (zeta * P[i])
    149:     dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
    150:   end
    151: */
    152: static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
    153: {
    154:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    155:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    156:   PetscInt     i, j;
    157:   PetscScalar  sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;

    159:   PetscFunctionBegin;
    160:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
    161:   if (lsb->phi == 0.0) {
    162:     PetscCall(MatMult_LMVMBFGS(B, X, Z));
    163:     PetscFunctionReturn(PETSC_SUCCESS);
    164:   }
    165:   if (lsb->phi == 1.0) {
    166:     PetscCall(MatMult_LMVMDFP(B, X, Z));
    167:     PetscFunctionReturn(PETSC_SUCCESS);
    168:   }

    170:   VecCheckSameSize(X, 2, Z, 3);
    171:   VecCheckMatCompatible(B, X, 2, Z, 3);

    173:   if (lsb->needP) {
    174:     /* Start the loop for (P[k] = (B_k) * S[k]) */
    175:     for (i = 0; i <= lmvm->k; ++i) {
    176:       PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
    177:       /* Compute the necessary dot products */
    178:       PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
    179:       for (j = 0; j < i; ++j) {
    180:         yjtsi = lsb->workscalar[j];
    181:         PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
    182:         /* Compute the pure BFGS component of the forward product */
    183:         PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
    184:         /* Tack on the convexly scaled extras to the forward product */
    185:         if (lsb->phi > 0.0) {
    186:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
    187:           PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
    188:           PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
    189:         }
    190:       }
    191:       PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
    192:       lsb->stp[i] = PetscRealPart(stp);
    193:     }
    194:     lsb->needP = PETSC_FALSE;
    195:   }

    197:   /* Start the outer iterations for (B * X) */
    198:   PetscCall(MatSymBrdnApplyJ0Fwd(B, X, Z));
    199:   /* Get all the dot products we need */
    200:   PetscCall(VecMDot(X, lmvm->k + 1, lmvm->Y, lsb->workscalar));
    201:   for (i = 0; i <= lmvm->k; ++i) {
    202:     ytx = lsb->workscalar[i];
    203:     PetscCall(VecDot(lmvm->S[i], Z, &stz));
    204:     /* Compute the pure BFGS component */
    205:     PetscCall(VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]));
    206:     /* Tack on the convexly scaled extras */
    207:     PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]));
    208:     PetscCall(VecDot(lsb->work, X, &wtx));
    209:     PetscCall(VecAXPY(Z, lsb->phi * lsb->stp[i] * PetscRealPart(wtx), lsb->work));
    210:   }
    211:   PetscFunctionReturn(PETSC_SUCCESS);
    212: }

    214: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
    215: {
    216:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    217:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    218:   Mat_LMVM     *dbase;
    219:   Mat_DiagBrdn *dctx;
    220:   PetscInt      old_k, i;
    221:   PetscReal     curvtol, ytytmp;
    222:   PetscScalar   curvature, ststmp;

    224:   PetscFunctionBegin;
    225:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
    226:   if (lmvm->prev_set) {
    227:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
    228:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
    229:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

    231:     /* Test if the updates can be accepted */
    232:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ytytmp));
    233:     if (ytytmp < lmvm->eps) curvtol = 0.0;
    234:     else curvtol = lmvm->eps * ytytmp;

    236:     if (PetscRealPart(curvature) > curvtol) {
    237:       /* Update is good, accept it */
    238:       lsb->watchdog = 0;
    239:       lsb->needP = lsb->needQ = PETSC_TRUE;
    240:       old_k                   = lmvm->k;
    241:       PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
    242:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
    243:       if (old_k == lmvm->k) {
    244:         for (i = 0; i <= lmvm->k - 1; ++i) {
    245:           lsb->yts[i] = lsb->yts[i + 1];
    246:           lsb->yty[i] = lsb->yty[i + 1];
    247:           lsb->sts[i] = lsb->sts[i + 1];
    248:         }
    249:       }
    250:       /* Update history of useful scalars */
    251:       PetscCall(VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &ststmp));
    252:       lsb->yts[lmvm->k] = PetscRealPart(curvature);
    253:       lsb->yty[lmvm->k] = ytytmp;
    254:       lsb->sts[lmvm->k] = PetscRealPart(ststmp);
    255:       /* Compute the scalar scale if necessary */
    256:       if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
    257:     } else {
    258:       /* Update is bad, skip it */
    259:       ++lmvm->nrejects;
    260:       ++lsb->watchdog;
    261:     }
    262:   } else {
    263:     switch (lsb->scale_type) {
    264:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    265:       dbase = (Mat_LMVM *)lsb->D->data;
    266:       dctx  = (Mat_DiagBrdn *)dbase->ctx;
    267:       PetscCall(VecSet(dctx->invD, lsb->delta));
    268:       break;
    269:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    270:       lsb->sigma = lsb->delta;
    271:       break;
    272:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    273:       lsb->sigma = 1.0;
    274:       break;
    275:     default:
    276:       break;
    277:     }
    278:   }

    280:   /* Update the scaling */
    281:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(lsb->D, X, F));

    283:   if (lsb->watchdog > lsb->max_seq_rejects) {
    284:     PetscCall(MatLMVMReset(B, PETSC_FALSE));
    285:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
    286:   }

    288:   /* Save the solution and function to be used in the next update */
    289:   PetscCall(VecCopy(X, lmvm->Xprev));
    290:   PetscCall(VecCopy(F, lmvm->Fprev));
    291:   lmvm->prev_set = PETSC_TRUE;
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
    296: {
    297:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
    298:   Mat_SymBrdn *blsb  = (Mat_SymBrdn *)bdata->ctx;
    299:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
    300:   Mat_SymBrdn *mlsb  = (Mat_SymBrdn *)mdata->ctx;
    301:   PetscInt     i;

    303:   PetscFunctionBegin;
    304:   mlsb->phi   = blsb->phi;
    305:   mlsb->needP = blsb->needP;
    306:   mlsb->needQ = blsb->needQ;
    307:   for (i = 0; i <= bdata->k; ++i) {
    308:     mlsb->stp[i] = blsb->stp[i];
    309:     mlsb->ytq[i] = blsb->ytq[i];
    310:     mlsb->yts[i] = blsb->yts[i];
    311:     mlsb->psi[i] = blsb->psi[i];
    312:     PetscCall(VecCopy(blsb->P[i], mlsb->P[i]));
    313:     PetscCall(VecCopy(blsb->Q[i], mlsb->Q[i]));
    314:   }
    315:   mlsb->scale_type      = blsb->scale_type;
    316:   mlsb->alpha           = blsb->alpha;
    317:   mlsb->beta            = blsb->beta;
    318:   mlsb->rho             = blsb->rho;
    319:   mlsb->delta           = blsb->delta;
    320:   mlsb->sigma_hist      = blsb->sigma_hist;
    321:   mlsb->watchdog        = blsb->watchdog;
    322:   mlsb->max_seq_rejects = blsb->max_seq_rejects;
    323:   switch (blsb->scale_type) {
    324:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    325:     mlsb->sigma = blsb->sigma;
    326:     break;
    327:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    328:     PetscCall(MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN));
    329:     break;
    330:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    331:     mlsb->sigma = 1.0;
    332:     break;
    333:   default:
    334:     break;
    335:   }
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
    340: {
    341:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
    342:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    343:   Mat_LMVM     *dbase;
    344:   Mat_DiagBrdn *dctx;

    346:   PetscFunctionBegin;
    347:   lsb->watchdog = 0;
    348:   lsb->needP = lsb->needQ = PETSC_TRUE;
    349:   if (lsb->allocated) {
    350:     if (destructive) {
    351:       PetscCall(VecDestroy(&lsb->work));
    352:       PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
    353:       PetscCall(PetscFree(lsb->psi));
    354:       PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
    355:       PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
    356:       switch (lsb->scale_type) {
    357:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    358:         PetscCall(MatLMVMReset(lsb->D, PETSC_TRUE));
    359:         break;
    360:       default:
    361:         break;
    362:       }
    363:       lsb->allocated = PETSC_FALSE;
    364:     } else {
    365:       PetscCall(PetscMemzero(lsb->psi, lmvm->m));
    366:       switch (lsb->scale_type) {
    367:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    368:         lsb->sigma = lsb->delta;
    369:         break;
    370:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    371:         PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
    372:         dbase = (Mat_LMVM *)lsb->D->data;
    373:         dctx  = (Mat_DiagBrdn *)dbase->ctx;
    374:         PetscCall(VecSet(dctx->invD, lsb->delta));
    375:         break;
    376:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    377:         lsb->sigma = 1.0;
    378:         break;
    379:       default:
    380:         break;
    381:       }
    382:     }
    383:   }
    384:   PetscCall(MatReset_LMVM(B, destructive));
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }

    388: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
    389: {
    390:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    391:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    393:   PetscFunctionBegin;
    394:   PetscCall(MatAllocate_LMVM(B, X, F));
    395:   if (!lsb->allocated) {
    396:     PetscCall(VecDuplicate(X, &lsb->work));
    397:     if (lmvm->m > 0) {
    398:       PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
    399:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
    400:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->P));
    401:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->Q));
    402:     }
    403:     switch (lsb->scale_type) {
    404:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    405:       PetscCall(MatLMVMAllocate(lsb->D, X, F));
    406:       break;
    407:     default:
    408:       break;
    409:     }
    410:     lsb->allocated = PETSC_TRUE;
    411:   }
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

    415: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
    416: {
    417:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    418:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    420:   PetscFunctionBegin;
    421:   if (lsb->allocated) {
    422:     PetscCall(VecDestroy(&lsb->work));
    423:     PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
    424:     PetscCall(PetscFree(lsb->psi));
    425:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
    426:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
    427:     lsb->allocated = PETSC_FALSE;
    428:   }
    429:   PetscCall(MatDestroy(&lsb->D));
    430:   PetscCall(PetscFree(lmvm->ctx));
    431:   PetscCall(MatDestroy_LMVM(B));
    432:   PetscFunctionReturn(PETSC_SUCCESS);
    433: }

    435: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
    436: {
    437:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    438:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    439:   PetscInt     n, N;

    441:   PetscFunctionBegin;
    442:   PetscCall(MatSetUp_LMVM(B));
    443:   if (!lsb->allocated) {
    444:     PetscCall(VecDuplicate(lmvm->Xprev, &lsb->work));
    445:     if (lmvm->m > 0) {
    446:       PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
    447:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
    448:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P));
    449:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q));
    450:     }
    451:     switch (lsb->scale_type) {
    452:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    453:       PetscCall(MatGetLocalSize(B, &n, &n));
    454:       PetscCall(MatGetSize(B, &N, &N));
    455:       PetscCall(MatSetSizes(lsb->D, n, n, N, N));
    456:       PetscCall(MatSetUp(lsb->D));
    457:       break;
    458:     default:
    459:       break;
    460:     }
    461:     lsb->allocated = PETSC_TRUE;
    462:   }
    463:   PetscFunctionReturn(PETSC_SUCCESS);
    464: }

    466: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
    467: {
    468:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    469:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    470:   PetscBool    isascii;

    472:   PetscFunctionBegin;
    473:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
    474:   if (isascii) {
    475:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale type: %s\n", MatLMVMSymBroydenScaleTypes[lsb->scale_type]));
    476:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", lsb->sigma_hist));
    477:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)lsb->alpha, (double)lsb->beta, (double)lsb->rho));
    478:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factors: phi=%g, theta=%g\n", (double)lsb->phi, (double)lsb->theta));
    479:   }
    480:   PetscCall(MatView_LMVM(B, pv));
    481:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatView(lsb->D, pv));
    482:   PetscFunctionReturn(PETSC_SUCCESS);
    483: }

    485: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
    486: {
    487:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    488:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    490:   PetscFunctionBegin;
    491:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
    492:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
    493:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL, 0.0, 1.0));
    494:   PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
    495:   PetscOptionsHeadEnd();
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }

    499: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(Mat B, PetscOptionItems *PetscOptionsObject)
    500: {
    501:   Mat_LMVM                  *lmvm = (Mat_LMVM *)B->data;
    502:   Mat_SymBrdn               *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    503:   Mat_LMVM                  *dbase;
    504:   Mat_DiagBrdn              *dctx;
    505:   MatLMVMSymBroydenScaleType stype = lsb->scale_type;
    506:   PetscBool                  flg;

    508:   PetscFunctionBegin;
    509:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL));
    510:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL, 0.0, 1.0));
    511:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL, 0.0, 1.0));
    512:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL, 0.0, 1.0));
    513:   PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", lsb->sigma_hist, &lsb->sigma_hist, NULL, 1));
    514:   PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
    515:   if (flg) PetscCall(MatLMVMSymBroydenSetScaleType(B, stype));
    516:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
    517:     const char *prefix;

    519:     PetscCall(MatGetOptionsPrefix(B, &prefix));
    520:     PetscCall(MatSetOptionsPrefix(lsb->D, prefix));
    521:     PetscCall(MatAppendOptionsPrefix(lsb->D, "J0_"));
    522:     PetscCall(MatSetFromOptions(lsb->D));
    523:     dbase            = (Mat_LMVM *)lsb->D->data;
    524:     dctx             = (Mat_DiagBrdn *)dbase->ctx;
    525:     dctx->delta_min  = lsb->delta_min;
    526:     dctx->delta_max  = lsb->delta_max;
    527:     dctx->theta      = lsb->theta;
    528:     dctx->rho        = lsb->rho;
    529:     dctx->alpha      = lsb->alpha;
    530:     dctx->beta       = lsb->beta;
    531:     dctx->sigma_hist = lsb->sigma_hist;
    532:   }
    533:   PetscFunctionReturn(PETSC_SUCCESS);
    534: }

    536: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
    537: {
    538:   Mat_LMVM    *lmvm;
    539:   Mat_SymBrdn *lsb;

    541:   PetscFunctionBegin;
    542:   PetscCall(MatCreate_LMVM(B));
    543:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN));
    544:   PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
    545:   PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
    546:   B->ops->view           = MatView_LMVMSymBrdn;
    547:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
    548:   B->ops->setup          = MatSetUp_LMVMSymBrdn;
    549:   B->ops->destroy        = MatDestroy_LMVMSymBrdn;

    551:   lmvm                = (Mat_LMVM *)B->data;
    552:   lmvm->square        = PETSC_TRUE;
    553:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
    554:   lmvm->ops->reset    = MatReset_LMVMSymBrdn;
    555:   lmvm->ops->update   = MatUpdate_LMVMSymBrdn;
    556:   lmvm->ops->mult     = MatMult_LMVMSymBrdn;
    557:   lmvm->ops->solve    = MatSolve_LMVMSymBrdn;
    558:   lmvm->ops->copy     = MatCopy_LMVMSymBrdn;

    560:   PetscCall(PetscNew(&lsb));
    561:   lmvm->ctx      = (void *)lsb;
    562:   lsb->allocated = PETSC_FALSE;
    563:   lsb->needP = lsb->needQ = PETSC_TRUE;
    564:   lsb->phi                = 0.125;
    565:   lsb->theta              = 0.125;
    566:   lsb->alpha              = 1.0;
    567:   lsb->rho                = 1.0;
    568:   lsb->beta               = 0.5;
    569:   lsb->sigma              = 1.0;
    570:   lsb->delta              = 1.0;
    571:   lsb->delta_min          = 1e-7;
    572:   lsb->delta_max          = 100.0;
    573:   lsb->sigma_hist         = 1;
    574:   lsb->scale_type         = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
    575:   lsb->watchdog           = 0;
    576:   lsb->max_seq_rejects    = lmvm->m / 2;

    578:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &lsb->D));
    579:   PetscCall(MatSetType(lsb->D, MATLMVMDIAGBROYDEN));
    580:   PetscFunctionReturn(PETSC_SUCCESS);
    581: }

    583: /*@
    584:   MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
    585:   in the SymBrdn approximations (also works for BFGS and DFP).

    587:   Input Parameters:
    588: + B     - `MATLMVM` matrix
    589: - delta - initial value for diagonal scaling

    591:   Level: intermediate

    593: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`
    594: @*/
    595: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
    596: {
    597:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
    598:   PetscBool is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn, is_dbfgs, is_ddfp, is_dqn;
    599:   PetscReal del_min, del_max, del_buf;

    601:   PetscFunctionBegin;
    602:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs));
    603:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
    604:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
    605:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
    606:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp));
    607:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn));
    608:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn));

    610:   if (is_bfgs || is_dfp || is_symbrdn || is_symbadbrdn) {
    611:     Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;

    613:     lsb     = (Mat_SymBrdn *)lmvm->ctx;
    614:     del_min = lsb->delta_min;
    615:     del_max = lsb->delta_max;
    616:   } else if (is_dbfgs || is_ddfp || is_dqn) {
    617:     Mat_DQN      *lqn     = (Mat_DQN *)lmvm->ctx;
    618:     Mat_LMVM     *dbase   = (Mat_LMVM *)lqn->diag_qn->data;
    619:     Mat_DiagBrdn *diagctx = (Mat_DiagBrdn *)dbase->ctx;

    621:     del_min = diagctx->delta_min;
    622:     del_max = diagctx->delta_max;
    623:   } else {
    624:     SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling only available for SymBrdn-derived types (DBFGS, BFGS, DDFP, DFP, SymBrdn, SymBadBrdn");
    625:   }

    627:   del_buf = PetscAbsReal(PetscRealPart(delta));
    628:   del_buf = PetscMin(del_buf, del_max);
    629:   del_buf = PetscMax(del_buf, del_min);
    630:   if (is_dbfgs || is_ddfp || is_dqn) {
    631:     Mat_DQN      *lqn     = (Mat_DQN *)lmvm->ctx;
    632:     Mat_LMVM     *dbase   = (Mat_LMVM *)lqn->diag_qn->data;
    633:     Mat_DiagBrdn *diagctx = (Mat_DiagBrdn *)dbase->ctx;

    635:     diagctx->delta = del_buf;
    636:   } else {
    637:     Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;

    639:     lsb        = (Mat_SymBrdn *)lmvm->ctx;
    640:     lsb->delta = del_buf;
    641:   }
    642:   PetscFunctionReturn(PETSC_SUCCESS);
    643: }

    645: /*@
    646:   MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.

    648:   Input Parameters:
    649: + B     - the `MATLMVM` matrix
    650: - stype - scale type, see `MatLMVMSymBroydenScaleType`

    652:   Options Database Key:
    653: . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type

    655:   Level: intermediate

    657:   MatLMVMSymBrdnScaleTypes\:
    658: +   `MAT_LMVM_SYMBROYDEN_SCALE_NONE` - initial Hessian is the identity matrix
    659: .   `MAT_LMVM_SYMBROYDEN_SCALE_SCALAR` - use the Shanno scalar as the initial Hessian
    660: -   `MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL` - use a diagonalized BFGS update as the initial Hessian

    662: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`, `MatLMVMSymBroydenScaleType`
    663: @*/
    664: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
    665: {
    666:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    667:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    669:   PetscFunctionBegin;
    671:   lsb->scale_type = stype;
    672:   PetscFunctionReturn(PETSC_SUCCESS);
    673: }

    675: /*@
    676:   MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
    677:   for approximating Jacobians.

    679:   Collective

    681:   Input Parameters:
    682: + comm - MPI communicator, set to `PETSC_COMM_SELF`
    683: . n    - number of local rows for storage vectors
    684: - N    - global size of the storage vectors

    686:   Output Parameter:
    687: . B - the matrix

    689:   Options Database Keys:
    690: + -mat_lmvm_phi        - (developer) convex ratio between BFGS and DFP components of the update
    691: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
    692: . -mat_lmvm_theta      - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    693: . -mat_lmvm_rho        - (developer) update limiter for the J0 scaling
    694: . -mat_lmvm_alpha      - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    695: . -mat_lmvm_beta       - (developer) exponential factor for the diagonal J0 scaling
    696: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

    698:   Level: intermediate

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

    704:   L-SymBrdn is a convex combination of L-DFP and
    705:   L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
    706:   phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
    707:   to be symmetric positive-definite.

    709:   To use the L-SymBrdn matrix with other vector types, the matrix must be
    710:   created using MatCreate() and MatSetType(), followed by `MatLMVMAllocate()`.
    711:   This ensures that the internal storage and work vectors are duplicated from the
    712:   correct type of vector.

    714: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
    715:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
    716: @*/
    717: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    718: {
    719:   PetscFunctionBegin;
    720:   PetscCall(KSPInitializePackage());
    721:   PetscCall(MatCreate(comm, B));
    722:   PetscCall(MatSetSizes(*B, n, n, N, N));
    723:   PetscCall(MatSetType(*B, MATLMVMSYMBROYDEN));
    724:   PetscCall(MatSetUp(*B));
    725:   PetscFunctionReturn(PETSC_SUCCESS);
    726: }

    728: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
    729: {
    730:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    731:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    733:   PetscFunctionBegin;
    734:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
    735:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
    736:     PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
    737:   } else {
    738:     switch (lsb->scale_type) {
    739:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    740:       PetscCall(VecAXPBY(Z, 1.0 / lsb->sigma, 0.0, X));
    741:       break;
    742:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    743:       PetscCall(MatMult(lsb->D, X, Z));
    744:       break;
    745:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    746:     default:
    747:       PetscCall(VecCopy(X, Z));
    748:       break;
    749:     }
    750:   }
    751:   PetscFunctionReturn(PETSC_SUCCESS);
    752: }

    754: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
    755: {
    756:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    757:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

    759:   PetscFunctionBegin;
    760:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
    761:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
    762:     PetscCall(MatLMVMApplyJ0Inv(B, F, dX));
    763:   } else {
    764:     switch (lsb->scale_type) {
    765:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    766:       PetscCall(VecAXPBY(dX, lsb->sigma, 0.0, F));
    767:       break;
    768:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
    769:       PetscCall(MatSolve(lsb->D, F, dX));
    770:       break;
    771:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
    772:     default:
    773:       PetscCall(VecCopy(F, dX));
    774:       break;
    775:     }
    776:   }
    777:   PetscFunctionReturn(PETSC_SUCCESS);
    778: }

    780: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
    781: {
    782:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    783:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
    784:   PetscInt     i, start;
    785:   PetscReal    a, b, c, sig1, sig2, signew;

    787:   PetscFunctionBegin;
    788:   if (lsb->sigma_hist == 0) {
    789:     signew = 1.0;
    790:   } else {
    791:     start  = PetscMax(0, lmvm->k - lsb->sigma_hist + 1);
    792:     signew = 0.0;
    793:     if (lsb->alpha == 1.0) {
    794:       for (i = start; i <= lmvm->k; ++i) signew += lsb->yts[i] / lsb->yty[i];
    795:     } else if (lsb->alpha == 0.5) {
    796:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yty[i];
    797:       signew = PetscSqrtReal(signew);
    798:     } else if (lsb->alpha == 0.0) {
    799:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yts[i];
    800:     } else {
    801:       /* compute coefficients of the quadratic */
    802:       a = b = c = 0.0;
    803:       for (i = start; i <= lmvm->k; ++i) {
    804:         a += lsb->yty[i];
    805:         b += lsb->yts[i];
    806:         c += lsb->sts[i];
    807:       }
    808:       a *= lsb->alpha;
    809:       b *= -(2.0 * lsb->alpha - 1.0);
    810:       c *= lsb->alpha - 1.0;
    811:       /* use quadratic formula to find roots */
    812:       sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
    813:       sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
    814:       /* accept the positive root as the scalar */
    815:       if (sig1 > 0.0) {
    816:         signew = sig1;
    817:       } else if (sig2 > 0.0) {
    818:         signew = sig2;
    819:       } else {
    820:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
    821:       }
    822:     }
    823:   }
    824:   lsb->sigma = lsb->rho * signew + (1.0 - lsb->rho) * lsb->sigma;
    825:   PetscFunctionReturn(PETSC_SUCCESS);
    826: }