1 #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h> /*I "petscksp.h" I*/
2
3 /* The DFP update can be written as
4
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
6
7 So B_{k+1}x can be computed in the following way
8
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 */
DFPKernel_Recursive(Mat B,MatLMVMMode mode,Vec X,Vec BX)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;
22
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;
36
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;
50
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 }
67
DFPKernel_CompactDense(Mat B,MatLMVMMode mode,Vec X,Vec BX)68 PETSC_INTERN PetscErrorCode DFPKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
69 {
70 PetscInt oldest, next;
71
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;
83
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));
89
90 PetscCall(MatLMVMGetWorkRow(B, &StB0X));
91 PetscCall(MatLMVMGetWorkRow(B, &YtX));
92 PetscCall(MatLMVMGetWorkRow(B, &u));
93 PetscCall(MatLMVMGetWorkRow(B, &v));
94
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));
98
99 PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
100 PetscCall(LMProductsSolve(YtS, oldest, next, YtX, YtX, /* ^H */ PETSC_FALSE));
101
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));
106
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));
110
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 }
118
DFPKernel_Dense(Mat B,MatLMVMMode mode,Vec X,Vec BX)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;
126
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);
132
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));
139
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 }
158
MatMult_LMVMDFP_Recursive(Mat B,Vec X,Vec BX)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 }
165
MatMult_LMVMDFP_CompactDense(Mat B,Vec X,Vec BX)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 }
172
MatMult_LMVMDFP_Dense(Mat B,Vec X,Vec BX)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 }
179
MatSolve_LMVMDFP_Recursive(Mat B,Vec X,Vec HX)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 }
186
MatSolve_LMVMDFP_CompactDense(Mat B,Vec X,Vec HX)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 }
193
MatSetFromOptions_LMVMDFP(Mat B,PetscOptionItems PetscOptionsObject)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;
198
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 }
206
MatLMVMSetMultAlgorithm_DFP(Mat B)207 static PetscErrorCode MatLMVMSetMultAlgorithm_DFP(Mat B)
208 {
209 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
210
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 }
230
MatCreate_LMVMDFP(Mat B)231 PetscErrorCode MatCreate_LMVMDFP(Mat B)
232 {
233 Mat_LMVM *lmvm;
234 Mat_SymBrdn *ldfp;
235
236 PetscFunctionBegin;
237 PetscCall(MatCreate_LMVMSymBrdn(B));
238 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP));
239 B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
240
241 lmvm = (Mat_LMVM *)B->data;
242 lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_DFP;
243 PetscCall(MatLMVMSetMultAlgorithm_DFP(B));
244
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 }
251
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.
256
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.
261
262 Collective
263
264 Input Parameters:
265 + comm - MPI communicator
266 . n - number of local rows for storage vectors
267 - N - global size of the storage vectors
268
269 Output Parameter:
270 . B - the matrix
271
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
279
280 Level: intermediate
281
282 Note:
283 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
284 paradigm instead of this routine directly.
285
286 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDFP`, `MatCreateLMVMBFGS()`, `MatCreateLMVMSR1()`,
287 `MatCreateLMVMBroyden()`, `MatCreateLMVMBadBroyden()`, `MatCreateLMVMSymBroyden()`
288 @*/
MatCreateLMVMDFP(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)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 }
299