Actual source code: lmvmutils.c

      1: #include <petscdevice.h>
      2: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
      3: #include <petsc/private/deviceimpl.h>
      4: #include <petscblaslapack.h>

      6: /*@
      7:   MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to a `MATLMVM` matrix.

      9:   Input Parameters:
     10: + B - A `MATLMVM` matrix
     11: . X - Solution vector
     12: - F - Function vector

     14:   Level: intermediate

     16:   Notes:

     18:   The first time this function is called for a `MATLMVM` matrix, no update is applied, but the given X and F vectors
     19:   are stored for use as Xprev and Fprev in the next update.

     21:   If the user has provided another `MATLMVM` matrix for the reference Jacobian (using `MatLMVMSetJ0()`, for example),
     22:   that matrix is also updated recursively.

     24:   If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMUpdate()` is
     25:   called, the row size and layout of `B` will be set to match `F` and the column size and layout of `B` will be set to
     26:   match `X`, and these sizes will be final.

     28: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMAllocate()`
     29: @*/
     30: PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
     31: {
     32:   Mat_LMVM *lmvm;
     33:   PetscBool same;

     35:   PetscFunctionBegin;
     39:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
     40:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
     41:   /* If B has specified layouts, this will check X and F are compatible;
     42:      if B does not have specified layouts, this will adopt them, so that
     43:      this pattern is okay

     45:        MatCreate(comm, &B);
     46:        MatLMVMSetType(B, MATLMVMBFGS);
     47:        MatLMVMUpdate(B, X, F);
     48:    */
     49:   PetscCall(MatLMVMUseVecLayoutsIfCompatible(B, X, F));
     50:   MatCheckPreallocated(B, 1);
     51:   PetscCall(PetscLogEventBegin(MATLMVM_Update, NULL, NULL, NULL, NULL));
     52:   lmvm = (Mat_LMVM *)B->data;
     53:   PetscCall(MatLMVMUpdate(lmvm->J0, X, F));
     54:   PetscCall((*lmvm->ops->update)(B, X, F));
     55:   PetscCall(PetscLogEventEnd(MATLMVM_Update, NULL, NULL, NULL, NULL));
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode MatLMVMCreateJ0(Mat B, Mat *J0)
     60: {
     61:   PetscLayout rmap, cmap;
     62:   VecType     vec_type;
     63:   const char *prefix;

     65:   PetscFunctionBegin;
     66:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), J0));
     67:   PetscCall(MatGetLayouts(B, &rmap, &cmap));
     68:   PetscCall(MatSetLayouts(*J0, rmap, cmap));
     69:   PetscCall(MatGetVecType(B, &vec_type));
     70:   PetscCall(MatSetVecType(*J0, vec_type));
     71:   PetscCall(MatGetOptionsPrefix(B, &prefix));
     72:   PetscCall(MatSetOptionsPrefix(*J0, prefix));
     73:   PetscCall(MatAppendOptionsPrefix(*J0, "mat_lmvm_J0_"));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: static PetscErrorCode MatLMVMCreateJ0KSP(Mat B, KSP *ksp)
     78: {
     79:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
     80:   const char *prefix;

     82:   PetscFunctionBegin;
     83:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)B), ksp));
     84:   PetscCall(KSPSetOperators(*ksp, lmvm->J0, lmvm->J0));
     85:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)B, (PetscObject)*ksp, 1));
     86:   PetscCall(MatGetOptionsPrefix(B, &prefix));
     87:   PetscCall(KSPSetOptionsPrefix(*ksp, prefix));
     88:   PetscCall(KSPAppendOptionsPrefix(*ksp, "mat_lmvm_J0_"));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: static PetscErrorCode MatLMVMCreateJ0KSP_ExactInverse(Mat B, KSP *ksp)
     93: {
     94:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
     95:   PC        pc;

     97:   PetscFunctionBegin;
     98:   PetscCall(MatLMVMCreateJ0KSP(B, ksp));
     99:   PetscCall(KSPSetType(*ksp, KSPPREONLY));
    100:   PetscCall(KSPGetPC(*ksp, &pc));
    101:   PetscCall(PCSetType(pc, PCMAT));
    102:   PetscCall(PCMatSetApplyOperation(pc, MATOP_SOLVE));
    103:   lmvm->disable_ksp_viewers = PETSC_TRUE;
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: /*@
    108:   MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
    109:   an identity matrix (scale = 1.0).

    111:   Input Parameter:
    112: . B - A `MATLMVM` matrix

    114:   Level: advanced

    116: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
    117: @*/
    118: PetscErrorCode MatLMVMClearJ0(Mat B)
    119: {
    120:   Mat_LMVM *lmvm;
    121:   PetscBool same;

    123:   PetscFunctionBegin;
    125:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    126:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    127:   lmvm = (Mat_LMVM *)B->data;
    128:   PetscCall(MatDestroy(&lmvm->J0));
    129:   PetscCall(KSPDestroy(&lmvm->J0ksp));
    130:   PetscCall(MatLMVMCreateJ0(B, &lmvm->J0));
    131:   PetscCall(MatSetType(lmvm->J0, MATCONSTANTDIAGONAL));
    132:   PetscCall(MatZeroEntries(lmvm->J0));
    133:   PetscCall(MatShift(lmvm->J0, 1.0));
    134:   PetscCall(MatLMVMCreateJ0KSP_ExactInverse(B, &lmvm->J0ksp));
    135:   lmvm->created_J0    = PETSC_TRUE;
    136:   lmvm->created_J0ksp = PETSC_TRUE;
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: /*@
    141:   MatLMVMSetJ0Scale - Allows the user to define a scalar value
    142:   mu such that J0 = mu*I.

    144:   Input Parameters:
    145: + B     - A `MATLMVM` matrix
    146: - scale - Scalar value mu that defines the initial Jacobian

    148:   Level: advanced

    150: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
    151: @*/
    152: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
    153: {
    154:   Mat_LMVM *lmvm;
    155:   PetscBool same;
    156:   PetscBool isconstant;

    158:   PetscFunctionBegin;
    160:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    161:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    162:   lmvm = (Mat_LMVM *)B->data;
    163:   PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATCONSTANTDIAGONAL, &isconstant));
    164:   if (!isconstant) PetscCall(MatLMVMClearJ0(B));
    165:   PetscCall(MatZeroEntries(lmvm->J0));
    166:   PetscCall(MatShift(lmvm->J0, scale));
    167:   PetscFunctionReturn(PETSC_SUCCESS);
    168: }

    170: #if PetscDefined(USE_DEBUG)
    172:     do { \
    173:       if (!(layout)->setupcalled) { \
    174:         PetscMPIInt global[2]; \
    175:         global[0] = (PetscMPIInt)(v); \
    176:         global[1] = -global[0]; \
    177:         PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &global[0], 2, MPI_INT, MPI_MIN, ((layout)->comm))); \
    178:         PetscCheck(global[1] == -global[0], ((layout)->comm), PETSC_ERR_ARG_WRONGSTATE, "PetscLayout has size == PETSC_DECIDE and local size == PETSC_DETERMINE on only some processes"); \
    179:       } \
    180:     } while (0)
    181: #else
    183:     do { \
    184:       (void)(comm); \
    185:       (void)(v); \
    186:     } while (0)
    187: #endif

    189: static PetscErrorCode MatLMVMCheckArgumentLayout(PetscLayout b, PetscLayout a)
    190: {
    191:   PetscBool   b_is_unspecified, a_is_specified, are_compatible;
    192:   PetscLayout b_setup = NULL, a_setup = NULL;

    194:   PetscFunctionBegin;
    195:   if (b == a) PetscFunctionReturn(PETSC_SUCCESS); // a layout is compatible with itself
    196:   if (b->setupcalled && a->setupcalled) {
    197:     // run the standard checks that are guaranteed to error on at least one process if the layouts are incompatible
    198:     PetscCheck(b->N == a->N, b->comm, PETSC_ERR_ARG_SIZ, "argument layout (size %" PetscInt_FMT ") is incompatible with MatLMVM layout (size %" PetscInt_FMT ")", a->N, b->N);
    199:     PetscCheck(b->n == a->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "argument layout (local size %" PetscInt_FMT ") is incompatible with MatLMVM layout (local size %" PetscInt_FMT ")", a->n, b->n);
    200:     PetscFunctionReturn(PETSC_SUCCESS);
    201:   }
    202:   a_is_specified = (a->n >= 0) || (a->N >= 0) ? PETSC_TRUE : PETSC_FALSE;
    204:   PetscCheck(a_is_specified, a->comm, PETSC_ERR_ARG_WRONGSTATE, "argument layout has n == PETSC_DETERMINE and N == PETSC_DECIDE, size must be specified first");
    205:   b_is_unspecified = (b->n < 0) && (b->N < 0) ? PETSC_TRUE : PETSC_FALSE;
    207:   if (b_is_unspecified) PetscFunctionReturn(PETSC_SUCCESS); // any layout can replace an unspecified layout
    208:   // we don't want to change the setup states in this check, so make duplicates if they have not been setup
    209:   if (!b->setupcalled) {
    210:     PetscCall(PetscLayoutDuplicate(b, &b_setup));
    211:     PetscCall(PetscLayoutSetUp(b_setup));
    212:   } else PetscCall(PetscLayoutReference(b, &b_setup));
    213:   if (!a->setupcalled) {
    214:     PetscCall(PetscLayoutDuplicate(a, &a_setup));
    215:     PetscCall(PetscLayoutSetUp(a_setup));
    216:   } else PetscCall(PetscLayoutReference(a, &a_setup));
    217:   PetscCall(PetscLayoutCompare(b_setup, a_setup, &are_compatible));
    218:   PetscCall(PetscLayoutDestroy(&a_setup));
    219:   PetscCall(PetscLayoutDestroy(&b_setup));
    220:   PetscCheck(are_compatible, b->comm, PETSC_ERR_ARG_SIZ, "argument layout (size %" PetscInt_FMT ") is incompatible with MatLMVM layout (size %" PetscInt_FMT ")", a->N, b->N);
    221:   PetscFunctionReturn(PETSC_SUCCESS);
    222: }

    224: static PetscErrorCode MatLMVMUseJ0LayoutsIfCompatible(Mat B, Mat J0)
    225: {
    226:   PetscFunctionBegin;
    227:   PetscCall(MatLMVMCheckArgumentLayout(B->rmap, J0->rmap));
    228:   PetscCall(MatLMVMCheckArgumentLayout(B->cmap, J0->cmap));
    229:   PetscCall(PetscLayoutSetUp(J0->rmap));
    230:   PetscCall(PetscLayoutSetUp(J0->cmap));
    231:   PetscCall(PetscLayoutReference(J0->rmap, &B->rmap));
    232:   PetscCall(PetscLayoutReference(J0->cmap, &B->cmap));
    233:   PetscFunctionReturn(PETSC_SUCCESS);
    234: }

    236: static PetscErrorCode MatLMVMUseJ0DiagLayoutsIfCompatible(Mat B, Vec J0_diag)
    237: {
    238:   PetscFunctionBegin;
    239:   PetscCall(MatLMVMCheckArgumentLayout(B->rmap, J0_diag->map));
    240:   PetscCall(MatLMVMCheckArgumentLayout(B->cmap, J0_diag->map));
    241:   PetscCall(PetscLayoutSetUp(J0_diag->map));
    242:   PetscCall(PetscLayoutReference(J0_diag->map, &B->rmap));
    243:   PetscCall(PetscLayoutReference(J0_diag->map, &B->cmap));
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

    247: PETSC_INTERN PetscErrorCode MatLMVMUseVecLayoutsIfCompatible(Mat B, Vec X, Vec F)
    248: {
    249:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    251:   PetscFunctionBegin;
    252:   PetscCall(MatLMVMCheckArgumentLayout(B->rmap, F->map));
    253:   PetscCall(MatLMVMCheckArgumentLayout(B->cmap, X->map));
    254:   PetscCall(PetscLayoutSetUp(F->map));
    255:   PetscCall(PetscLayoutSetUp(X->map));
    256:   PetscCall(PetscLayoutReference(F->map, &B->rmap));
    257:   PetscCall(PetscLayoutReference(X->map, &B->cmap));
    258:   if (lmvm->created_J0) {
    259:     PetscCall(PetscLayoutReference(B->rmap, &lmvm->J0->rmap));
    260:     PetscCall(PetscLayoutReference(B->cmap, &lmvm->J0->cmap));
    261:   }
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }

    265: /*@
    266:   MatLMVMSetJ0Diag - Allows the user to define a vector
    267:   V such that J0 = diag(V).

    269:   Input Parameters:
    270: + B - An LMVM-type matrix
    271: - V - Vector that defines the diagonal of the initial Jacobian: values are copied, V is not referenced

    273:   Level: advanced

    275:   Note:
    276:   If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0Diag()` is
    277:   called, the rows and columns of `B` will each have the size and layout of `V`, and these sizes will be final.

    279: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
    280: @*/
    281: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
    282: {
    283:   Mat       J0diag;
    284:   PetscBool same;
    285:   VecType   vec_type;

    287:   PetscFunctionBegin;
    290:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    291:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    292:   PetscCheckSameComm(B, 1, V, 2);
    293:   PetscCall(MatLMVMUseJ0DiagLayoutsIfCompatible(B, V));
    294:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &J0diag));
    295:   PetscCall(MatSetLayouts(J0diag, V->map, V->map));
    296:   PetscCall(VecGetType(V, &vec_type));
    297:   PetscCall(MatSetVecType(J0diag, vec_type));
    298:   PetscCall(MatSetType(J0diag, MATDIAGONAL));
    299:   PetscCall(MatDiagonalSet(J0diag, V, INSERT_VALUES));
    300:   PetscCall(MatLMVMSetJ0(B, J0diag));
    301:   PetscCall(MatDestroy(&J0diag));
    302:   PetscFunctionReturn(PETSC_SUCCESS);
    303: }

    305: PETSC_INTERN PetscErrorCode MatLMVMGetJ0InvDiag(Mat B, Vec *V)
    306: {
    307:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
    308:   PetscBool isvdiag;

    310:   PetscFunctionBegin;
    311:   PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATDIAGONAL, &isvdiag));
    312:   if (!isvdiag) {
    313:     PetscCall(MatLMVMClearJ0(B));
    314:     PetscCall(MatSetType(lmvm->J0, MATDIAGONAL));
    315:     PetscCall(MatZeroEntries(lmvm->J0));
    316:     PetscCall(MatShift(lmvm->J0, 1.0));
    317:   }
    318:   PetscCall(MatDiagonalGetInverseDiagonal(lmvm->J0, V));
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: PETSC_INTERN PetscErrorCode MatLMVMRestoreJ0InvDiag(Mat B, Vec *V)
    323: {
    324:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    326:   PetscFunctionBegin;
    327:   PetscCall(MatDiagonalRestoreInverseDiagonal(lmvm->J0, V));
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }

    331: PETSC_INTERN PetscErrorCode MatLMVMJ0KSPIsExact(Mat B, PetscBool *is_exact)
    332: {
    333:   PetscBool    is_preonly, is_pcmat, has_pmat;
    334:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    335:   Mat          pc_pmat;
    336:   PC           pc;
    337:   MatOperation matop;

    339:   PetscFunctionBegin;
    340:   *is_exact = PETSC_FALSE;
    341:   PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0ksp, KSPPREONLY, &is_preonly));
    342:   if (!is_preonly) PetscFunctionReturn(PETSC_SUCCESS);
    343:   PetscCall(KSPGetPC(lmvm->J0ksp, &pc));
    344:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMAT, &is_pcmat));
    345:   if (!is_pcmat) PetscFunctionReturn(PETSC_SUCCESS);
    346:   PetscCall(PCGetOperatorsSet(pc, NULL, &has_pmat));
    347:   if (!has_pmat) PetscFunctionReturn(PETSC_SUCCESS);
    348:   PetscCall(PCGetOperators(pc, NULL, &pc_pmat));
    349:   if (pc_pmat != lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
    350:   PetscCall(PCMatGetApplyOperation(pc, &matop));
    351:   *is_exact = (matop == MATOP_SOLVE) ? PETSC_TRUE : PETSC_FALSE;
    352:   PetscFunctionReturn(PETSC_SUCCESS);
    353: }

    355: /*@
    356:   MatLMVMSetJ0 - Allows the user to define the initial Jacobian matrix from which the LMVM-type approximation is built
    357:   up.

    359:   Input Parameters:
    360: + B  - An LMVM-type matrix
    361: - J0 - The initial Jacobian matrix, will be referenced by B.

    363:   Level: advanced

    365:   Notes:
    366:   A KSP is created for inverting J0 with prefix `-mat_lmvm_J0_` and J0 is set to both operators in `KSPSetOperators()`.
    367:   If you want to use a separate preconditioning matrix, use `MatLMVMSetJ0KSP()` directly.

    369:   If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0()` is
    370:   called, then `B` will adopt the sizes and layouts of `J0`, and these sizes will be final.

    372: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
    373: @*/
    374: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
    375: {
    376:   Mat_LMVM *lmvm;
    377:   PetscBool same;
    378:   PetscBool J0_has_solve;

    380:   PetscFunctionBegin;
    383:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    384:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    385:   lmvm = (Mat_LMVM *)B->data;
    386:   if (J0 == lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
    387:   PetscCheckSameComm(B, 1, J0, 2);
    388:   PetscCall(MatLMVMUseJ0LayoutsIfCompatible(B, J0));
    389:   PetscCall(PetscObjectReference((PetscObject)J0));
    390:   PetscCall(MatDestroy(&lmvm->J0));
    391:   lmvm->J0         = J0;
    392:   lmvm->created_J0 = PETSC_FALSE;
    393:   PetscCall(MatHasOperation(J0, MATOP_SOLVE, &J0_has_solve));
    394:   if (J0_has_solve) {
    395:     PetscCall(KSPDestroy(&lmvm->J0ksp));
    396:     PetscCall(MatLMVMCreateJ0KSP_ExactInverse(B, &lmvm->J0ksp));
    397:     lmvm->created_J0ksp = PETSC_TRUE;
    398:   } else {
    399:     if (lmvm->created_J0ksp) {
    400:       PetscBool is_preonly, is_pcmat = PETSC_FALSE, is_pcmat_solve = PETSC_FALSE;
    401:       PC        pc;

    403:       PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0ksp, KSPPREONLY, &is_preonly));
    404:       PetscCall(KSPGetPC(lmvm->J0ksp, &pc));
    405:       if (pc) {
    406:         PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMAT, &is_pcmat));
    407:         if (is_pcmat) {
    408:           MatOperation matop;

    410:           PetscCall(PCMatGetApplyOperation(pc, &matop));
    411:           if (matop == MATOP_SOLVE) is_pcmat_solve = PETSC_TRUE;
    412:         }
    413:       }
    414:       if (is_preonly && is_pcmat_solve) {
    415:         /* The KSP is one created by LMVM for a mat that has a MatSolve() implementation.  Because this new J0 doesn't, change it to
    416:            a default KSP */
    417:         PetscCall(KSPDestroy(&lmvm->J0ksp));
    418:         PetscCall(MatLMVMCreateJ0KSP(B, &lmvm->J0ksp));
    419:       }
    420:     }
    421:     PetscCall(KSPSetOperators(lmvm->J0ksp, J0, J0));
    422:   }
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: /*@
    427:   MatLMVMSetJ0PC - Allows the user to define a `PC` object that acts as the initial inverse-Jacobian matrix.

    429:   Input Parameters:
    430: + B    - A `MATLMVM` matrix
    431: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0

    433:   Level: advanced

    435:   Notes:
    436:   `J0pc` should already contain all the operators necessary for its application.  The `MATLMVM` matrix only calls
    437:   `PCApply()` without changing any other options.

    439:   If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0PC()` is
    440:   called, then `B` will adopt the sizes and layouts of the operators of `J0pc`, and these sizes will be final.

    442: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
    443: @*/
    444: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
    445: {
    446:   Mat_LMVM *lmvm;
    447:   PetscBool same, mat_set, pmat_set;
    448:   PC        current_pc;
    449:   Mat       J0 = NULL;

    451:   PetscFunctionBegin;
    454:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    455:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    456:   lmvm = (Mat_LMVM *)B->data;
    457:   PetscCall(PCGetOperatorsSet(J0pc, &mat_set, &pmat_set));
    458:   PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
    459:   if (mat_set) PetscCall(PCGetOperators(J0pc, &J0, NULL));
    460:   else PetscCall(PCGetOperators(J0pc, NULL, &J0));
    461:   PetscCall(KSPGetPC(lmvm->J0ksp, &current_pc));
    462:   if (J0pc == current_pc && J0 == lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
    463:   PetscCall(MatLMVMSetJ0(B, J0));
    464:   PetscCall(KSPSetPC(lmvm->J0ksp, J0pc));
    465:   PetscFunctionReturn(PETSC_SUCCESS);
    466: }

    468: /*@
    469:   MatLMVMSetJ0KSP - Allows the user to provide a pre-configured KSP solver for the initial inverse-Jacobian
    470:   approximation.

    472:   Input Parameters:
    473: + B     - A `MATLMVM` matrix
    474: - J0ksp - `KSP` solver for the initial inverse-Jacobian application

    476:   Level: advanced

    478:   Note:
    479:   The `KSP` solver should already contain all the operators necessary to perform the inversion. The `MATLMVM` matrix
    480:   only calls `KSPSolve()` without changing any other options.

    482:   If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0KSP()` is
    483:   called, then `B` will adopt the sizes and layouts of the operators of `J0ksp`, and these sizes will be final.

    485: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
    486: @*/
    487: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
    488: {
    489:   Mat_LMVM *lmvm;
    490:   PetscBool same, mat_set, pmat_set;
    491:   Mat       J0;

    493:   PetscFunctionBegin;
    496:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    497:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    498:   lmvm = (Mat_LMVM *)B->data;
    499:   PetscCall(KSPGetOperatorsSet(J0ksp, &mat_set, &pmat_set));
    500:   PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
    501:   if (mat_set) PetscCall(KSPGetOperators(J0ksp, &J0, NULL));
    502:   else PetscCall(KSPGetOperators(J0ksp, NULL, &J0));
    503:   if (J0ksp == lmvm->J0ksp && lmvm->J0 == J0) PetscFunctionReturn(PETSC_SUCCESS);
    504:   PetscCall(MatLMVMSetJ0(B, J0));
    505:   if (J0ksp != lmvm->J0ksp) {
    506:     lmvm->created_J0ksp       = PETSC_FALSE;
    507:     lmvm->disable_ksp_viewers = PETSC_FALSE; // if the user supplies a more complicated KSP, don't turn off viewers
    508:   }
    509:   PetscCall(PetscObjectReference((PetscObject)J0ksp));
    510:   PetscCall(KSPDestroy(&lmvm->J0ksp));
    511:   lmvm->J0ksp = J0ksp;
    512:   PetscFunctionReturn(PETSC_SUCCESS);
    513: }

    515: /*@
    516:   MatLMVMGetJ0 - Returns a pointer to the internal `J0` matrix.

    518:   Input Parameter:
    519: . B - A `MATLMVM` matrix

    521:   Output Parameter:
    522: . J0 - `Mat` object for defining the initial Jacobian

    524:   Level: advanced

    526:   Note:

    528:   If `J0` was created by `B` it will have the options prefix `-mat_lmvm_J0_`.

    530: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
    531: @*/
    532: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
    533: {
    534:   Mat_LMVM *lmvm;
    535:   PetscBool same;

    537:   PetscFunctionBegin;
    539:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    540:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    541:   lmvm = (Mat_LMVM *)B->data;
    542:   *J0  = lmvm->J0;
    543:   PetscFunctionReturn(PETSC_SUCCESS);
    544: }

    546: /*@
    547:   MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
    548:   associated with the initial Jacobian.

    550:   Input Parameter:
    551: . B - A `MATLMVM` matrix

    553:   Output Parameter:
    554: . J0pc - `PC` object for defining the initial inverse-Jacobian

    556:   Level: advanced

    558: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
    559: @*/
    560: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
    561: {
    562:   Mat_LMVM *lmvm;
    563:   PetscBool same;

    565:   PetscFunctionBegin;
    567:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    568:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    569:   lmvm = (Mat_LMVM *)B->data;
    570:   PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
    571:   PetscFunctionReturn(PETSC_SUCCESS);
    572: }

    574: /*@
    575:   MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
    576:   associated with the initial Jacobian.

    578:   Input Parameter:
    579: . B - A `MATLMVM` matrix

    581:   Output Parameter:
    582: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian

    584:   Level: advanced

    586: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
    587: @*/
    588: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
    589: {
    590:   Mat_LMVM *lmvm;
    591:   PetscBool same;

    593:   PetscFunctionBegin;
    595:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    596:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    597:   lmvm   = (Mat_LMVM *)B->data;
    598:   *J0ksp = lmvm->J0ksp;
    599:   PetscFunctionReturn(PETSC_SUCCESS);
    600: }

    602: /*@
    603:   MatLMVMApplyJ0Fwd - Applies an approximation of the forward
    604:   matrix-vector product with the initial Jacobian.

    606:   Input Parameters:
    607: + B - A `MATLMVM` matrix
    608: - X - vector to multiply with J0

    610:   Output Parameter:
    611: . Y - resulting vector for the operation

    613:   Level: advanced

    615: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
    616:           `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
    617: @*/
    618: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
    619: {
    620:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    622:   PetscFunctionBegin;
    623:   PetscCall(MatMult(lmvm->J0, X, Y));
    624:   PetscFunctionReturn(PETSC_SUCCESS);
    625: }

    627: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0HermitianTranspose(Mat B, Vec X, Vec Y)
    628: {
    629:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    631:   PetscFunctionBegin;
    632:   PetscCall(MatMultHermitianTranspose(lmvm->J0, X, Y));
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: /*@
    637:   MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
    638:   inverse to the given vector.

    640:   Input Parameters:
    641: + B - A `MATLMVM` matrix
    642: - X - vector to "multiply" with J0^{-1}

    644:   Output Parameter:
    645: . Y - resulting vector for the operation

    647:   Level: advanced

    649:   Note:
    650:   The specific form of the application
    651:   depends on whether the user provided a scaling factor, a J0 matrix,
    652:   a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
    653:   provided, the function simply does an identity matrix application
    654:   (vector copy).

    656: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
    657:           `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
    658: @*/
    659: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
    660: {
    661:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    663:   PetscFunctionBegin;
    664:   if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
    665:   PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
    666:   if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
    667:   PetscFunctionReturn(PETSC_SUCCESS);
    668: }

    670: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvTranspose(Mat B, Vec X, Vec Y)
    671: {
    672:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    674:   PetscFunctionBegin;
    675:   if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
    676:   PetscCall(KSPSolveTranspose(lmvm->J0ksp, X, Y));
    677:   if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
    678:   PetscFunctionReturn(PETSC_SUCCESS);
    679: }

    681: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvHermitianTranspose(Mat B, Vec X, Vec Y)
    682: {
    683:   PetscFunctionBegin;
    684:   if (!PetscDefined(USE_COMPLEX)) {
    685:     PetscCall(MatLMVMApplyJ0InvTranspose(B, X, Y));
    686:   } else {
    687:     Vec X_conj;

    689:     PetscCall(VecDuplicate(X, &X_conj));
    690:     PetscCall(VecCopy(X, X_conj));
    691:     PetscCall(VecConjugate(X_conj));
    692:     PetscCall(MatLMVMApplyJ0InvTranspose(B, X_conj, Y));
    693:     PetscCall(VecConjugate(Y));
    694:     PetscCall(VecDestroy(&X_conj));
    695:   }
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: /*@
    700:   MatLMVMIsAllocated - Returns a boolean flag that shows whether
    701:   the necessary data structures for the underlying matrix is allocated.

    703:   Input Parameter:
    704: . B - A `MATLMVM` matrix

    706:   Output Parameter:
    707: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise

    709:   Level: intermediate

    711: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
    712: @*/
    713: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
    714: {
    715:   PetscBool same;

    717:   PetscFunctionBegin;
    719:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    720:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    721:   *flg = B->preallocated;
    722:   PetscFunctionReturn(PETSC_SUCCESS);
    723: }

    725: /*@
    726:   MatLMVMAllocate - Produces all necessary common memory for
    727:   LMVM approximations based on the solution and function vectors
    728:   provided.

    730:   Input Parameters:
    731: + B - A `MATLMVM` matrix
    732: . X - Solution vector
    733: - F - Function vector

    735:   Level: intermediate

    737:   Note:
    738:   If `MatSetSizes()` and `MatSetUp()` have not been called
    739:   before `MatLMVMAllocate()`, the row layout of `B` will be set to match `F`
    740:   and the column layout of `B` will be set to match `X`.

    742: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
    743: @*/
    744: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
    745: {
    746:   Mat_LMVM *lmvm;
    747:   PetscBool same;

    749:   PetscFunctionBegin;
    753:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    754:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    755:   lmvm = (Mat_LMVM *)B->data;
    756:   PetscCall(MatAllocate_LMVM(B, X, F));
    757:   PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
    758:   PetscFunctionReturn(PETSC_SUCCESS);
    759: }

    761: /*@
    762:   MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.

    764:   Input Parameter:
    765: . B - A `MATLMVM` matrix

    767:   Level: intermediate

    769: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
    770: @*/
    771: PetscErrorCode MatLMVMResetShift(Mat B)
    772: {
    773:   Mat_LMVM *lmvm;
    774:   PetscBool same;

    776:   PetscFunctionBegin;
    778:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    779:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    780:   lmvm        = (Mat_LMVM *)B->data;
    781:   lmvm->shift = 0.0;
    782:   PetscFunctionReturn(PETSC_SUCCESS);
    783: }

    785: PETSC_INTERN PetscErrorCode MatLMVMReset_Internal(Mat B, MatLMVMResetMode mode)
    786: {
    787:   Mat_LMVM *lmvm;
    788:   PetscBool same;

    790:   PetscFunctionBegin;
    792:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    793:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    794:   lmvm = (Mat_LMVM *)B->data;
    795:   PetscCall(MatLMVMReset_Internal(lmvm->J0, mode));
    796:   if (lmvm->ops->reset) PetscCall((*lmvm->ops->reset)(B, mode));
    797:   PetscCall(MatReset_LMVM(B, mode));
    798:   PetscFunctionReturn(PETSC_SUCCESS);
    799: }

    801: /*@
    802:   MatLMVMReset - Flushes all of the accumulated updates out of
    803:   the `MATLMVM` approximation.

    805:   Input Parameters:
    806: + B           - A `MATLMVM` matrix
    807: - destructive - flag for enabling destruction of data structures

    809:   Level: intermediate

    811:   Note:
    812:   In practice, this will not actually
    813:   destroy the data associated with the updates. It simply resets
    814:   counters, which leads to existing data being overwritten, and
    815:   `MatSolve()` being applied as if there are no updates. A boolean
    816:   flag is available to force destruction of the update vectors.

    818:   If the user has provided another LMVM matrix as J0, the J0
    819:   matrix is also reset to the identity matrix in this function.

    821: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
    822: @*/
    823: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
    824: {
    825:   Mat_LMVM *lmvm;
    826:   PetscBool same;

    828:   PetscFunctionBegin;
    830:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    831:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    832:   lmvm = (Mat_LMVM *)B->data;
    833:   PetscCall(PetscInfo(B, "Resetting %s after %" PetscInt_FMT " iterations\n", ((PetscObject)B)->type_name, lmvm->k));
    834:   PetscCall(MatLMVMReset_Internal(B, destructive ? MAT_LMVM_RESET_ALL : MAT_LMVM_RESET_HISTORY));
    835:   ++lmvm->nresets;
    836:   PetscFunctionReturn(PETSC_SUCCESS);
    837: }

    839: /*@
    840:   MatLMVMSetHistorySize - Set the number of past iterates to be
    841:   stored for the construction of the limited-memory quasi-Newton update.

    843:   Input Parameters:
    844: + B         - A `MATLMVM` matrix
    845: - hist_size - number of past iterates (default 5)

    847:   Options Database Key:
    848: . -mat_lmvm_hist_size <m> - set number of past iterates

    850:   Level: beginner

    852: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
    853: @*/
    854: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
    855: {
    856:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
    857:   PetscBool same;

    859:   PetscFunctionBegin;
    861:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    862:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    863:   PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
    864:   if (lmvm->m != hist_size) PetscCall(MatLMVMReset_Internal(B, MAT_LMVM_RESET_BASES));
    865:   lmvm->m = hist_size;
    866:   PetscFunctionReturn(PETSC_SUCCESS);
    867: }

    869: PetscErrorCode MatLMVMGetHistorySize(Mat B, PetscInt *hist_size)
    870: {
    871:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
    872:   PetscBool same;

    874:   PetscFunctionBegin;
    876:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    877:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    878:   *hist_size = lmvm->m;
    879:   PetscFunctionReturn(PETSC_SUCCESS);
    880: }

    882: /*@
    883:   MatLMVMGetUpdateCount - Returns the number of accepted updates.

    885:   Input Parameter:
    886: . B - A `MATLMVM` matrix

    888:   Output Parameter:
    889: . nupdates - number of accepted updates

    891:   Level: intermediate

    893:   Note:
    894:   This number may be greater than the total number of update vectors
    895:   stored in the matrix (`MatLMVMGetHistorySize()`). The counters are reset when `MatLMVMReset()`
    896:   is called.

    898: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
    899: @*/
    900: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
    901: {
    902:   Mat_LMVM *lmvm;
    903:   PetscBool same;

    905:   PetscFunctionBegin;
    907:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    908:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    909:   lmvm      = (Mat_LMVM *)B->data;
    910:   *nupdates = lmvm->nupdates;
    911:   PetscFunctionReturn(PETSC_SUCCESS);
    912: }

    914: /*@
    915:   MatLMVMGetRejectCount - Returns the number of rejected updates.
    916:   The counters are reset when `MatLMVMReset()` is called.

    918:   Input Parameter:
    919: . B - A `MATLMVM` matrix

    921:   Output Parameter:
    922: . nrejects - number of rejected updates

    924:   Level: intermediate

    926: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
    927: @*/
    928: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
    929: {
    930:   Mat_LMVM *lmvm;
    931:   PetscBool same;

    933:   PetscFunctionBegin;
    935:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
    936:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
    937:   lmvm      = (Mat_LMVM *)B->data;
    938:   *nrejects = lmvm->nrejects;
    939:   PetscFunctionReturn(PETSC_SUCCESS);
    940: }

    942: PETSC_INTERN PetscErrorCode MatLMVMGetJ0Scalar(Mat B, PetscBool *is_scalar, PetscScalar *scale)
    943: {
    944:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    946:   PetscFunctionBegin;
    947:   PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATCONSTANTDIAGONAL, is_scalar));
    948:   if (*is_scalar) PetscCall(MatConstantDiagonalGetConstant(lmvm->J0, scale));
    949:   PetscFunctionReturn(PETSC_SUCCESS);
    950: }

    952: static PetscErrorCode MatLMVMUpdateOpVecs(Mat B, LMBasis X, LMBasis OpX, PetscErrorCode (*op)(Mat, Vec, Vec))
    953: {
    954:   Mat_LMVM        *lmvm = (Mat_LMVM *)B->data;
    955:   PetscObjectId    J0_id;
    956:   PetscObjectState J0_state;
    957:   PetscInt         oldest, next;

    959:   PetscFunctionBegin;
    960:   PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
    961:   PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
    962:   PetscCall(LMBasisGetRange(X, &oldest, &next));
    963:   if (OpX->operator_id != J0_id || OpX->operator_state != J0_state) {
    964:     // invalidate OpX
    965:     OpX->k              = oldest;
    966:     OpX->operator_id    = J0_id;
    967:     OpX->operator_state = J0_state;
    968:     PetscCall(LMBasisSetCachedProduct(OpX, NULL, NULL));
    969:   }
    970:   OpX->k = PetscMax(OpX->k, oldest);
    971:   for (PetscInt i = OpX->k; i < next; i++) {
    972:     Vec x_i, op_x_i;

    974:     PetscCall(LMBasisGetVecRead(X, i, &x_i));
    975:     PetscCall(LMBasisGetNextVec(OpX, &op_x_i));
    976:     PetscCall(op(B, x_i, op_x_i));
    977:     PetscCall(LMBasisRestoreNextVec(OpX, &op_x_i));
    978:     PetscCall(LMBasisRestoreVecRead(X, i, &x_i));
    979:   }
    980:   PetscAssert(OpX->k == X->k && OpX->operator_id == J0_id && OpX->operator_state == J0_state, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Invalid state for operator-updated LMBasis");
    981:   PetscFunctionReturn(PETSC_SUCCESS);
    982: }

    984: static PetscErrorCode MatLMVMUpdateOpDiffVecs(Mat B, LMBasis Y, PetscScalar alpha, LMBasis OpX, LMBasis YmalphaOpX)
    985: {
    986:   Mat_LMVM        *lmvm = (Mat_LMVM *)B->data;
    987:   PetscInt         start;
    988:   PetscObjectId    J0_id;
    989:   PetscObjectState J0_state;
    990:   PetscInt         oldest, next;

    992:   PetscFunctionBegin;
    993:   PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
    994:   PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
    995:   PetscAssert(Y->m == OpX->m, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Incompatible Y and OpX in MatLMVMUpdateOpDiffVecs()");
    996:   PetscAssert(Y->k == OpX->k, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Stale OpX in MatLMVMUpdateOpDiffVecs()");
    997:   PetscCall(LMBasisGetRange(Y, &oldest, &next));
    998:   if (YmalphaOpX->operator_id != J0_id || YmalphaOpX->operator_state != J0_state) {
    999:     // invalidate OpX
    1000:     YmalphaOpX->k              = oldest;
    1001:     YmalphaOpX->operator_id    = J0_id;
    1002:     YmalphaOpX->operator_state = J0_state;
    1003:     PetscCall(LMBasisSetCachedProduct(YmalphaOpX, NULL, NULL));
    1004:   }
    1005:   YmalphaOpX->k = PetscMax(YmalphaOpX->k, oldest);
    1006:   start         = YmalphaOpX->k;
    1007:   if (next - start == Y->m) { // full matrix AXPY
    1008:     PetscCall(MatCopy(Y->vecs, YmalphaOpX->vecs, SAME_NONZERO_PATTERN));
    1009:     PetscCall(MatAXPY(YmalphaOpX->vecs, -alpha, OpX->vecs, SAME_NONZERO_PATTERN));
    1010:     YmalphaOpX->k = Y->k;
    1011:   } else {
    1012:     for (PetscInt i = start; i < next; i++) {
    1013:       Vec y_i, op_x_i, y_m_op_x_i;

    1015:       PetscCall(LMBasisGetVecRead(Y, i, &y_i));
    1016:       PetscCall(LMBasisGetVecRead(OpX, i, &op_x_i));
    1017:       PetscCall(LMBasisGetNextVec(YmalphaOpX, &y_m_op_x_i));
    1018:       PetscCall(VecAXPBYPCZ(y_m_op_x_i, 1.0, -alpha, 0.0, y_i, op_x_i));
    1019:       PetscCall(LMBasisRestoreNextVec(YmalphaOpX, &y_m_op_x_i));
    1020:       PetscCall(LMBasisRestoreVecRead(OpX, i, &op_x_i));
    1021:       PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
    1022:     }
    1023:   }
    1024:   PetscAssert(YmalphaOpX->k == Y->k && YmalphaOpX->operator_id == J0_id && YmalphaOpX->operator_state == J0_state, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Invalid state for operator-updated LMBasis");
    1025:   PetscFunctionReturn(PETSC_SUCCESS);
    1026: }

    1028: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedBasis(Mat B, MatLMVMBasisType type, LMBasis *basis_p, MatLMVMBasisType *returned_type, PetscScalar *scale)
    1029: {
    1030:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
    1031:   LMBasis     basis;
    1032:   PetscBool   is_scalar;
    1033:   PetscScalar scale_;

    1035:   PetscFunctionBegin;
    1036:   switch (type) {
    1037:   case LMBASIS_S:
    1038:   case LMBASIS_Y:
    1039:     *basis_p = lmvm->basis[type];
    1040:     if (returned_type) *returned_type = type;
    1041:     if (scale) *scale = 1.0;
    1042:     break;
    1043:   case LMBASIS_B0S:
    1044:   case LMBASIS_H0Y:
    1045:     // if B_0 = gamma * I, do not actually compute these bases
    1046:     PetscAssertPointer(returned_type, 4);
    1047:     PetscAssertPointer(scale, 5);
    1048:     PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
    1049:     if (is_scalar) {
    1050:       *basis_p       = lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y];
    1051:       *returned_type = (type == LMBASIS_B0S) ? LMBASIS_S : LMBASIS_Y;
    1052:       *scale         = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
    1053:     } else {
    1054:       LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];

    1056:       *returned_type = type;
    1057:       *scale         = 1.0;
    1058:       if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
    1059:       basis = lmvm->basis[type];
    1060:       PetscCall(MatLMVMUpdateOpVecs(B, orig_basis, basis, (type == LMBASIS_B0S) ? MatLMVMApplyJ0Fwd : MatLMVMApplyJ0Inv));
    1061:       *basis_p = basis;
    1062:     }
    1063:     break;
    1064:   case LMBASIS_S_MINUS_H0Y:
    1065:   case LMBASIS_Y_MINUS_B0S: {
    1066:     MatLMVMBasisType op_basis_t = (type == LMBASIS_S_MINUS_H0Y) ? LMBASIS_H0Y : LMBASIS_B0S;
    1067:     LMBasis          op_basis;

    1069:     if (returned_type) *returned_type = type;
    1070:     if (scale) *scale = 1.0;
    1071:     if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
    1072:     basis = lmvm->basis[type];
    1073:     PetscCall(MatLMVMGetUpdatedBasis(B, op_basis_t, &op_basis, &op_basis_t, &scale_));
    1074:     PetscCall(MatLMVMUpdateOpDiffVecs(B, lmvm->basis[MatLMVMBasisSizeOf(type)], scale_, op_basis, basis));
    1075:     *basis_p = basis;
    1076:   } break;
    1077:   default:
    1078:     PetscUnreachable();
    1079:   }
    1080:   basis = *basis_p;
    1081:   PetscFunctionReturn(PETSC_SUCCESS);
    1082: }

    1084: PETSC_INTERN PetscErrorCode MatLMVMBasisGetVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
    1085: {
    1086:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
    1087:   PetscBool   is_scalar;
    1088:   PetscScalar scale_;

    1090:   PetscFunctionBegin;
    1091:   switch (type) {
    1092:   case LMBASIS_B0S:
    1093:   case LMBASIS_H0Y:
    1094:     // if B_0 = gamma * I, do not actually compute these bases
    1095:     PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
    1096:     if (is_scalar) {
    1097:       *scale = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
    1098:       PetscCall(LMBasisGetVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
    1099:     } else if (lmvm->do_not_cache_J0_products) {
    1100:       Vec     tmp;
    1101:       Vec     w;
    1102:       LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];
    1103:       LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];

    1105:       PetscCall(LMBasisGetVecRead(orig_basis, i, &w));
    1106:       PetscCall(LMBasisGetWorkVec(size_basis, &tmp));
    1107:       if (type == LMBASIS_B0S) {
    1108:         PetscCall(MatLMVMApplyJ0Fwd(B, w, tmp));
    1109:       } else {
    1110:         PetscCall(MatLMVMApplyJ0Inv(B, w, tmp));
    1111:       }
    1112:       PetscCall(LMBasisRestoreVecRead(orig_basis, i, &w));
    1113:       *scale = 1.0;
    1114:       *y     = tmp;
    1115:     } else {
    1116:       LMBasis     basis;
    1117:       PetscScalar dummy;

    1119:       PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &type, &dummy));
    1120:       PetscCall(LMBasisGetVecRead(basis, i, y));
    1121:       *scale = 1.0;
    1122:     }
    1123:     break;
    1124:   default:
    1125:     SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisGetVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y.  Use MatLMVMGetUpdatedBasis() and LMBasisGetVecRead().");
    1126:   }
    1127:   PetscFunctionReturn(PETSC_SUCCESS);
    1128: }

    1130: PETSC_INTERN PetscErrorCode MatLMVMBasisRestoreVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
    1131: {
    1132:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
    1133:   PetscBool   is_scalar;
    1134:   PetscScalar scale_;

    1136:   PetscFunctionBegin;
    1137:   switch (type) {
    1138:   case LMBASIS_B0S:
    1139:   case LMBASIS_H0Y:
    1140:     // if B_0 = gamma * I, do not actually compute these bases
    1141:     PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
    1142:     if (is_scalar) {
    1143:       PetscCall(LMBasisRestoreVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
    1144:     } else if (lmvm->do_not_cache_J0_products) {
    1145:       LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];

    1147:       PetscCall(LMBasisRestoreWorkVec(size_basis, y));
    1148:     } else {
    1149:       PetscCall(LMBasisRestoreVecRead(lmvm->basis[type], i, y));
    1150:     }
    1151:     break;
    1152:   default:
    1153:     SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisRestoreVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y.  Use MatLMVMGetUpdatedBasis() and LMBasisRestoreVecRead().");
    1154:   }
    1155:   PetscFunctionReturn(PETSC_SUCCESS);
    1156: }

    1158: PETSC_INTERN PetscErrorCode MatLMVMGetRange(Mat B, PetscInt *oldest, PetscInt *next)
    1159: {
    1160:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    1162:   PetscFunctionBegin;
    1163:   PetscCall(LMBasisGetRange(lmvm->basis[LMBASIS_S], oldest, next));
    1164:   PetscFunctionReturn(PETSC_SUCCESS);
    1165: }

    1167: PETSC_INTERN PetscErrorCode MatLMVMGetWorkRow(Mat B, Vec *array_p)
    1168: {
    1169:   LMBasis basis;

    1171:   PetscFunctionBegin;
    1172:   PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
    1173:   PetscCall(LMBasisGetWorkRow(basis, array_p));
    1174:   PetscFunctionReturn(PETSC_SUCCESS);
    1175: }

    1177: PETSC_INTERN PetscErrorCode MatLMVMRestoreWorkRow(Mat B, Vec *array_p)
    1178: {
    1179:   LMBasis basis;

    1181:   PetscFunctionBegin;
    1182:   PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
    1183:   PetscCall(LMBasisRestoreWorkRow(basis, array_p));
    1184:   PetscFunctionReturn(PETSC_SUCCESS);
    1185: }

    1187: static PetscErrorCode MatLMVMApplyOpThenVecs(PetscScalar alpha, Mat B, PetscInt oldest, PetscInt next, MatLMVMBasisType type_S, PetscErrorCode (*op)(Mat, Vec, Vec), Vec x, PetscScalar beta, Vec y)
    1188: {
    1189:   LMBasis S;
    1190:   Vec     B0x;

    1192:   PetscFunctionBegin;
    1193:   PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
    1194:   PetscCall(LMBasisGetWorkVec(S, &B0x));
    1195:   PetscCall(op(B, x, B0x));
    1196:   PetscCall(LMBasisGEMVH(S, oldest, next, alpha, B0x, beta, y));
    1197:   PetscCall(LMBasisRestoreWorkVec(S, &B0x));
    1198:   PetscFunctionReturn(PETSC_SUCCESS);
    1199: }

    1201: static PetscErrorCode MatLMVMApplyVecsThenOp(PetscScalar alpha, Mat B, PetscInt oldest, PetscInt next, MatLMVMBasisType type_S, MatLMVMBasisType type_Y, PetscErrorCode (*op)(Mat, Vec, Vec), Vec x, PetscScalar beta, Vec y)
    1202: {
    1203:   LMBasis S, Y;
    1204:   Vec     S_x;
    1205:   Vec     B0S_x;

    1207:   PetscFunctionBegin;
    1208:   PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
    1209:   PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, NULL, NULL));
    1210:   PetscCall(LMBasisGetWorkVec(S, &S_x));
    1211:   PetscCall(LMBasisGEMV(S, oldest, next, alpha, x, 0.0, S_x));
    1212:   PetscCall(LMBasisGetWorkVec(Y, &B0S_x));
    1213:   PetscCall(op(B, S_x, B0S_x));
    1214:   PetscCall(VecAYPX(y, beta, B0S_x));
    1215:   PetscCall(LMBasisRestoreWorkVec(Y, &B0S_x));
    1216:   PetscCall(LMBasisRestoreWorkVec(S, &S_x));
    1217:   PetscFunctionReturn(PETSC_SUCCESS);
    1218: }

    1220: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMVH(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec v, PetscScalar beta, Vec array)
    1221: {
    1222:   Mat_LMVM        *lmvm              = (Mat_LMVM *)B->data;
    1223:   PetscBool        cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
    1224:   LMBasis          basis;
    1225:   MatLMVMBasisType basis_t;
    1226:   PetscScalar      gamma;

    1228:   PetscFunctionBegin;
    1229:   if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
    1230:     PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &basis_t, &gamma));
    1231:     PetscCall(LMBasisGEMVH(basis, oldest, next, alpha * gamma, v, beta, array));
    1232:   } else {
    1233:     switch (type) {
    1234:     case LMBASIS_B0S:
    1235:       PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_S, MatLMVMApplyJ0HermitianTranspose, v, beta, array));
    1236:       break;
    1237:     case LMBASIS_H0Y:
    1238:       PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_Y, MatLMVMApplyJ0InvHermitianTranspose, v, beta, array));
    1239:       break;
    1240:     case LMBASIS_Y_MINUS_B0S:
    1241:       PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_Y], oldest, next, alpha, v, beta, array));
    1242:       PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_B0S, oldest, next, -alpha, v, 1.0, array));
    1243:       break;
    1244:     case LMBASIS_S_MINUS_H0Y:
    1245:       PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_S], oldest, next, alpha, v, beta, array));
    1246:       PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_H0Y, oldest, next, -alpha, v, 1.0, array));
    1247:       break;
    1248:     default:
    1249:       PetscUnreachable();
    1250:     }
    1251:   }
    1252:   PetscFunctionReturn(PETSC_SUCCESS);
    1253: }

    1255: // x must come from MatLMVMGetRowWork()
    1256: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMV(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec x, PetscScalar beta, Vec y)
    1257: {
    1258:   Mat_LMVM *lmvm              = (Mat_LMVM *)B->data;
    1259:   PetscBool cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
    1260:   LMBasis   basis;

    1262:   PetscFunctionBegin;
    1263:   if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
    1264:     PetscScalar      gamma;
    1265:     MatLMVMBasisType base_type;

    1267:     PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &base_type, &gamma));
    1268:     PetscCall(LMBasisGEMV(basis, oldest, next, alpha * gamma, x, beta, y));
    1269:   } else {
    1270:     switch (type) {
    1271:     case LMBASIS_B0S:
    1272:       PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_S, LMBASIS_Y, MatLMVMApplyJ0Fwd, x, beta, y));
    1273:       break;
    1274:     case LMBASIS_H0Y:
    1275:       PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_Y, LMBASIS_S, MatLMVMApplyJ0Inv, x, beta, y));
    1276:       break;
    1277:     case LMBASIS_Y_MINUS_B0S:
    1278:       PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_Y], oldest, next, alpha, x, beta, y));
    1279:       PetscCall(MatLMVMBasisGEMV(B, LMBASIS_B0S, oldest, next, -alpha, x, 1.0, y));
    1280:       break;
    1281:     case LMBASIS_S_MINUS_H0Y:
    1282:       PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_S], oldest, next, alpha, x, beta, y));
    1283:       PetscCall(MatLMVMBasisGEMV(B, LMBASIS_H0Y, oldest, next, -alpha, x, 1.0, y));
    1284:       break;
    1285:     default:
    1286:       PetscUnreachable();
    1287:     }
    1288:   }
    1289:   PetscFunctionReturn(PETSC_SUCCESS);
    1290: }

    1292: PETSC_INTERN PetscErrorCode MatLMVMCreateProducts(Mat B, LMBlockType block_type, LMProducts *products)
    1293: {
    1294:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    1296:   PetscFunctionBegin;
    1297:   PetscCall(LMProductsCreate(lmvm->basis[LMBASIS_S], block_type, products));
    1298:   (*products)->debug = lmvm->debug;
    1299:   PetscFunctionReturn(PETSC_SUCCESS);
    1300: }

    1302: static PetscErrorCode MatLMVMProductsUpdate(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type)
    1303: {
    1304:   Mat_LMVM        *lmvm = (Mat_LMVM *)B->data;
    1305:   LMBasis          X, Y;
    1306:   MatLMVMBasisType true_type_X, true_type_Y;
    1307:   PetscScalar      alpha_X, alpha_Y;
    1308:   PetscInt         oldest, next;
    1309:   LMProducts       G;

    1311:   PetscFunctionBegin;
    1312:   PetscCall(MatLMVMGetUpdatedBasis(B, type_X, &X, &true_type_X, &alpha_X));
    1313:   PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, &true_type_Y, &alpha_Y));
    1314:   if (!lmvm->products[block_type][true_type_X][true_type_Y]) PetscCall(MatLMVMCreateProducts(B, block_type, &lmvm->products[block_type][true_type_X][true_type_Y]));
    1315:   PetscCall(LMProductsUpdate(lmvm->products[block_type][true_type_X][true_type_Y], X, Y));
    1316:   if (true_type_X == type_X && true_type_Y == type_Y) PetscFunctionReturn(PETSC_SUCCESS);
    1317:   if (!lmvm->products[block_type][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, block_type, &lmvm->products[block_type][type_X][type_Y]));
    1318:   G = lmvm->products[block_type][type_X][type_Y];
    1319:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    1320:   PetscCall(LMProductsPrepare(G, lmvm->J0, oldest, next));
    1321:   if (G->k < lmvm->k) {
    1322:     PetscCall(LMProductsCopy(lmvm->products[block_type][true_type_X][true_type_Y], lmvm->products[block_type][type_X][type_Y]));
    1323:     if (alpha_X * alpha_Y != 1.0) PetscCall(LMProductsScale(lmvm->products[block_type][type_X][type_Y], alpha_X * alpha_Y));
    1324:   }
    1325:   PetscFunctionReturn(PETSC_SUCCESS);
    1326: }

    1328: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedProducts(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type, LMProducts *lmwd)
    1329: {
    1330:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    1332:   PetscFunctionBegin;
    1333:   PetscCall(MatLMVMProductsUpdate(B, type_X, type_Y, block_type));
    1334:   *lmwd = lmvm->products[block_type][type_X][type_Y];
    1335:   PetscFunctionReturn(PETSC_SUCCESS);
    1336: }

    1338: PETSC_INTERN PetscErrorCode MatLMVMProductsInsertDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar v)
    1339: {
    1340:   Mat_LMVM  *lmvm = (Mat_LMVM *)B->data;
    1341:   LMProducts products;

    1343:   PetscFunctionBegin;
    1344:   if (!lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]));
    1345:   products = lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y];
    1346:   PetscCall(LMProductsInsertNextDiagonalValue(products, idx, v));
    1347:   PetscFunctionReturn(PETSC_SUCCESS);
    1348: }

    1350: PETSC_INTERN PetscErrorCode MatLMVMProductsGetDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar *v)
    1351: {
    1352:   LMProducts products = NULL;

    1354:   PetscFunctionBegin;
    1355:   PetscCall(MatLMVMGetUpdatedProducts(B, type_X, type_Y, LMBLOCK_DIAGONAL, &products));
    1356:   PetscCall(LMProductsGetDiagonalValue(products, idx, v));
    1357:   PetscFunctionReturn(PETSC_SUCCESS);
    1358: }