Actual source code: symbadbrdn.c

      1: #include "symbrdn.h" /*I "petscksp.h" I*/

      3: static PetscErrorCode MatMult_LMVMSymBadBrdn_Recursive(Mat B, Vec X, Vec Y)
      4: {
      5:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
      6:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

      8:   PetscFunctionBegin;
      9:   PetscCheck(lsb->psi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Psi must first be set using MatLMVMSymBadBroydenSetPsi()");
     10:   if (lsb->psi_scalar == 0.0) {
     11:     PetscCall(DFPKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y));
     12:   } else if (lsb->psi_scalar == 1.0) {
     13:     PetscCall(BFGSKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y));
     14:   } else {
     15:     PetscCall(SymBroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y, PETSC_TRUE));
     16:   }
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: static PetscErrorCode MatSolve_LMVMSymBadBrdn_Recursive(Mat B, Vec X, Vec Y)
     21: {
     22:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
     23:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

     25:   PetscFunctionBegin;
     26:   PetscCheck(lsb->psi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Psi must first be set using MatLMVMSymBadBroydenSetPsi()");
     27:   if (lsb->psi_scalar == 0.0) {
     28:     PetscCall(BFGSKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y));
     29:   } else if (lsb->psi_scalar == 1.0) {
     30:     PetscCall(DFPKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y));
     31:   } else {
     32:     PetscCall(SymBroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y, PETSC_FALSE));
     33:   }
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode MatMult_LMVMSymBadBrdn_CompactDense(Mat B, Vec X, Vec Y)
     38: {
     39:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
     40:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

     42:   PetscFunctionBegin;
     43:   PetscCheck(lsb->psi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Psi must first be set using MatLMVMSymBadBroydenSetPsi()");
     44:   if (lsb->psi_scalar == 0.0) {
     45:     PetscCall(DFPKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Y));
     46:   } else if (lsb->psi_scalar == 1.0) {
     47:     PetscCall(BFGSKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Y));
     48:   } else {
     49:     PetscCall(SymBroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Y, PETSC_TRUE));
     50:   }
     51:   PetscFunctionReturn(PETSC_SUCCESS);
     52: }

     54: static PetscErrorCode MatSolve_LMVMSymBadBrdn_CompactDense(Mat B, Vec X, Vec Y)
     55: {
     56:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
     57:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

     59:   PetscFunctionBegin;
     60:   PetscCheck(lsb->psi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Psi must first be set using MatLMVMSymBadBroydenSetPsi()");
     61:   if (lsb->psi_scalar == 0.0) {
     62:     PetscCall(BFGSKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, Y));
     63:   } else if (lsb->psi_scalar == 1.0) {
     64:     PetscCall(DFPKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, Y));
     65:   } else {
     66:     PetscCall(SymBroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, Y, PETSC_FALSE));
     67:   }
     68:   PetscFunctionReturn(PETSC_SUCCESS);
     69: }

     71: static PetscErrorCode MatLMVMSetMultAlgorithm_SymBadBrdn(Mat B)
     72: {
     73:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;

     75:   PetscFunctionBegin;
     76:   switch (lmvm->mult_alg) {
     77:   case MAT_LMVM_MULT_RECURSIVE:
     78:     lmvm->ops->mult  = MatMult_LMVMSymBadBrdn_Recursive;
     79:     lmvm->ops->solve = MatSolve_LMVMSymBadBrdn_Recursive;
     80:     break;
     81:   case MAT_LMVM_MULT_DENSE:
     82:   case MAT_LMVM_MULT_COMPACT_DENSE:
     83:     lmvm->ops->mult  = MatMult_LMVMSymBadBrdn_CompactDense;
     84:     lmvm->ops->solve = MatSolve_LMVMSymBadBrdn_CompactDense;
     85:     break;
     86:   }
     87:   lmvm->ops->multht  = lmvm->ops->mult;
     88:   lmvm->ops->solveht = lmvm->ops->solve;
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: static PetscErrorCode MatSetFromOptions_LMVMSymBadBrdn(Mat B, PetscOptionItems PetscOptionsObject)
     93: {
     94:   Mat_LMVM                  *lmvm = (Mat_LMVM *)B->data;
     95:   Mat_SymBrdn               *lsb  = (Mat_SymBrdn *)lmvm->ctx;
     96:   MatLMVMSymBroydenScaleType stype;

     98:   PetscFunctionBegin;
     99:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
    100:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Bad Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBADBROYDEN)");
    101:   PetscCall(PetscOptionsReal("-mat_lmvm_psi", "convex ratio between DFP and BFGS components of the update", "", lsb->psi_scalar, &lsb->psi_scalar, NULL));
    102:   PetscCheck(lsb->psi_scalar >= 0.0 && lsb->psi_scalar <= 1.0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
    103:   PetscCall(SymBroydenRescaleSetFromOptions(B, lsb->rescale, PetscOptionsObject));
    104:   PetscOptionsHeadEnd();
    105:   PetscCall(SymBroydenRescaleGetType(lsb->rescale, &stype));
    106:   if (stype == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(SymBroydenRescaleSetDiagonalMode(lsb->rescale, PETSC_TRUE));
    107:   PetscFunctionReturn(PETSC_SUCCESS);
    108: }

    110: PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat B)
    111: {
    112:   Mat_LMVM    *lmvm;
    113:   Mat_SymBrdn *lsb;

    115:   PetscFunctionBegin;
    116:   PetscCall(MatCreate_LMVMSymBrdn(B));
    117:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBROYDEN));
    118:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;

    120:   lmvm                        = (Mat_LMVM *)B->data;
    121:   lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_SymBadBrdn;

    123:   lsb                   = (Mat_SymBrdn *)lmvm->ctx;
    124:   lsb->psi_scalar       = 0.875;
    125:   lsb->phi_scalar       = PETSC_DETERMINE;
    126:   lsb->rescale->forward = PETSC_FALSE;
    127:   lsb->rescale->theta   = 1.0 - lsb->psi_scalar;

    129:   PetscCall(MatLMVMSetMultAlgorithm_SymBadBrdn(B));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*@
    134:   MatCreateLMVMSymBadBroyden - Creates a limited-memory Symmetric "Bad" Broyden-type matrix used
    135:   for approximating Jacobians.

    137:   Collective

    139:   Input Parameters:
    140: + comm - MPI communicator
    141: . n    - number of local rows for storage vectors
    142: - N    - global size of the storage vectors

    144:   Output Parameter:
    145: . B - the matrix

    147:   Options Database Keys:
    148: + -mat_lmvm_hist_size         - the number of history vectors to keep
    149: . -mat_lmvm_psi               - convex ratio between BFGS and DFP components of the update
    150: . -mat_lmvm_scale_type        - type of scaling applied to J0 (none, scalar, diagonal)
    151: . -mat_lmvm_mult_algorithm    - the algorithm to use for multiplication (recursive, dense, compact_dense)
    152: . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed
    153: . -mat_lmvm_eps               - (developer) numerical zero tolerance for testing when an update should be skipped
    154: . -mat_lmvm_debug             - (developer) perform internal debugging checks
    155: . -mat_lmvm_theta             - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
    156: . -mat_lmvm_rho               - (developer) update limiter for the J0 scaling
    157: . -mat_lmvm_alpha             - (developer) coefficient factor for the quadratic subproblem in J0 scaling
    158: . -mat_lmvm_beta              - (developer) exponential factor for the diagonal J0 scaling
    159: - -mat_lmvm_sigma_hist        - (developer) number of past updates to use in J0 scaling

    161:   Level: intermediate

    163:   Notes:
    164:   L-SymBadBrdn is a convex combination of L-DFP and L-BFGS such that $B^{-1} = (1 - \psi)*B_{\text{DFP}}^{-1} +
    165:   \psi*B_{\text{BFGS}}^{-1}$. The combination factor $\psi$ is restricted to the range $[0, 1]$, where the L-SymBadBrdn matrix
    166:   is guaranteed to be symmetric positive-definite. Note that this combination is on the inverses and not on the
    167:   forwards. For forward convex combinations, use the L-SymBrdn matrix (`MATLMVMSYMBROYDEN`).

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

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

    176: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
    177:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMBadBroyden()`
    178: @*/
    179: PetscErrorCode MatCreateLMVMSymBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
    180: {
    181:   PetscFunctionBegin;
    182:   PetscCall(KSPInitializePackage());
    183:   PetscCall(MatCreate(comm, B));
    184:   PetscCall(MatSetSizes(*B, n, n, N, N));
    185:   PetscCall(MatSetType(*B, MATLMVMSYMBADBROYDEN));
    186:   PetscCall(MatSetUp(*B));
    187:   PetscFunctionReturn(PETSC_SUCCESS);
    188: }