1 #pragma once 2 #include <petscksp.h> 3 #include <petsc/private/matimpl.h> 4 #include <petsc/private/vecimpl.h> 5 #include "lmbasis.h" 6 #include "lmproducts.h" 7 8 PETSC_INTERN PetscLogEvent MATLMVM_Update; 9 10 /* 11 MATLMVM format - a base matrix-type that represents Limited-Memory 12 Variable Metric (LMVM) approximations of a Jacobian. 13 14 LMVM approximations can be symmetric, symmetric positive-definite, 15 rectangular, or otherwise square with no determinable properties. 16 Each derived LMVM type should automatically set its matrix properties 17 if its construction can guarantee symmetry (MAT_SYMMETRIC) or symmetric 18 positive-definiteness (MAT_SPD). 19 */ 20 21 /* MatLMVMReset(Mat, PetscBool) has a simple boolean for destructive/nondestructive reset, 22 but internally it is helpful to have more control 23 */ 24 enum { 25 MAT_LMVM_RESET_HISTORY = 0x0, 26 MAT_LMVM_RESET_BASES = 0x1, 27 MAT_LMVM_RESET_J0 = 0x2, 28 MAT_LMVM_RESET_VECS = 0x4, 29 MAT_LMVM_RESET_ALL = 0xf, 30 }; 31 32 typedef PetscInt MatLMVMResetMode; 33 34 #define MatLMVMResetClearsBases(mode) ((mode) & MAT_LMVM_RESET_BASES) 35 #define MatLMVMResetClearsJ0(mode) ((mode) & MAT_LMVM_RESET_J0) 36 #define MatLMVMResetClearsVecs(mode) ((mode) & MAT_LMVM_RESET_VECS) 37 #define MatLMVMResetClearsAll(mode) ((mode) == MAT_LMVM_RESET_ALL) 38 39 typedef struct _MatOps_LMVM *MatOps_LMVM; 40 struct _MatOps_LMVM { 41 PetscErrorCode (*update)(Mat, Vec, Vec); 42 PetscErrorCode (*reset)(Mat, MatLMVMResetMode); 43 PetscErrorCode (*mult)(Mat, Vec, Vec); 44 PetscErrorCode (*multht)(Mat, Vec, Vec); 45 PetscErrorCode (*solve)(Mat, Vec, Vec); 46 PetscErrorCode (*solveht)(Mat, Vec, Vec); 47 PetscErrorCode (*copy)(Mat, Mat, MatStructure); 48 PetscErrorCode (*setmultalgorithm)(Mat); 49 }; 50 51 /* Identifies vector bases used internally. 52 53 - bases are listed in "dual pairs": e.g. if an algorithm for MatMult() uses basis X somewhere, then the dual 54 algorithm for MatSolve will use the dual basis X' in the same place 55 56 - bases are stored in the `basis[]` array in Mat_LMVM, rather than as struct members with names, to enable "modal" 57 access: code of the form `basis[LMVMModeMap(LMBASIS_S, mode)]` can be used for 58 the primal algorithm (`mode == MATLMVM_MODE_PRIMAL`) and the dual algorithm (`mode == MATLMMV_MODE_DUAL`). 59 */ 60 enum { 61 LMBASIS_S = 0, // differences between solutions, S_i = (X_{i+1} - X_i) 62 LMBASIS_Y = 1, // differences in function values, Y_i = (F_{i+1} - F_i) 63 LMBASIS_H0Y = 2, // H_0 = J_0^{-1} 64 LMBASIS_B0S = 3, // B_0 is the symbol used instead of J_0 in many textbooks and papers, we use it internally 65 LMBASIS_S_MINUS_H0Y = 4, 66 LMBASIS_Y_MINUS_B0S = 5, 67 LMBASIS_END 68 }; 69 70 typedef PetscInt MatLMVMBasisType; 71 72 typedef enum { 73 MATLMVM_MODE_PRIMAL = 0, 74 MATLMVM_MODE_DUAL = 1, 75 } MatLMVMMode; 76 77 #define LMVMModeMap(a, mode) ((a) ^ (PetscInt)(mode)) 78 #define MatLMVMApplyJ0Mode(mode) ((mode) == MATLMVM_MODE_PRIMAL ? MatLMVMApplyJ0Fwd : MatLMVMApplyJ0Inv) 79 #define MatLMVMApplyJ0HermitianTransposeMode(mode) ((mode) == MATLMVM_MODE_PRIMAL ? MatLMVMApplyJ0HermitianTranspose : MatLMVMApplyJ0InvHermitianTranspose) 80 #define MatLMVMBasisSizeOf(type) ((type) & LMBASIS_Y) 81 82 typedef struct { 83 /* Core data structures for stored updates */ 84 struct _MatOps_LMVM ops[1]; 85 PetscBool prev_set; 86 PetscInt m, k, nupdates, nrejects, nresets; 87 LMBasis basis[LMBASIS_END]; 88 LMProducts products[LMBLOCK_END][LMBASIS_END][LMBASIS_END]; 89 Vec Xprev, Fprev; 90 91 /* User-defined initial Jacobian tools */ 92 PetscScalar shift; 93 Mat J0; 94 KSP J0ksp; 95 PetscBool disable_ksp_viewers; 96 PetscBool created_J0; 97 PetscBool created_J0ksp; 98 PetscBool do_not_cache_J0_products; 99 PetscBool cache_gradient_products; 100 101 /* Miscellenous parameters */ 102 PetscReal eps; /* (default: PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0)) */ 103 MatLMVMMultAlgorithm mult_alg; 104 void *ctx; /* implementation specific context */ 105 PetscBool debug; 106 } Mat_LMVM; 107 108 /* Shared internal functions for LMVM matrices */ 109 PETSC_INTERN PetscErrorCode MatUpdateKernel_LMVM(Mat, Vec, Vec); 110 PETSC_INTERN PetscErrorCode MatUpdate_LMVM(Mat, Vec, Vec); 111 PETSC_INTERN PetscErrorCode MatAllocate_LMVM(Mat, Vec, Vec); 112 PETSC_INTERN PetscErrorCode MatLMVMAllocateBases(Mat); 113 PETSC_INTERN PetscErrorCode MatLMVMAllocateVecs(Mat); 114 PETSC_INTERN PetscErrorCode MatLMVMReset_Internal(Mat, MatLMVMResetMode); 115 PETSC_INTERN PetscErrorCode MatReset_LMVM(Mat, MatLMVMResetMode); 116 117 /* LMVM implementations of core Mat functionality */ 118 PETSC_INTERN PetscErrorCode MatSetFromOptions_LMVM(Mat, PetscOptionItems PetscOptionsObject); 119 PETSC_INTERN PetscErrorCode MatSetUp_LMVM(Mat); 120 PETSC_INTERN PetscErrorCode MatView_LMVM(Mat, PetscViewer); 121 PETSC_INTERN PetscErrorCode MatDestroy_LMVM(Mat); 122 PETSC_INTERN PetscErrorCode MatCreate_LMVM(Mat); 123 124 /* Create functions for derived LMVM types */ 125 PETSC_EXTERN PetscErrorCode MatCreate_LMVMDFP(Mat); 126 PETSC_EXTERN PetscErrorCode MatCreate_LMVMDDFP(Mat); 127 PETSC_EXTERN PetscErrorCode MatCreate_LMVMBFGS(Mat); 128 PETSC_EXTERN PetscErrorCode MatCreate_LMVMDBFGS(Mat); 129 PETSC_EXTERN PetscErrorCode MatCreate_LMVMDQN(Mat); 130 PETSC_EXTERN PetscErrorCode MatCreate_LMVMSR1(Mat); 131 PETSC_EXTERN PetscErrorCode MatCreate_LMVMBrdn(Mat); 132 PETSC_EXTERN PetscErrorCode MatCreate_LMVMBadBrdn(Mat); 133 PETSC_EXTERN PetscErrorCode MatCreate_LMVMSymBrdn(Mat); 134 PETSC_EXTERN PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat); 135 PETSC_EXTERN PetscErrorCode MatCreate_LMVMDiagBrdn(Mat); 136 137 PETSC_INTERN PetscErrorCode MatLMVMGetJ0InvDiag(Mat, Vec *); 138 PETSC_INTERN PetscErrorCode MatLMVMRestoreJ0InvDiag(Mat, Vec *); 139 140 PETSC_INTERN PetscErrorCode MatLMVMGetRange(Mat, PetscInt *, PetscInt *); 141 142 PETSC_INTERN PetscErrorCode MatLMVMApplyJ0HermitianTranspose(Mat, Vec, Vec); 143 PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvHermitianTranspose(Mat, Vec, Vec); 144 145 PETSC_INTERN PetscErrorCode MatLMVMGetJ0Scalar(Mat, PetscBool *, PetscScalar *); 146 PETSC_INTERN PetscErrorCode MatLMVMJ0KSPIsExact(Mat, PetscBool *); 147 148 PETSC_INTERN PetscErrorCode MatLMVMUseVecLayoutsIfCompatible(Mat, Vec, Vec); 149 150 PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedBasis(Mat, MatLMVMBasisType, LMBasis *, MatLMVMBasisType *, PetscScalar *); 151 PETSC_INTERN PetscErrorCode MatLMVMGetWorkRow(Mat, Vec *); 152 PETSC_INTERN PetscErrorCode MatLMVMRestoreWorkRow(Mat, Vec *); 153 PETSC_INTERN PetscErrorCode MatLMVMBasisGetVecRead(Mat, MatLMVMBasisType, PetscInt, Vec *, PetscScalar *); 154 PETSC_INTERN PetscErrorCode MatLMVMBasisRestoreVecRead(Mat, MatLMVMBasisType, PetscInt, Vec *, PetscScalar *); 155 PETSC_INTERN PetscErrorCode MatLMVMBasisGEMVH(Mat, MatLMVMBasisType, PetscInt, PetscInt, PetscScalar, Vec, PetscScalar, Vec); 156 PETSC_INTERN PetscErrorCode MatLMVMBasisGEMV(Mat, MatLMVMBasisType, PetscInt, PetscInt, PetscScalar, Vec, PetscScalar, Vec); 157 158 PETSC_INTERN PetscErrorCode MatLMVMCreateProducts(Mat, LMBlockType, LMProducts *); 159 PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedProducts(Mat, MatLMVMBasisType, MatLMVMBasisType, LMBlockType, LMProducts *); 160 PETSC_INTERN PetscErrorCode MatLMVMProductsInsertDiagonalValue(Mat, MatLMVMBasisType, MatLMVMBasisType, PetscInt, PetscScalar); 161 PETSC_INTERN PetscErrorCode MatLMVMProductsGetDiagonalValue(Mat, MatLMVMBasisType, MatLMVMBasisType, PetscInt, PetscScalar *); 162 163 PETSC_INTERN PetscBool ByrdNocedalSchnabelCite; 164 PETSC_INTERN const char ByrdNocedalSchnabelCitation[]; 165