Actual source code: lmvm.h
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"
8: PETSC_INTERN PetscLogEvent MATLMVM_Update;
10: /*
11: MATLMVM format - a base matrix-type that represents Limited-Memory
12: Variable Metric (LMVM) approximations of a Jacobian.
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: */
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: };
32: typedef PetscInt MatLMVMResetMode;
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)
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: };
51: /* Identifies vector bases used internally.
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
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 vaues, 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: };
70: typedef PetscInt MatLMVMBasisType;
72: typedef enum {
73: MATLMVM_MODE_PRIMAL = 0,
74: MATLMVM_MODE_DUAL = 1,
75: } MatLMVMMode;
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)
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;
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;
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;
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);
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);
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);
137: PETSC_INTERN PetscErrorCode MatLMVMGetJ0InvDiag(Mat, Vec *);
138: PETSC_INTERN PetscErrorCode MatLMVMRestoreJ0InvDiag(Mat, Vec *);
140: PETSC_INTERN PetscErrorCode MatLMVMGetRange(Mat, PetscInt *, PetscInt *);
142: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0HermitianTranspose(Mat, Vec, Vec);
143: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvHermitianTranspose(Mat, Vec, Vec);
145: PETSC_INTERN PetscErrorCode MatLMVMGetJ0Scalar(Mat, PetscBool *, PetscScalar *);
146: PETSC_INTERN PetscErrorCode MatLMVMJ0KSPIsExact(Mat, PetscBool *);
148: PETSC_INTERN PetscErrorCode MatLMVMUseVecLayoutsIfCompatible(Mat, Vec, Vec);
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);
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 *);
163: PETSC_INTERN PetscBool ByrdNocedalSchnabelCite;
164: PETSC_INTERN const char ByrdNocedalSchnabelCitation[];