Actual source code: dfp.c

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

      3: /* The DFP update can be written as

      5:      B_{k+1} = (I - y_k (y_k^T s_k)^{-1} s_k^T) B_k (I - s_k (y_k^T s_k)^{-1} y_k^T) + y_k (y_k^T s_k)^{-1} y_k^T

      7:    So B_{k+1}x can be computed in the following way

      9:      a_k = (y_k^T s_k)^{-1} y_k^T x
     10:      g = B_k (x - a s_k) <--- recursion
     11:      B_{k+1}x = g + y_k(a - (y_k^T s_k)^{-1} s_k^T g)
     12:  */
     13: PETSC_INTERN PetscErrorCode DFPKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BX)
     14: {
     15:   MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
     16:   MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
     17:   LMBasis          S, Y;
     18:   LMProducts       YtS;
     19:   Vec              G     = X;
     20:   PetscScalar     *alpha = NULL;
     21:   PetscInt         oldest, next;

     23:   PetscFunctionBegin;
     24:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
     25:   PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
     26:   PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
     27:   PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_DIAGONAL, &YtS));
     28:   if (next > oldest) {
     29:     PetscCall(LMBasisGetWorkVec(S, &G));
     30:     PetscCall(VecCopy(X, G));
     31:     PetscCall(PetscMalloc1(next - oldest, &alpha));
     32:   }
     33:   for (PetscInt i = next - 1; i >= oldest; i--) {
     34:     Vec         s_i, y_i;
     35:     PetscScalar yitsi, yitx, a;

     37:     PetscCall(LMBasisGetVecRead(Y, i, &y_i));
     38:     PetscCall(VecDot(G, y_i, &yitx));
     39:     PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
     40:     PetscCall(LMProductsGetDiagonalValue(YtS, i, &yitsi));
     41:     alpha[i - oldest] = a = yitx / yitsi;
     42:     PetscCall(LMBasisGetVecRead(S, i, &s_i));
     43:     PetscCall(VecAXPY(G, -a, s_i));
     44:     PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
     45:   }
     46:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, BX));
     47:   for (PetscInt i = oldest; i < next; i++) {
     48:     Vec         s_i, y_i;
     49:     PetscScalar yitsi, sitbx, a, b;

     51:     PetscCall(LMBasisGetVecRead(S, i, &s_i));
     52:     PetscCall(VecDot(BX, s_i, &sitbx));
     53:     PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
     54:     PetscCall(LMProductsGetDiagonalValue(YtS, i, &yitsi));
     55:     a = alpha[i - oldest];
     56:     b = sitbx / yitsi;
     57:     PetscCall(LMBasisGetVecRead(Y, i, &y_i));
     58:     PetscCall(VecAXPY(BX, a - b, y_i));
     59:     PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
     60:   }
     61:   if (next > oldest) {
     62:     PetscCall(PetscFree(alpha));
     63:     PetscCall(LMBasisRestoreWorkVec(S, &G));
     64:   }
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: PETSC_INTERN PetscErrorCode DFPKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
     69: {
     70:   PetscInt oldest, next;

     72:   PetscFunctionBegin;
     73:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX));
     74:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
     75:   if (next > oldest) {
     76:     MatLMVMBasisType S_t   = LMVMModeMap(LMBASIS_S, mode);
     77:     MatLMVMBasisType Y_t   = LMVMModeMap(LMBASIS_Y, mode);
     78:     MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
     79:     Vec              StB0X, YtX, u, v;
     80:     PetscBool        use_B0S;
     81:     LMBasis          S, Y;
     82:     LMProducts       YtS, StB0S, D;

     84:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
     85:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
     86:     PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_UPPER_TRIANGLE, &YtS));
     87:     PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_DIAGONAL, &D));
     88:     PetscCall(MatLMVMGetUpdatedProducts(B, S_t, B0S_t, LMBLOCK_UPPER_TRIANGLE, &StB0S));

     90:     PetscCall(MatLMVMGetWorkRow(B, &StB0X));
     91:     PetscCall(MatLMVMGetWorkRow(B, &YtX));
     92:     PetscCall(MatLMVMGetWorkRow(B, &u));
     93:     PetscCall(MatLMVMGetWorkRow(B, &v));

     95:     PetscCall(SymBroydenCompactDenseKernelUseB0S(B, mode, X, &use_B0S));
     96:     if (use_B0S) PetscCall(MatLMVMBasisGEMVH(B, B0S_t, oldest, next, 1.0, X, 0.0, StB0X));
     97:     else PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, BX, 0.0, StB0X));

     99:     PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
    100:     PetscCall(LMProductsSolve(YtS, oldest, next, YtX, YtX, /* ^H */ PETSC_FALSE));

    102:     PetscCall(VecAXPBY(u, -1.0, 0.0, YtX));
    103:     PetscCall(LMProductsMult(D, oldest, next, 1.0, YtX, 0.0, v, /* ^H */ PETSC_FALSE));
    104:     PetscCall(LMProductsMultHermitian(StB0S, oldest, next, 1.0, YtX, 1.0, v));
    105:     PetscCall(VecAXPY(v, -1.0, StB0X));

    107:     PetscCall(LMProductsSolve(YtS, oldest, next, v, v, /* ^H */ PETSC_TRUE));
    108:     PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, v, 1.0, BX));
    109:     PetscCall(MatLMVMBasisGEMV(B, B0S_t, oldest, next, 1.0, u, 1.0, BX));

    111:     PetscCall(MatLMVMRestoreWorkRow(B, &v));
    112:     PetscCall(MatLMVMRestoreWorkRow(B, &u));
    113:     PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
    114:     PetscCall(MatLMVMRestoreWorkRow(B, &StB0X));
    115:   }
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: PETSC_INTERN PetscErrorCode DFPKernel_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
    120: {
    121:   Vec        G   = X;
    122:   Vec        YtX = NULL, u = NULL;
    123:   PetscInt   oldest, next;
    124:   LMProducts YtS = NULL, D = NULL;
    125:   LMBasis    S = NULL, Y = NULL;

    127:   PetscFunctionBegin;
    128:   PetscCall(MatLMVMGetRange(B, &oldest, &next));
    129:   if (next > oldest) {
    130:     MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
    131:     MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);

    133:     PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
    134:     PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
    135:     PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_UPPER_TRIANGLE, &YtS));
    136:     PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_DIAGONAL, &D));
    137:     PetscCall(LMBasisGetWorkVec(Y, &G));
    138:     PetscCall(VecCopy(X, G));

    140:     PetscCall(MatLMVMGetWorkRow(B, &YtX));
    141:     PetscCall(MatLMVMGetWorkRow(B, &u));
    142:     PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
    143:     PetscCall(LMProductsSolve(YtS, oldest, next, YtX, YtX, /* ^H */ PETSC_FALSE));
    144:     PetscCall(LMBasisGEMV(S, oldest, next, -1.0, YtX, 1.0, G));
    145:   }
    146:   PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, BX));
    147:   if (next > oldest) {
    148:     PetscCall(LMProductsMult(D, oldest, next, 1.0, YtX, 0.0, u, /* ^H */ PETSC_FALSE));
    149:     PetscCall(LMBasisGEMVH(S, oldest, next, -1.0, BX, 1.0, u));
    150:     PetscCall(LMProductsSolve(YtS, oldest, next, u, u, /* ^H */ PETSC_TRUE));
    151:     PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, u, 1.0, BX));
    152:     PetscCall(MatLMVMRestoreWorkRow(B, &u));
    153:     PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
    154:     PetscCall(LMBasisRestoreWorkVec(Y, &G));
    155:   }
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: static PetscErrorCode MatMult_LMVMDFP_Recursive(Mat B, Vec X, Vec BX)
    160: {
    161:   PetscFunctionBegin;
    162:   PetscCall(DFPKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, BX));
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }

    166: static PetscErrorCode MatMult_LMVMDFP_CompactDense(Mat B, Vec X, Vec BX)
    167: {
    168:   PetscFunctionBegin;
    169:   PetscCall(DFPKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, BX));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: static PetscErrorCode MatMult_LMVMDFP_Dense(Mat B, Vec X, Vec BX)
    174: {
    175:   PetscFunctionBegin;
    176:   PetscCall(DFPKernel_Dense(B, MATLMVM_MODE_PRIMAL, X, BX));
    177:   PetscFunctionReturn(PETSC_SUCCESS);
    178: }

    180: static PetscErrorCode MatSolve_LMVMDFP_Recursive(Mat B, Vec X, Vec HX)
    181: {
    182:   PetscFunctionBegin;
    183:   PetscCall(BFGSKernel_Recursive(B, MATLMVM_MODE_DUAL, X, HX));
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: static PetscErrorCode MatSolve_LMVMDFP_CompactDense(Mat B, Vec X, Vec HX)
    188: {
    189:   PetscFunctionBegin;
    190:   PetscCall(BFGSKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, HX));
    191:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }

    194: static PetscErrorCode MatSetFromOptions_LMVMDFP(Mat B, PetscOptionItems PetscOptionsObject)
    195: {
    196:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
    197:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;

    199:   PetscFunctionBegin;
    200:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
    201:   PetscOptionsHeadBegin(PetscOptionsObject, "DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
    202:   PetscCall(SymBroydenRescaleSetFromOptions(B, ldfp->rescale, PetscOptionsObject));
    203:   PetscOptionsHeadEnd();
    204:   PetscFunctionReturn(PETSC_SUCCESS);
    205: }

    207: static PetscErrorCode MatLMVMSetMultAlgorithm_DFP(Mat B)
    208: {
    209:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

    211:   PetscFunctionBegin;
    212:   switch (lmvm->mult_alg) {
    213:   case MAT_LMVM_MULT_RECURSIVE:
    214:     lmvm->ops->mult  = MatMult_LMVMDFP_Recursive;
    215:     lmvm->ops->solve = MatSolve_LMVMDFP_Recursive;
    216:     break;
    217:   case MAT_LMVM_MULT_DENSE:
    218:     lmvm->ops->mult  = MatMult_LMVMDFP_Dense;
    219:     lmvm->ops->solve = MatSolve_LMVMDFP_CompactDense;
    220:     break;
    221:   case MAT_LMVM_MULT_COMPACT_DENSE:
    222:     lmvm->ops->mult  = MatMult_LMVMDFP_CompactDense;
    223:     lmvm->ops->solve = MatSolve_LMVMDFP_CompactDense;
    224:     break;
    225:   }
    226:   lmvm->ops->multht  = lmvm->ops->mult;
    227:   lmvm->ops->solveht = lmvm->ops->solve;
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: PetscErrorCode MatCreate_LMVMDFP(Mat B)
    232: {
    233:   Mat_LMVM    *lmvm;
    234:   Mat_SymBrdn *ldfp;

    236:   PetscFunctionBegin;
    237:   PetscCall(MatCreate_LMVMSymBrdn(B));
    238:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP));
    239:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;

    241:   lmvm                        = (Mat_LMVM *)B->data;
    242:   lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_DFP;
    243:   PetscCall(MatLMVMSetMultAlgorithm_DFP(B));

    245:   ldfp             = (Mat_SymBrdn *)lmvm->ctx;
    246:   ldfp->phi_scalar = 1.0;
    247:   ldfp->psi_scalar = 0.0;
    248:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBadBroydenSetPsi_C", NULL));
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: /*@
    253:   MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
    254:   used for approximating Jacobians. L-DFP is symmetric positive-definite by
    255:   construction, and is the dual of L-BFGS where Y and S vectors swap roles.

    257:   To use the L-DFP matrix with other vector types, the matrix must be
    258:   created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
    259:   This ensures that the internal storage and work vectors are duplicated from the
    260:   correct type of vector.

    262:   Collective

    264:   Input Parameters:
    265: + comm - MPI communicator
    266: . n    - number of local rows for storage vectors
    267: - N    - global size of the storage vectors

    269:   Output Parameter:
    270: . B - the matrix

    272:   Options Database Keys:
    273: + -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
    274: . -mat_lmvm_theta      - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    275: . -mat_lmvm_rho        - (developer) update limiter for the J0 scaling
    276: . -mat_lmvm_alpha      - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    277: . -mat_lmvm_beta       - (developer) exponential factor for the diagonal J0 scaling
    278: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

    280:   Level: intermediate

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

    286: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDFP`, `MatCreateLMVMBFGS()`, `MatCreateLMVMSR1()`,
    287:           `MatCreateLMVMBroyden()`, `MatCreateLMVMBadBroyden()`, `MatCreateLMVMSymBroyden()`
    288: @*/
    289: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    290: {
    291:   PetscFunctionBegin;
    292:   PetscCall(KSPInitializePackage());
    293:   PetscCall(MatCreate(comm, B));
    294:   PetscCall(MatSetSizes(*B, n, n, N, N));
    295:   PetscCall(MatSetType(*B, MATLMVMDFP));
    296:   PetscCall(MatSetUp(*B));
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }