1 #include "symbrdn.h" /*I "petscksp.h" I*/
2
MatMult_LMVMSymBadBrdn_Recursive(Mat B,Vec X,Vec Y)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;
7
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 }
19
MatSolve_LMVMSymBadBrdn_Recursive(Mat B,Vec X,Vec Y)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;
24
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 }
36
MatMult_LMVMSymBadBrdn_CompactDense(Mat B,Vec X,Vec Y)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;
41
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 }
53
MatSolve_LMVMSymBadBrdn_CompactDense(Mat B,Vec X,Vec Y)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;
58
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 }
70
MatLMVMSetMultAlgorithm_SymBadBrdn(Mat B)71 static PetscErrorCode MatLMVMSetMultAlgorithm_SymBadBrdn(Mat B)
72 {
73 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
74
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 }
91
MatSetFromOptions_LMVMSymBadBrdn(Mat B,PetscOptionItems PetscOptionsObject)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;
97
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 }
109
MatCreate_LMVMSymBadBrdn(Mat B)110 PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat B)
111 {
112 Mat_LMVM *lmvm;
113 Mat_SymBrdn *lsb;
114
115 PetscFunctionBegin;
116 PetscCall(MatCreate_LMVMSymBrdn(B));
117 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBROYDEN));
118 B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;
119
120 lmvm = (Mat_LMVM *)B->data;
121 lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_SymBadBrdn;
122
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;
128
129 PetscCall(MatLMVMSetMultAlgorithm_SymBadBrdn(B));
130 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 /*@
134 MatCreateLMVMSymBadBroyden - Creates a limited-memory Symmetric "Bad" Broyden-type matrix used
135 for approximating Jacobians.
136
137 Collective
138
139 Input Parameters:
140 + comm - MPI communicator
141 . n - number of local rows for storage vectors
142 - N - global size of the storage vectors
143
144 Output Parameter:
145 . B - the matrix
146
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
160
161 Level: intermediate
162
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`).
168
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.
172
173 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()` paradigm instead of this
174 routine directly.
175
176 .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
177 `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMBadBroyden()`
178 @*/
MatCreateLMVMSymBadBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)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 }
189