Actual source code: brdn.c

      1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>
      2: #include <petscblaslapack.h>

      4: // Broyden is dual to bad Broyden: MatSolve routines are dual to bad Broyden MatMult routines

      6: static PetscErrorCode MatSolve_LMVMBrdn_Recursive(Mat B, Vec F, Vec dX)
      7: {
      8:   PetscFunctionBegin;
      9:   PetscCall(BadBroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
     10:   PetscFunctionReturn(PETSC_SUCCESS);
     11: }

     13: static PetscErrorCode MatSolveHermitianTranspose_LMVMBrdn_Recursive(Mat B, Vec F, Vec dX)
     14: {
     15:   PetscFunctionBegin;
     16:   PetscCall(BadBroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: static PetscErrorCode MatSolve_LMVMBrdn_CompactDense(Mat B, Vec F, Vec dX)
     21: {
     22:   PetscFunctionBegin;
     23:   PetscCall(BadBroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: static PetscErrorCode MatSolveHermitianTranspose_LMVMBrdn_CompactDense(Mat B, Vec F, Vec dX)
     28: {
     29:   PetscFunctionBegin;
     30:   PetscCall(BadBroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
     31:   PetscFunctionReturn(PETSC_SUCCESS);
     32: }

     34: /* The Broyden kernel can be written as

     36:      B_k = B_{k-1} + (y_{k-1} - B_{k-1} s_{k-1}) (s_{k-1}^H s_{k-1})^{-1} s_{k-1}^H

     38:    This means that we can write B_k x using an intermediate variable alpha_k as

     40:      alpha_k = (s_{k-1}^H s_{k-1})^{-1} s_{k-1}^H x
     41:      B_k x   = y_{k-1} alpha_k + B_{k-1}(x - s_{k-1} alpha_k)
     42:                                  \_____________ ____________/
     43:                                                V
     44:                                      recursive application
     45:  */

     47: PETSC_INTERN PetscErrorCode BroydenKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BX)
     48: {
     49:   Vec              G   = X;
     50:   Vec              W   = BX;
     51:   MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
     52:   MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
     53:   LMBasis          S = NULL, Y = NULL;
     54:   PetscInt         oldest, next;

     56:   PetscFunctionBegin;
     57:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
     58:   if (next > oldest) {
     59:     LMProducts StS;

     61:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
     62:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
     63:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_DIAGONAL, &StS));

     65:     PetscCall(LMBasisGetWorkVec(S, &G));
     66:     PetscCall(VecCopy(X, G));
     67:     PetscCall(VecZeroEntries(BX));

     69:     for (PetscInt i = next - 1; i >= oldest; i--) {
     70:       Vec         s_i, y_i;
     71:       PetscScalar sitg, sitsi, alphai;

     73:       PetscCall(LMBasisGetVecRead(S, i, &s_i));
     74:       PetscCall(LMBasisGetVecRead(Y, i, &y_i));
     75:       PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));

     77:       PetscCall(VecDot(G, s_i, &sitg));
     78:       alphai = sitg / sitsi;
     79:       PetscCall(VecAXPY(BX, alphai, y_i));
     80:       PetscCall(VecAXPY(G, -alphai, s_i));

     82:       PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
     83:       PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
     84:     }

     86:     PetscCall(LMBasisGetWorkVec(Y, &W));
     87:   }
     88:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, W));
     89:   if (next > oldest) {
     90:     PetscCall(VecAXPY(BX, 1.0, W));
     91:     PetscCall(LMBasisRestoreWorkVec(Y, &W));
     92:     PetscCall(LMBasisRestoreWorkVec(S, &G));
     93:   }
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: /* The adjoint of the kernel is

     99:      B_k^H = B_{k-1}^H + s_{k-1} (s_{k-1}^H s_{k-1})^{-1} (y_{k-1} - B_{k-1} s_{k-1})^H

    101:    This means that we can write B_k^H x using an intermediate variable w_k = B_{k-1}^H x

    103:      w_k = B_{k-1}^H x <-- recursive application
    104:      B_k^H x = w_k + s_{k-1} (s_{k-1}^H s_{k-1})^{-1} (y_{k-1}^H x - s_k^H w_k)
    105:  */

    107: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
    108: {
    109:   PetscInt oldest, next;

    111:   PetscFunctionBegin;
    112:   PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
    113:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    114:   if (next > oldest) {
    115:     MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
    116:     MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
    117:     LMBasis          S, Y;
    118:     LMProducts       StS;

    120:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_DIAGONAL, &StS));
    121:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    122:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));

    124:     for (PetscInt i = oldest; i < next; i++) {
    125:       Vec         s_i, y_i;
    126:       PetscScalar sitBHX, sitsi, yitx;

    128:       PetscCall(LMBasisGetVecRead(S, i, &s_i));
    129:       PetscCall(LMBasisGetVecRead(Y, i, &y_i));
    130:       PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));

    132:       PetscCall(VecDotBegin(BHX, s_i, &sitBHX));
    133:       PetscCall(VecDotBegin(X, y_i, &yitx));
    134:       PetscCall(VecDotEnd(BHX, s_i, &sitBHX));
    135:       PetscCall(VecDotEnd(X, y_i, &yitx));
    136:       PetscCall(VecAXPY(BHX, (yitx - sitBHX) / sitsi, s_i));
    137:       PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
    138:       PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
    139:     }
    140:   }
    141:   PetscFunctionReturn(PETSC_SUCCESS);
    142: }

    144: /*
    145:    The Broyden kernel can be written as

    147:      B_k = B_0 + (Y - B_0 S) (triu(S^H S))^{-1} S^H

    149:    where triu is the upper triangular component.  We solve by back substitution each time
    150:    we apply
    151:  */

    153: PETSC_INTERN PetscErrorCode BroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
    154: {
    155:   PetscInt oldest, next;

    157:   PetscFunctionBegin;
    158:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX));
    159:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    160:   if (next > oldest) {
    161:     MatLMVMBasisType S_t           = LMVMModeMap(LMBASIS_S, mode);
    162:     MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
    163:     LMBasis          S;
    164:     LMProducts       StS;
    165:     Vec              StX;

    167:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    168:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
    169:     PetscCall(MatLMVMGetWorkRow(B, &StX));
    170:     PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, X, 0.0, StX));
    171:     PetscCall(LMProductsSolve(StS, oldest, next, StX, StX, /* ^H */ PETSC_FALSE));
    172:     PetscCall(MatLMVMBasisGEMV(B, Y_minus_B0S_t, oldest, next, 1.0, StX, 1.0, BX));
    173:     PetscCall(MatLMVMRestoreWorkRow(B, &StX));
    174:   }
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: /*
    179:    The adjoint of the above formula is

    181:      B_k^H = B_0^H + S^H (triu(S^H S))^{-H} (Y - B_0 S)^H
    182:  */

    184: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
    185: {
    186:   PetscInt oldest, next;

    188:   PetscFunctionBegin;
    189:   PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
    190:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    191:   if (next > oldest) {
    192:     MatLMVMBasisType S_t           = LMVMModeMap(LMBASIS_S, mode);
    193:     MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
    194:     LMProducts       StS;
    195:     LMBasis          S;
    196:     Vec              YmB0StX;

    198:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    199:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
    200:     PetscCall(MatLMVMGetWorkRow(B, &YmB0StX));
    201:     PetscCall(MatLMVMBasisGEMVH(B, Y_minus_B0S_t, oldest, next, 1.0, X, 0.0, YmB0StX));
    202:     PetscCall(LMProductsSolve(StS, oldest, next, YmB0StX, YmB0StX, /* ^H */ PETSC_TRUE));
    203:     PetscCall(LMBasisGEMV(S, oldest, next, 1.0, YmB0StX, 1.0, BHX));
    204:     PetscCall(MatLMVMRestoreWorkRow(B, &YmB0StX));
    205:   }
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: /*
    210:    The dense Broyden formula use for CompactDense can be written in a non-update form as

    212:      B_k = [B_0 | Y] [ I | -S ] [          I           ]
    213:                      [---+----] [----------------------]
    214:                      [ 0 |  I ] [ triu(S^H S)^{-1} S^H ]

    216:    The advantage of this form is B_0 appears only once and the (Y - B_0 S) vectors are not needed
    217:  */

    219: PETSC_INTERN PetscErrorCode BroydenKernel_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
    220: {
    221:   Vec              G   = X;
    222:   Vec              W   = BX;
    223:   MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
    224:   MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
    225:   LMBasis          S = NULL, Y = NULL;
    226:   PetscInt         oldest, next;

    228:   PetscFunctionBegin;
    229:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    230:   if (next > oldest) {
    231:     LMProducts StS;
    232:     Vec        StX;

    234:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    235:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
    236:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
    237:     PetscCall(LMBasisGetWorkVec(S, &G));
    238:     PetscCall(MatLMVMGetWorkRow(B, &StX));
    239:     PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, X, 0.0, StX));
    240:     PetscCall(LMProductsSolve(StS, oldest, next, StX, StX, PETSC_FALSE));
    241:     PetscCall(VecCopy(X, G));
    242:     PetscCall(LMBasisGEMV(S, oldest, next, -1.0, StX, 1.0, G));
    243:     PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, StX, 0.0, BX));
    244:     PetscCall(MatLMVMRestoreWorkRow(B, &StX));
    245:     PetscCall(LMBasisGetWorkVec(Y, &W));
    246:   }
    247:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, W));
    248:   if (next > oldest) {
    249:     PetscCall(VecAXPY(BX, 1.0, W));
    250:     PetscCall(LMBasisRestoreWorkVec(Y, &W));
    251:     PetscCall(LMBasisRestoreWorkVec(S, &G));
    252:   }
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }

    256: /*
    257:    The adoint of the above formula is

    259:      B_k^h = [I | S triu(S^H S)^{-H} ] [  I   | 0 ] [ B_0^H ]
    260:                                        [------+---] [-------]
    261:                                        [ -S^H | I ] [  Y^H  ]
    262:  */

    264: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
    265: {
    266:   PetscInt oldest, next;

    268:   PetscFunctionBegin;
    269:   PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
    270:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    271:   if (next > oldest) {
    272:     MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
    273:     MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
    274:     LMBasis          S, Y;
    275:     LMProducts       StS;
    276:     Vec              YtX, StBHX;

    278:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
    279:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    280:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
    281:     PetscCall(MatLMVMGetWorkRow(B, &YtX));
    282:     PetscCall(MatLMVMGetWorkRow(B, &StBHX));
    283:     PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, BHX, 0.0, StBHX));
    284:     PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
    285:     PetscCall(VecAXPY(YtX, -1.0, StBHX));
    286:     PetscCall(LMProductsSolve(StS, oldest, next, YtX, YtX, PETSC_TRUE));
    287:     PetscCall(LMBasisGEMV(S, oldest, next, 1.0, YtX, 1.0, BHX));
    288:     PetscCall(MatLMVMRestoreWorkRow(B, &StBHX));
    289:     PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
    290:   }
    291:   PetscFunctionReturn(PETSC_SUCCESS);
    292: }

    294: static PetscErrorCode MatMult_LMVMBrdn_Recursive(Mat B, Vec X, Vec Z)
    295: {
    296:   PetscFunctionBegin;
    297:   PetscCall(BroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Z));
    298:   PetscFunctionReturn(PETSC_SUCCESS);
    299: }

    301: static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_Recursive(Mat B, Vec X, Vec Z)
    302: {
    303:   PetscFunctionBegin;
    304:   PetscCall(BroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_PRIMAL, X, Z));
    305:   PetscFunctionReturn(PETSC_SUCCESS);
    306: }

    308: static PetscErrorCode MatMult_LMVMBrdn_CompactDense(Mat B, Vec X, Vec Z)
    309: {
    310:   PetscFunctionBegin;
    311:   PetscCall(BroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Z));
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }

    315: static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_CompactDense(Mat B, Vec X, Vec Z)
    316: {
    317:   PetscFunctionBegin;
    318:   PetscCall(BroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Z));
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: PETSC_UNUSED static PetscErrorCode MatMult_LMVMBrdn_Dense(Mat B, Vec X, Vec Z)
    323: {
    324:   PetscFunctionBegin;
    325:   PetscCall(BroydenKernel_Dense(B, MATLMVM_MODE_PRIMAL, X, Z));
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: PETSC_UNUSED static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_Dense(Mat B, Vec X, Vec Z)
    330: {
    331:   PetscFunctionBegin;
    332:   PetscCall(BroydenKernelHermitianTranspose_Dense(B, MATLMVM_MODE_PRIMAL, X, Z));
    333:   PetscFunctionReturn(PETSC_SUCCESS);
    334: }

    336: static PetscErrorCode MatUpdate_LMVMBrdn(Mat B, Vec X, Vec F)
    337: {
    338:   Mat_LMVM *lmvm          = (Mat_LMVM *)B->data;
    339:   Mat_Brdn *brdn          = (Mat_Brdn *)lmvm->ctx;
    340:   PetscBool cache_YtFprev = (lmvm->mult_alg != MAT_LMVM_MULT_RECURSIVE) ? lmvm->cache_gradient_products : PETSC_FALSE;

    342:   PetscFunctionBegin;
    343:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
    344:   if (lmvm->prev_set) {
    345:     LMBasis    Y;
    346:     PetscInt   oldest, next;
    347:     Vec        Fprev_old   = NULL;
    348:     Vec        YtFprev_old = NULL;
    349:     LMProducts YtY         = NULL;

    351:     PetscCall(MatLMVMGetRange(B, &oldest, &next));
    352:     PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));

    354:     if (cache_YtFprev) {
    355:       if (!brdn->YtFprev) PetscCall(LMBasisCreateRow(Y, &brdn->YtFprev));
    356:       PetscCall(LMBasisGetWorkVec(Y, &Fprev_old));
    357:       PetscCall(MatLMVMGetUpdatedProducts(B, LMBASIS_Y, LMBASIS_Y, LMBLOCK_UPPER_TRIANGLE, &YtY));
    358:       PetscCall(LMProductsGetNextColumn(YtY, &YtFprev_old));
    359:       PetscCall(VecCopy(lmvm->Fprev, Fprev_old));
    360:       PetscCall(VecCopy(brdn->YtFprev, YtFprev_old));
    361:     }
    362:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
    363:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
    364:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
    365:     /* Accept the update */
    366:     PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
    367:     if (cache_YtFprev) {
    368:       PetscInt oldest_new, next_new;

    370:       PetscCall(MatLMVMGetRange(B, &oldest_new, &next_new));
    371:       // Compute the one new Y_i^T Fprev_old value
    372:       PetscCall(LMBasisGEMVH(Y, next, next_new, 1.0, Fprev_old, 0.0, YtFprev_old));
    373:       PetscCall(LMBasisGEMVH(Y, oldest_new, next_new, 1.0, F, 0.0, brdn->YtFprev));
    374:       PetscCall(LMBasisSetCachedProduct(Y, F, brdn->YtFprev));
    375:       PetscCall(VecAXPBY(YtFprev_old, 1.0, -1.0, brdn->YtFprev));
    376:       PetscCall(LMProductsRestoreNextColumn(YtY, &YtFprev_old));
    377:       PetscCall(LMBasisRestoreWorkVec(Y, &Fprev_old));
    378:     }
    379:   } else if (cache_YtFprev) {
    380:     LMBasis Y;

    382:     PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
    383:     if (!brdn->YtFprev) PetscCall(LMBasisCreateRow(Y, &brdn->YtFprev));
    384:   }
    385:   /* Save the solution and function to be used in the next update */
    386:   PetscCall(VecCopy(X, lmvm->Xprev));
    387:   PetscCall(VecCopy(F, lmvm->Fprev));
    388:   lmvm->prev_set = PETSC_TRUE;
    389:   PetscFunctionReturn(PETSC_SUCCESS);
    390: }

    392: static PetscErrorCode MatReset_LMVMBrdn(Mat B, MatLMVMResetMode mode)
    393: {
    394:   Mat_LMVM *lmvm  = (Mat_LMVM *)B->data;
    395:   Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;

    397:   PetscFunctionBegin;
    398:   if (MatLMVMResetClearsBases(mode)) {
    399:     for (PetscInt i = 0; i < BROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisDestroy(&lbrdn->basis[i]));
    400:     for (PetscInt i = 0; i < BROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsDestroy(&lbrdn->products[i]));
    401:     PetscCall(VecDestroy(&lbrdn->YtFprev));
    402:   } else {
    403:     for (PetscInt i = 0; i < BROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisReset(lbrdn->basis[i]));
    404:     for (PetscInt i = 0; i < BROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsReset(lbrdn->products[i]));
    405:     if (lbrdn->YtFprev) PetscCall(VecZeroEntries(lbrdn->YtFprev));
    406:   }
    407:   PetscFunctionReturn(PETSC_SUCCESS);
    408: }

    410: static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
    411: {
    412:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    414:   PetscFunctionBegin;
    415:   PetscCall(MatReset_LMVMBrdn(B, MAT_LMVM_RESET_ALL));
    416:   PetscCall(PetscFree(lmvm->ctx));
    417:   PetscCall(MatDestroy_LMVM(B));
    418:   PetscFunctionReturn(PETSC_SUCCESS);
    419: }

    421: static PetscErrorCode MatLMVMSetMultAlgorithm_Brdn(Mat B)
    422: {
    423:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    425:   PetscFunctionBegin;
    426:   switch (lmvm->mult_alg) {
    427:   case MAT_LMVM_MULT_RECURSIVE:
    428:     lmvm->ops->mult    = MatMult_LMVMBrdn_Recursive;
    429:     lmvm->ops->multht  = MatMultHermitianTranspose_LMVMBrdn_Recursive;
    430:     lmvm->ops->solve   = MatSolve_LMVMBrdn_Recursive;
    431:     lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_Recursive;
    432:     break;
    433:   case MAT_LMVM_MULT_DENSE:
    434:     lmvm->ops->mult    = MatMult_LMVMBrdn_Dense;
    435:     lmvm->ops->multht  = MatMultHermitianTranspose_LMVMBrdn_Dense;
    436:     lmvm->ops->solve   = MatSolve_LMVMBrdn_CompactDense;
    437:     lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_CompactDense;
    438:     break;
    439:   case MAT_LMVM_MULT_COMPACT_DENSE:
    440:     lmvm->ops->mult    = MatMult_LMVMBrdn_CompactDense;
    441:     lmvm->ops->multht  = MatMultHermitianTranspose_LMVMBrdn_CompactDense;
    442:     lmvm->ops->solve   = MatSolve_LMVMBrdn_CompactDense;
    443:     lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_CompactDense;
    444:     break;
    445:   }
    446:   PetscFunctionReturn(PETSC_SUCCESS);
    447: }

    449: PetscErrorCode MatCreate_LMVMBrdn(Mat B)
    450: {
    451:   Mat_LMVM *lmvm;
    452:   Mat_Brdn *lbrdn;

    454:   PetscFunctionBegin;
    455:   PetscCall(MatCreate_LMVM(B));
    456:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBROYDEN));
    457:   B->ops->destroy = MatDestroy_LMVMBrdn;

    459:   lmvm                        = (Mat_LMVM *)B->data;
    460:   lmvm->ops->reset            = MatReset_LMVMBrdn;
    461:   lmvm->ops->update           = MatUpdate_LMVMBrdn;
    462:   lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_Brdn;

    464:   PetscCall(PetscNew(&lbrdn));
    465:   lmvm->ctx = (void *)lbrdn;

    467:   PetscCall(MatLMVMSetMultAlgorithm_Brdn(B));
    468:   PetscFunctionReturn(PETSC_SUCCESS);
    469: }

    471: /*@
    472:   MatCreateLMVMBroyden - Creates a limited-memory "good" Broyden-type approximation
    473:   matrix used for a Jacobian. L-Brdn is not guaranteed to be symmetric or
    474:   positive-definite.

    476:   To use the L-Brdn matrix with other vector types, the matrix must be
    477:   created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
    478:   This ensures that the internal storage and work vectors are duplicated from the
    479:   correct type of vector.

    481:   Collective

    483:   Input Parameters:
    484: + comm - MPI communicator
    485: . n    - number of local rows for storage vectors
    486: - N    - global size of the storage vectors

    488:   Output Parameter:
    489: . B - the matrix

    491:   Options Database Keys:
    492: + -mat_lmvm_hist_size         - the number of history vectors to keep
    493: . -mat_lmvm_mult_algorithm    - the algorithm to use for multiplication (recursive, dense, compact_dense)
    494: . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed
    495: - -mat_lmvm_debug             - (developer) perform internal debugging checks

    497:   Level: intermediate

    499:   Note:
    500:   It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
    501:   paradigm instead of this routine directly.

    503: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
    504:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBadBroyden()`, `MatCreateLMVMSymBroyden()`
    505: @*/
    506: PetscErrorCode MatCreateLMVMBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    507: {
    508:   PetscFunctionBegin;
    509:   PetscCall(KSPInitializePackage());
    510:   PetscCall(MatCreate(comm, B));
    511:   PetscCall(MatSetSizes(*B, n, n, N, N));
    512:   PetscCall(MatSetType(*B, MATLMVMBROYDEN));
    513:   PetscCall(MatSetUp(*B));
    514:   PetscFunctionReturn(PETSC_SUCCESS);
    515: }