Actual source code: symbrdnrescale.c

      1: #include <petscdevice.h>
      2: #include "symbrdnrescale.h"

      4: PetscLogEvent SBRDN_Rescale;

      6: const char *const MatLMVMSymBroydenScaleTypes[] = {"none", "scalar", "diagonal", "user", "decide", "MatLMVMSymBroydenScaleType", "MAT_LMVM_SYMBROYDEN_SCALING_", NULL};

      8: static PetscErrorCode SymBroydenRescaleUpdateScalar(Mat B, SymBroydenRescale ldb)
      9: {
     10:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
     11:   PetscReal a, b, c, signew;
     12:   PetscReal sigma_inv, sigma;
     13:   PetscInt  oldest, next;

     15:   PetscFunctionBegin;
     16:   next   = ldb->k;
     17:   oldest = PetscMax(0, ldb->k - ldb->sigma_hist);
     18:   PetscCall(MatNorm(lmvm->J0, NORM_INFINITY, &sigma_inv));
     19:   sigma = 1.0 / sigma_inv;
     20:   if (ldb->sigma_hist == 0) {
     21:     signew = 1.0;
     22:   } else {
     23:     signew = 0.0;
     24:     if (ldb->alpha == 1.0) {
     25:       for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->yts[i] / ldb->yty[i];
     26:     } else if (ldb->alpha == 0.5) {
     27:       for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->sts[i] / ldb->yty[i];
     28:       signew = PetscSqrtReal(signew);
     29:     } else if (ldb->alpha == 0.0) {
     30:       for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->sts[i] / ldb->yts[i];
     31:     } else {
     32:       /* compute coefficients of the quadratic */
     33:       a = b = c = 0.0;
     34:       for (PetscInt i = 0; i < next - oldest; ++i) {
     35:         a += ldb->yty[i];
     36:         b += ldb->yts[i];
     37:         c += ldb->sts[i];
     38:       }
     39:       a *= ldb->alpha;
     40:       b *= -(2.0 * ldb->alpha - 1.0);
     41:       c *= ldb->alpha - 1.0;
     42:       /* use quadratic formula to find roots */
     43:       PetscReal sqrtdisc = PetscSqrtReal(b * b - 4 * a * c);
     44:       if (b >= 0.0) {
     45:         if (a >= 0.0) {
     46:           signew = (2 * c) / (-b - sqrtdisc);
     47:         } else {
     48:           signew = (-b - sqrtdisc) / (2 * a);
     49:         }
     50:       } else {
     51:         if (a >= 0.0) {
     52:           signew = (-b + sqrtdisc) / (2 * a);
     53:         } else {
     54:           signew = (2 * c) / (-b + sqrtdisc);
     55:         }
     56:       }
     57:       PetscCheck(signew > 0.0, PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
     58:     }
     59:   }
     60:   sigma = ldb->rho * signew + (1.0 - ldb->rho) * sigma;
     61:   PetscCall(MatLMVMSetJ0Scale(B, 1.0 / sigma));
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: static PetscErrorCode DiagonalUpdate(SymBroydenRescale ldb, Vec D, Vec s, Vec y, Vec V, Vec W, Vec BFGS, Vec DFP, PetscReal theta, PetscReal yts)
     66: {
     67:   PetscFunctionBegin;
     68:   /*  V = |y o y| */
     69:   PetscCall(VecPointwiseMult(V, y, y));
     70:   if (PetscDefined(USE_COMPLEX)) PetscCall(VecAbs(V));

     72:   /*  W = D o s */
     73:   PetscReal stDs;
     74:   PetscCall(VecPointwiseMult(W, D, s));
     75:   PetscCall(VecDotRealPart(W, s, &stDs));

     77:   PetscCall(VecAXPY(D, 1.0 / yts, ldb->V));

     79:   /*  Safeguard stDs */
     80:   stDs = PetscMax(stDs, ldb->tol);

     82:   if (theta != 1.0) {
     83:     /*  BFGS portion of the update */

     85:     /*  U = |(D o s) o (D o s)| */
     86:     PetscCall(VecPointwiseMult(BFGS, W, W));
     87:     if (PetscDefined(USE_COMPLEX)) PetscCall(VecAbs(BFGS));

     89:     /*  Assemble */
     90:     PetscCall(VecScale(BFGS, -1.0 / stDs));
     91:   }

     93:   if (theta != 0.0) {
     94:     /*  DFP portion of the update */
     95:     /*  U = Real(conj(y) o D o s) */
     96:     PetscCall(VecCopy(y, DFP));
     97:     PetscCall(VecConjugate(DFP));
     98:     PetscCall(VecPointwiseMult(DFP, DFP, W));
     99:     if (PetscDefined(USE_COMPLEX)) {
    100:       PetscCall(VecCopy(DFP, W));
    101:       PetscCall(VecConjugate(W));
    102:       PetscCall(VecAXPY(DFP, 1.0, W));
    103:     } else {
    104:       PetscCall(VecScale(DFP, 2.0));
    105:     }

    107:     /*  Assemble */
    108:     PetscCall(VecAXPBY(DFP, stDs / yts, -1.0, V));
    109:   }

    111:   if (theta == 0.0) {
    112:     PetscCall(VecAXPY(D, 1.0, BFGS));
    113:   } else if (theta == 1.0) {
    114:     PetscCall(VecAXPY(D, 1.0 / yts, DFP));
    115:   } else {
    116:     /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
    117:     PetscCall(VecAXPBYPCZ(D, 1.0 - theta, theta / yts, 1.0, BFGS, DFP));
    118:   }
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: static PetscErrorCode SymBroydenRescaleUpdateDiagonal(Mat B, SymBroydenRescale ldb)
    123: {
    124:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
    125:   PetscInt    oldest, next;
    126:   Vec         invD, s_last, y_last;
    127:   LMBasis     S, Y;
    128:   PetscScalar yts;
    129:   PetscReal   sigma;

    131:   PetscFunctionBegin;
    132:   next   = ldb->k;
    133:   oldest = PetscMax(0, ldb->k - ldb->sigma_hist);
    134:   PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, next - 1, &yts));
    135:   PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_S, &S, NULL, NULL));
    136:   PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
    137:   PetscCall(LMBasisGetVecRead(S, next - 1, &s_last));
    138:   PetscCall(LMBasisGetVecRead(Y, next - 1, &y_last));
    139:   PetscCall(MatLMVMGetJ0InvDiag(B, &invD));
    140:   if (ldb->forward) {
    141:     /* We are doing diagonal scaling of the forward Hessian B */
    142:     /*  BFGS = DFP = inv(D); */
    143:     PetscCall(VecCopy(invD, ldb->invDnew));
    144:     PetscCall(VecReciprocal(ldb->invDnew));
    145:     PetscCall(DiagonalUpdate(ldb, ldb->invDnew, s_last, y_last, ldb->V, ldb->W, ldb->BFGS, ldb->DFP, ldb->theta, PetscRealPart(yts)));
    146:     /*  Obtain inverse and ensure positive definite */
    147:     PetscCall(VecReciprocal(ldb->invDnew));
    148:   } else {
    149:     /* Inverse Hessian update instead. */
    150:     PetscCall(VecCopy(invD, ldb->invDnew));
    151:     PetscCall(DiagonalUpdate(ldb, ldb->invDnew, y_last, s_last, ldb->V, ldb->W, ldb->DFP, ldb->BFGS, 1.0 - ldb->theta, PetscRealPart(yts)));
    152:   }
    153:   PetscCall(VecAbs(ldb->invDnew));
    154:   PetscCall(LMBasisRestoreVecRead(Y, next - 1, &y_last));
    155:   PetscCall(LMBasisRestoreVecRead(S, next - 1, &s_last));

    157:   if (ldb->sigma_hist > 0) {
    158:     // We are computing the scaling factor sigma that minimizes
    159:     //
    160:     // Sum_i || sigma^(alpha) (D^(-beta) o y_i) - sigma^(alpha-1) (D^(1-beta) o s_i) ||_2^2
    161:     //                        `-------.-------'                   `--------.-------'
    162:     //                               v_i                                  w_i
    163:     //
    164:     // To do this we first have to compute the sums of the dot product terms
    165:     //
    166:     // yy_sum = Sum_i v_i^T v_i,
    167:     // ys_sum = Sum_i v_i^T w_i, and
    168:     // ss_sum = Sum_i w_i^T w_i.
    169:     //
    170:     // These appear in the quadratic equation for the optimality condition for sigma,
    171:     //
    172:     // [alpha yy_sum] sigma^2 - [(2 alpha - 1) ys_sum] * sigma + [(alpha - 1) * ss_sum] = 0
    173:     //
    174:     // which we solve for sigma.

    176:     PetscReal yy_sum = 0; /*  No safeguard required */
    177:     PetscReal ys_sum = 0; /*  No safeguard required */
    178:     PetscReal ss_sum = 0; /*  No safeguard required */
    179:     PetscInt  start  = PetscMax(oldest, lmvm->k - ldb->sigma_hist);

    181:     Vec D_minus_beta             = NULL;
    182:     Vec D_minus_beta_squared     = NULL;
    183:     Vec D_one_minus_beta         = NULL;
    184:     Vec D_one_minus_beta_squared = NULL;
    185:     if (ldb->beta == 0.5) {
    186:       D_minus_beta_squared = ldb->invDnew; // (D^(-0.5))^2 = D^-1

    188:       PetscCall(VecCopy(ldb->invDnew, ldb->U));
    189:       PetscCall(VecReciprocal(ldb->U));
    190:       D_one_minus_beta_squared = ldb->U; // (D^(1-0.5))^2 = D
    191:     } else if (ldb->beta == 0.0) {
    192:       PetscCall(VecCopy(ldb->invDnew, ldb->U));
    193:       PetscCall(VecReciprocal(ldb->U));
    194:       D_one_minus_beta = ldb->U; // D^1
    195:     } else if (ldb->beta == 1.0) {
    196:       D_minus_beta = ldb->invDnew; // D^-1
    197:     } else {
    198:       PetscCall(VecCopy(ldb->invDnew, ldb->DFP));
    199:       PetscCall(VecPow(ldb->DFP, ldb->beta));
    200:       D_minus_beta = ldb->DFP;

    202:       PetscCall(VecCopy(ldb->invDnew, ldb->BFGS));
    203:       PetscCall(VecPow(ldb->BFGS, ldb->beta - 1));
    204:       D_one_minus_beta = ldb->BFGS;
    205:     }
    206:     for (PetscInt i = start - oldest; i < next - oldest; ++i) {
    207:       Vec s_i, y_i;

    209:       PetscCall(LMBasisGetVecRead(S, oldest + i, &s_i));
    210:       PetscCall(LMBasisGetVecRead(Y, oldest + i, &y_i));
    211:       if (ldb->beta == 0.5) {
    212:         PetscReal ytDinvy, stDs;

    214:         PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta_squared));
    215:         PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta_squared));
    216:         PetscCall(VecDotRealPart(ldb->W, s_i, &stDs));
    217:         PetscCall(VecDotRealPart(ldb->V, y_i, &ytDinvy));
    218:         ss_sum += stDs;        // ||s||_{D^(2*(1-beta))}^2
    219:         ys_sum += ldb->yts[i]; // s^T D^(1 - 2*beta) y
    220:         yy_sum += ytDinvy;     // ||y||_{D^(-2*beta)}^2
    221:       } else if (ldb->beta == 0.0) {
    222:         PetscScalar ytDs_scalar;
    223:         PetscReal   stDsr;

    225:         PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta));
    226:         PetscCall(VecDotNorm2(y_i, ldb->W, &ytDs_scalar, &stDsr));
    227:         ss_sum += stDsr;                      // ||s||_{D^(2*(1-beta))}^2
    228:         ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
    229:         yy_sum += ldb->yty[i];                // ||y||_{D^(-2*beta)}^2
    230:       } else if (ldb->beta == 1.0) {
    231:         PetscScalar ytDs_scalar;
    232:         PetscReal   ytDyr;

    234:         PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta));
    235:         PetscCall(VecDotNorm2(s_i, ldb->V, &ytDs_scalar, &ytDyr));
    236:         ss_sum += ldb->sts[i];                // ||s||_{D^(2*(1-beta))}^2
    237:         ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
    238:         yy_sum += ytDyr;                      // ||y||_{D^(-2*beta)}^2
    239:       } else {
    240:         PetscScalar ytDs_scalar;
    241:         PetscReal   ytDyr, stDs;

    243:         PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta));
    244:         PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta));
    245:         PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs_scalar, &ytDyr));
    246:         PetscCall(VecDotRealPart(ldb->W, ldb->W, &stDs));
    247:         ss_sum += stDs;                       // ||s||_{D^(2*(1-beta))}^2
    248:         ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
    249:         yy_sum += ytDyr;                      // ||y||_{D^(-2*beta)}^2
    250:       }
    251:       PetscCall(LMBasisRestoreVecRead(Y, oldest + i, &y_i));
    252:       PetscCall(LMBasisRestoreVecRead(S, oldest + i, &s_i));
    253:     }

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

    259:       sigma = ss_sum / ys_sum;
    260:     } else if (1.0 == ldb->alpha) {
    261:       /* yy_sum is never 0; if it were, we'd be at the minimum */
    262:       sigma = ys_sum / yy_sum;
    263:     } else {
    264:       PetscReal a         = ldb->alpha * yy_sum;
    265:       PetscReal b         = -(2.0 * ldb->alpha - 1.0) * ys_sum;
    266:       PetscReal c         = (ldb->alpha - 1.0) * ss_sum;
    267:       PetscReal sqrt_disc = PetscSqrtReal(b * b - 4 * a * c);

    269:       // numerically stable computation of positive root
    270:       if (b >= 0.0) {
    271:         if (a >= 0) {
    272:           PetscReal denom = PetscMax(-b - sqrt_disc, ldb->tol);

    274:           sigma = (2 * c) / denom;
    275:         } else {
    276:           PetscReal denom = PetscMax(2 * a, ldb->tol);

    278:           sigma = (-b - sqrt_disc) / denom;
    279:         }
    280:       } else {
    281:         if (a >= 0) {
    282:           PetscReal denom = PetscMax(2 * a, ldb->tol);

    284:           sigma = (-b + sqrt_disc) / denom;
    285:         } else {
    286:           PetscReal denom = PetscMax(-b + sqrt_disc, ldb->tol);

    288:           sigma = (2 * c) / denom;
    289:         }
    290:       }
    291:     }
    292:   } else {
    293:     sigma = 1.0;
    294:   }
    295:   /*  If Q has small values, then Q^(r_beta - 1)
    296:       can have very large values.  Hence, ys_sum
    297:       and ss_sum can be infinity.  In this case,
    298:       sigma can either be not-a-number or infinity. */
    299:   if (PetscIsNormalReal(sigma)) PetscCall(VecScale(ldb->invDnew, sigma));
    300:   /* Combine the old diagonal and the new diagonal using a convex limiter */
    301:   if (ldb->rho == 1.0) PetscCall(VecCopy(ldb->invDnew, invD));
    302:   else if (ldb->rho) PetscCall(VecAXPBY(invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
    303:   PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD));
    304:   PetscFunctionReturn(PETSC_SUCCESS);
    305: }

    307: static PetscErrorCode SymBroydenRescaleUpdateJ0(Mat B, SymBroydenRescale ldb)
    308: {
    309:   PetscFunctionBegin;
    310:   if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(SymBroydenRescaleUpdateScalar(B, ldb));
    311:   else if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(SymBroydenRescaleUpdateDiagonal(B, ldb));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: PETSC_INTERN PetscErrorCode SymBroydenRescaleUpdate(Mat B, SymBroydenRescale ldb)
    316: {
    317:   PetscInt oldest, next;

    319:   PetscFunctionBegin;
    320:   PetscCall(SymBroydenRescaleSetUp(B, ldb));
    321:   if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_NONE || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_USER) PetscFunctionReturn(PETSC_SUCCESS);
    322:   PetscCall(PetscLogEventBegin(SBRDN_Rescale, NULL, NULL, NULL, NULL));
    323:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    324:   if (next > ldb->k) {
    325:     PetscInt new_oldest = PetscMax(0, next - ldb->sigma_hist);
    326:     PetscInt ldb_oldest = PetscMax(0, ldb->k - ldb->sigma_hist);

    328:     if (new_oldest > ldb_oldest) {
    329:       for (PetscInt i = new_oldest; i < ldb->k; i++) {
    330:         ldb->yty[i - new_oldest] = ldb->yty[i - ldb_oldest];
    331:         ldb->yts[i - new_oldest] = ldb->yts[i - ldb_oldest];
    332:         ldb->sts[i - new_oldest] = ldb->sts[i - ldb_oldest];
    333:       }
    334:     }
    335:     for (PetscInt i = PetscMax(new_oldest, ldb->k); i < next; i++) {
    336:       PetscScalar yty, sts, yts;

    338:       PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_Y, i, &yty));
    339:       PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, i, &yts));
    340:       PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_S, LMBASIS_S, i, &sts));
    341:       ldb->yty[i - new_oldest] = PetscRealPart(yty);
    342:       ldb->yts[i - new_oldest] = PetscRealPart(yts);
    343:       ldb->sts[i - new_oldest] = PetscRealPart(sts);
    344:     }
    345:     ldb->k = next;
    346:   }
    347:   PetscCall(SymBroydenRescaleUpdateJ0(B, ldb));
    348:   PetscCall(PetscLogEventEnd(SBRDN_Rescale, NULL, NULL, NULL, NULL));
    349:   PetscFunctionReturn(PETSC_SUCCESS);
    350: }

    352: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDelta(Mat B, SymBroydenRescale ldb, PetscReal delta)
    353: {
    354:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
    355:   PetscBool same;

    357:   PetscFunctionBegin;
    358:   same       = (delta == ldb->delta) ? PETSC_TRUE : PETSC_FALSE;
    359:   ldb->delta = delta;
    360:   ldb->delta = PetscMin(ldb->delta, ldb->delta_max);
    361:   ldb->delta = PetscMax(ldb->delta, ldb->delta_min);
    362:   // if there have been no updates yet, propagate delta to J0
    363:   if (!same && !lmvm->prev_set && (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL)) {
    364:     ldb->initialized = PETSC_FALSE;
    365:     B->preallocated  = PETSC_FALSE; // make sure SyBroydenInitializeJ0() is called in the next MatCheckPreallocated()
    366:   }
    367:   PetscFunctionReturn(PETSC_SUCCESS);
    368: }

    370: PETSC_INTERN PetscErrorCode SymBroydenRescaleCopy(SymBroydenRescale bctx, SymBroydenRescale mctx)
    371: {
    372:   PetscFunctionBegin;
    373:   mctx->scale_type = bctx->scale_type;
    374:   mctx->theta      = bctx->theta;
    375:   mctx->alpha      = bctx->alpha;
    376:   mctx->beta       = bctx->beta;
    377:   mctx->rho        = bctx->rho;
    378:   mctx->delta      = bctx->delta;
    379:   mctx->delta_min  = bctx->delta_min;
    380:   mctx->delta_max  = bctx->delta_max;
    381:   mctx->tol        = bctx->tol;
    382:   mctx->sigma_hist = bctx->sigma_hist;
    383:   mctx->forward    = bctx->forward;
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDiagonalMode(SymBroydenRescale ldb, PetscBool forward)
    388: {
    389:   PetscFunctionBegin;
    390:   ldb->forward = forward;
    391:   PetscFunctionReturn(PETSC_SUCCESS);
    392: }

    394: PETSC_INTERN PetscErrorCode SymBroydenRescaleGetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType *stype)
    395: {
    396:   PetscFunctionBegin;
    397:   *stype = ldb->scale_type;
    398:   PetscFunctionReturn(PETSC_SUCCESS);
    399: }

    401: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType stype)
    402: {
    403:   PetscFunctionBegin;
    404:   ldb->scale_type = stype;
    405:   PetscFunctionReturn(PETSC_SUCCESS);
    406: }

    408: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetFromOptions(Mat B, SymBroydenRescale ldb, PetscOptionItems PetscOptionsObject)
    409: {
    410:   MatLMVMSymBroydenScaleType stype = ldb->scale_type;
    411:   PetscBool                  flg;

    413:   PetscFunctionBegin;
    414:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for updating diagonal Jacobian approximation (MATLMVMDIAGBRDN)");
    415:   PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBroydenScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
    416:   PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL));
    417:   PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL));
    418:   PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
    419:   PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL, 0.0, 1.0));
    420:   PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
    421:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
    422:   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));
    423:   PetscOptionsHeadEnd();
    424:   PetscCheck(!(ldb->theta < 0.0) && !(ldb->theta > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
    425:   PetscCheck(!(ldb->alpha < 0.0) && !(ldb->alpha > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
    426:   PetscCheck(!(ldb->rho < 0.0) && !(ldb->rho > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
    427:   PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
    428:   if (flg) PetscCall(SymBroydenRescaleSetType(ldb, stype));
    429:   PetscFunctionReturn(PETSC_SUCCESS);
    430: }

    432: static PetscErrorCode SymBroydenRescaleAllocate(Mat B, SymBroydenRescale ldb)
    433: {
    434:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    436:   PetscFunctionBegin;
    437:   if (!ldb->allocated) {
    438:     PetscCall(PetscMalloc3(ldb->sigma_hist, &ldb->yty, ldb->sigma_hist, &ldb->yts, ldb->sigma_hist, &ldb->sts));
    439:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
    440:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
    441:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
    442:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
    443:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
    444:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
    445:     ldb->allocated = PETSC_TRUE;
    446:   }
    447:   PetscFunctionReturn(PETSC_SUCCESS);
    448: }

    450: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetUp(Mat B, SymBroydenRescale ldb)
    451: {
    452:   PetscFunctionBegin;
    453:   PetscCall(SymBroydenRescaleAllocate(B, ldb));
    454:   if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DECIDE) {
    455:     Mat       J0;
    456:     PetscBool is_constant_or_diagonal;

    458:     PetscCall(MatLMVMGetJ0(B, &J0));
    459:     PetscCall(PetscObjectTypeCompareAny((PetscObject)J0, &is_constant_or_diagonal, MATCONSTANTDIAGONAL, MATDIAGONAL, ""));
    460:     ldb->scale_type = is_constant_or_diagonal ? MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL : MAT_LMVM_SYMBROYDEN_SCALE_NONE;
    461:   }
    462:   PetscFunctionReturn(PETSC_SUCCESS);
    463: }

    465: PETSC_INTERN PetscErrorCode SymBroydenRescaleInitializeJ0(Mat B, SymBroydenRescale ldb)
    466: {
    467:   PetscFunctionBegin;
    468:   if (!ldb->initialized) {
    469:     PetscCall(SymBroydenRescaleSetUp(B, ldb));
    470:     switch (ldb->scale_type) {
    471:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL: {
    472:       Vec invD;

    474:       PetscCall(MatLMVMGetJ0InvDiag(B, &invD));
    475:       PetscCall(VecSet(invD, ldb->delta));
    476:       PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD));
    477:       break;
    478:     }
    479:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
    480:       PetscCall(MatLMVMSetJ0Scale(B, 1.0 / ldb->delta));
    481:       break;
    482:     default:
    483:       break;
    484:     }
    485:     ldb->initialized = PETSC_TRUE;
    486:   }
    487:   PetscFunctionReturn(PETSC_SUCCESS);
    488: }

    490: PETSC_INTERN PetscErrorCode SymBroydenRescaleView(SymBroydenRescale ldb, PetscViewer pv)
    491: {
    492:   PetscFunctionBegin;
    493:   PetscBool isascii;
    494:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
    495:   if (isascii) {
    496:     PetscCall(PetscViewerASCIIPrintf(pv, "Rescale type: %s\n", MatLMVMSymBroydenScaleTypes[ldb->scale_type]));
    497:     if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
    498:       PetscCall(PetscViewerASCIIPrintf(pv, "Rescale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
    499:       PetscCall(PetscViewerASCIIPrintf(pv, "Rescale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
    500:     }
    501:     if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(PetscViewerASCIIPrintf(pv, "Rescale convex factor: theta=%g\n", (double)ldb->theta));
    502:   }
    503:   PetscFunctionReturn(PETSC_SUCCESS);
    504: }

    506: PETSC_INTERN PetscErrorCode SymBroydenRescaleReset(Mat B, SymBroydenRescale ldb, MatLMVMResetMode mode)
    507: {
    508:   PetscFunctionBegin;
    509:   ldb->k = 0;
    510:   if (MatLMVMResetClearsBases(mode)) {
    511:     if (ldb->allocated) {
    512:       PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
    513:       PetscCall(VecDestroy(&ldb->invDnew));
    514:       PetscCall(VecDestroy(&ldb->BFGS));
    515:       PetscCall(VecDestroy(&ldb->DFP));
    516:       PetscCall(VecDestroy(&ldb->U));
    517:       PetscCall(VecDestroy(&ldb->V));
    518:       PetscCall(VecDestroy(&ldb->W));
    519:       ldb->allocated = PETSC_FALSE;
    520:     }
    521:   }
    522:   if (B && ldb->initialized && !MatLMVMResetClearsAll(mode)) PetscCall(SymBroydenRescaleInitializeJ0(B, ldb)); // eagerly reset J0 if we are rescaling
    523:   PetscFunctionReturn(PETSC_SUCCESS);
    524: }

    526: PETSC_INTERN PetscErrorCode SymBroydenRescaleDestroy(SymBroydenRescale *ldb)
    527: {
    528:   PetscFunctionBegin;
    529:   PetscCall(SymBroydenRescaleReset(NULL, *ldb, MAT_LMVM_RESET_ALL));
    530:   PetscCall(PetscFree(*ldb));
    531:   *ldb = NULL;
    532:   PetscFunctionReturn(PETSC_SUCCESS);
    533: }

    535: PETSC_INTERN PetscErrorCode SymBroydenRescaleCreate(SymBroydenRescale *ldb)
    536: {
    537:   PetscFunctionBegin;
    538:   PetscCall(PetscNew(ldb));
    539:   (*ldb)->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DECIDE;
    540:   (*ldb)->theta      = 0.0;
    541:   (*ldb)->alpha      = 1.0;
    542:   (*ldb)->rho        = 1.0;
    543:   (*ldb)->forward    = PETSC_TRUE;
    544:   (*ldb)->beta       = 0.5;
    545:   (*ldb)->delta      = 1.0;
    546:   (*ldb)->delta_min  = 1e-7;
    547:   (*ldb)->delta_max  = 100.0;
    548:   (*ldb)->tol        = 1e-8;
    549:   (*ldb)->sigma_hist = 1;
    550:   (*ldb)->allocated  = PETSC_FALSE;
    551:   PetscFunctionReturn(PETSC_SUCCESS);
    552: }