1 #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h> /*I "petscksp.h" I*/
2 #include <petscblaslapack.h>
3
4 // Broyden is dual to bad Broyden: MatSolve routines are dual to bad Broyden MatMult routines
5
MatSolve_LMVMBrdn_Recursive(Mat B,Vec F,Vec dX)6 static PetscErrorCode MatSolve_LMVMBrdn_Recursive(Mat B, Vec F, Vec dX)
7 {
8 PetscFunctionBegin;
9 PetscCall(BadBroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
10 PetscFunctionReturn(PETSC_SUCCESS);
11 }
12
MatSolveHermitianTranspose_LMVMBrdn_Recursive(Mat B,Vec F,Vec dX)13 static PetscErrorCode MatSolveHermitianTranspose_LMVMBrdn_Recursive(Mat B, Vec F, Vec dX)
14 {
15 PetscFunctionBegin;
16 PetscCall(BadBroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
MatSolve_LMVMBrdn_CompactDense(Mat B,Vec F,Vec dX)20 static PetscErrorCode MatSolve_LMVMBrdn_CompactDense(Mat B, Vec F, Vec dX)
21 {
22 PetscFunctionBegin;
23 PetscCall(BadBroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
MatSolveHermitianTranspose_LMVMBrdn_CompactDense(Mat B,Vec F,Vec dX)27 static PetscErrorCode MatSolveHermitianTranspose_LMVMBrdn_CompactDense(Mat B, Vec F, Vec dX)
28 {
29 PetscFunctionBegin;
30 PetscCall(BadBroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
34 /* The Broyden kernel can be written as
35
36 B_k = B_{k-1} + (y_{k-1} - B_{k-1} s_{k-1}) (s_{k-1}^H s_{k-1})^{-1} s_{k-1}^H
37
38 This means that we can write B_k x using an intermediate variable alpha_k as
39
40 alpha_k = (s_{k-1}^H s_{k-1})^{-1} s_{k-1}^H x
41 B_k x = y_{k-1} alpha_k + B_{k-1}(x - s_{k-1} alpha_k)
42 \_____________ ____________/
43 V
44 recursive application
45 */
46
BroydenKernel_Recursive(Mat B,MatLMVMMode mode,Vec X,Vec BX)47 PETSC_INTERN PetscErrorCode BroydenKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BX)
48 {
49 Vec G = X;
50 Vec W = BX;
51 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
52 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
53 LMBasis S = NULL, Y = NULL;
54 PetscInt oldest, next;
55
56 PetscFunctionBegin;
57 PetscCall(MatLMVMGetRange(B, &oldest, &next));
58 if (next > oldest) {
59 LMProducts StS;
60
61 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
62 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
63 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_DIAGONAL, &StS));
64
65 PetscCall(LMBasisGetWorkVec(S, &G));
66 PetscCall(VecCopy(X, G));
67 PetscCall(VecZeroEntries(BX));
68
69 for (PetscInt i = next - 1; i >= oldest; i--) {
70 Vec s_i, y_i;
71 PetscScalar sitg, sitsi, alphai;
72
73 PetscCall(LMBasisGetVecRead(S, i, &s_i));
74 PetscCall(LMBasisGetVecRead(Y, i, &y_i));
75 PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));
76
77 PetscCall(VecDot(G, s_i, &sitg));
78 alphai = sitg / sitsi;
79 PetscCall(VecAXPY(BX, alphai, y_i));
80 PetscCall(VecAXPY(G, -alphai, s_i));
81
82 PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
83 PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
84 }
85
86 PetscCall(LMBasisGetWorkVec(Y, &W));
87 }
88 PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, W));
89 if (next > oldest) {
90 PetscCall(VecAXPY(BX, 1.0, W));
91 PetscCall(LMBasisRestoreWorkVec(Y, &W));
92 PetscCall(LMBasisRestoreWorkVec(S, &G));
93 }
94 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
97 /* The adjoint of the kernel is
98
99 B_k^H = B_{k-1}^H + s_{k-1} (s_{k-1}^H s_{k-1})^{-1} (y_{k-1} - B_{k-1} s_{k-1})^H
100
101 This means that we can write B_k^H x using an intermediate variable w_k = B_{k-1}^H x
102
103 w_k = B_{k-1}^H x <-- recursive application
104 B_k^H x = w_k + s_{k-1} (s_{k-1}^H s_{k-1})^{-1} (y_{k-1}^H x - s_k^H w_k)
105 */
106
BroydenKernelHermitianTranspose_Recursive(Mat B,MatLMVMMode mode,Vec X,Vec BHX)107 PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
108 {
109 PetscInt oldest, next;
110
111 PetscFunctionBegin;
112 PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
113 PetscCall(MatLMVMGetRange(B, &oldest, &next));
114 if (next > oldest) {
115 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
116 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
117 LMBasis S, Y;
118 LMProducts StS;
119
120 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_DIAGONAL, &StS));
121 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
122 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
123
124 for (PetscInt i = oldest; i < next; i++) {
125 Vec s_i, y_i;
126 PetscScalar sitBHX, sitsi, yitx;
127
128 PetscCall(LMBasisGetVecRead(S, i, &s_i));
129 PetscCall(LMBasisGetVecRead(Y, i, &y_i));
130 PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));
131
132 PetscCall(VecDotBegin(BHX, s_i, &sitBHX));
133 PetscCall(VecDotBegin(X, y_i, &yitx));
134 PetscCall(VecDotEnd(BHX, s_i, &sitBHX));
135 PetscCall(VecDotEnd(X, y_i, &yitx));
136 PetscCall(VecAXPY(BHX, (yitx - sitBHX) / sitsi, s_i));
137 PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
138 PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
139 }
140 }
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
144 /*
145 The Broyden kernel can be written as
146
147 B_k = B_0 + (Y - B_0 S) (triu(S^H S))^{-1} S^H
148
149 where triu is the upper triangular component. We solve by back substitution each time
150 we apply
151 */
152
BroydenKernel_CompactDense(Mat B,MatLMVMMode mode,Vec X,Vec BX)153 PETSC_INTERN PetscErrorCode BroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
154 {
155 PetscInt oldest, next;
156
157 PetscFunctionBegin;
158 PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX));
159 PetscCall(MatLMVMGetRange(B, &oldest, &next));
160 if (next > oldest) {
161 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
162 MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
163 LMBasis S;
164 LMProducts StS;
165 Vec StX;
166
167 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
168 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
169 PetscCall(MatLMVMGetWorkRow(B, &StX));
170 PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, X, 0.0, StX));
171 PetscCall(LMProductsSolve(StS, oldest, next, StX, StX, /* ^H */ PETSC_FALSE));
172 PetscCall(MatLMVMBasisGEMV(B, Y_minus_B0S_t, oldest, next, 1.0, StX, 1.0, BX));
173 PetscCall(MatLMVMRestoreWorkRow(B, &StX));
174 }
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /*
179 The adjoint of the above formula is
180
181 B_k^H = B_0^H + S^H (triu(S^H S))^{-H} (Y - B_0 S)^H
182 */
183
BroydenKernelHermitianTranspose_CompactDense(Mat B,MatLMVMMode mode,Vec X,Vec BHX)184 PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
185 {
186 PetscInt oldest, next;
187
188 PetscFunctionBegin;
189 PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
190 PetscCall(MatLMVMGetRange(B, &oldest, &next));
191 if (next > oldest) {
192 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
193 MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
194 LMProducts StS;
195 LMBasis S;
196 Vec YmB0StX;
197
198 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
199 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
200 PetscCall(MatLMVMGetWorkRow(B, &YmB0StX));
201 PetscCall(MatLMVMBasisGEMVH(B, Y_minus_B0S_t, oldest, next, 1.0, X, 0.0, YmB0StX));
202 PetscCall(LMProductsSolve(StS, oldest, next, YmB0StX, YmB0StX, /* ^H */ PETSC_TRUE));
203 PetscCall(LMBasisGEMV(S, oldest, next, 1.0, YmB0StX, 1.0, BHX));
204 PetscCall(MatLMVMRestoreWorkRow(B, &YmB0StX));
205 }
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
209 /*
210 The dense Broyden formula use for CompactDense can be written in a non-update form as
211
212 B_k = [B_0 | Y] [ I | -S ] [ I ]
213 [---+----] [----------------------]
214 [ 0 | I ] [ triu(S^H S)^{-1} S^H ]
215
216 The advantage of this form is B_0 appears only once and the (Y - B_0 S) vectors are not needed
217 */
218
BroydenKernel_Dense(Mat B,MatLMVMMode mode,Vec X,Vec BX)219 PETSC_INTERN PetscErrorCode BroydenKernel_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
220 {
221 Vec G = X;
222 Vec W = BX;
223 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
224 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
225 LMBasis S = NULL, Y = NULL;
226 PetscInt oldest, next;
227
228 PetscFunctionBegin;
229 PetscCall(MatLMVMGetRange(B, &oldest, &next));
230 if (next > oldest) {
231 LMProducts StS;
232 Vec StX;
233
234 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
235 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
236 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
237 PetscCall(LMBasisGetWorkVec(S, &G));
238 PetscCall(MatLMVMGetWorkRow(B, &StX));
239 PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, X, 0.0, StX));
240 PetscCall(LMProductsSolve(StS, oldest, next, StX, StX, PETSC_FALSE));
241 PetscCall(VecCopy(X, G));
242 PetscCall(LMBasisGEMV(S, oldest, next, -1.0, StX, 1.0, G));
243 PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, StX, 0.0, BX));
244 PetscCall(MatLMVMRestoreWorkRow(B, &StX));
245 PetscCall(LMBasisGetWorkVec(Y, &W));
246 }
247 PetscCall(MatLMVMApplyJ0Mode(mode)(B, G, W));
248 if (next > oldest) {
249 PetscCall(VecAXPY(BX, 1.0, W));
250 PetscCall(LMBasisRestoreWorkVec(Y, &W));
251 PetscCall(LMBasisRestoreWorkVec(S, &G));
252 }
253 PetscFunctionReturn(PETSC_SUCCESS);
254 }
255
256 /*
257 The adoint of the above formula is
258
259 B_k^h = [I | S triu(S^H S)^{-H} ] [ I | 0 ] [ B_0^H ]
260 [------+---] [-------]
261 [ -S^H | I ] [ Y^H ]
262 */
263
BroydenKernelHermitianTranspose_Dense(Mat B,MatLMVMMode mode,Vec X,Vec BHX)264 PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
265 {
266 PetscInt oldest, next;
267
268 PetscFunctionBegin;
269 PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, X, BHX));
270 PetscCall(MatLMVMGetRange(B, &oldest, &next));
271 if (next > oldest) {
272 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
273 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
274 LMBasis S, Y;
275 LMProducts StS;
276 Vec YtX, StBHX;
277
278 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, S_t, LMBLOCK_UPPER_TRIANGLE, &StS));
279 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
280 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
281 PetscCall(MatLMVMGetWorkRow(B, &YtX));
282 PetscCall(MatLMVMGetWorkRow(B, &StBHX));
283 PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, BHX, 0.0, StBHX));
284 PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
285 PetscCall(VecAXPY(YtX, -1.0, StBHX));
286 PetscCall(LMProductsSolve(StS, oldest, next, YtX, YtX, PETSC_TRUE));
287 PetscCall(LMBasisGEMV(S, oldest, next, 1.0, YtX, 1.0, BHX));
288 PetscCall(MatLMVMRestoreWorkRow(B, &StBHX));
289 PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
290 }
291 PetscFunctionReturn(PETSC_SUCCESS);
292 }
293
MatMult_LMVMBrdn_Recursive(Mat B,Vec X,Vec Z)294 static PetscErrorCode MatMult_LMVMBrdn_Recursive(Mat B, Vec X, Vec Z)
295 {
296 PetscFunctionBegin;
297 PetscCall(BroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Z));
298 PetscFunctionReturn(PETSC_SUCCESS);
299 }
300
MatMultHermitianTranspose_LMVMBrdn_Recursive(Mat B,Vec X,Vec Z)301 static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_Recursive(Mat B, Vec X, Vec Z)
302 {
303 PetscFunctionBegin;
304 PetscCall(BroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_PRIMAL, X, Z));
305 PetscFunctionReturn(PETSC_SUCCESS);
306 }
307
MatMult_LMVMBrdn_CompactDense(Mat B,Vec X,Vec Z)308 static PetscErrorCode MatMult_LMVMBrdn_CompactDense(Mat B, Vec X, Vec Z)
309 {
310 PetscFunctionBegin;
311 PetscCall(BroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Z));
312 PetscFunctionReturn(PETSC_SUCCESS);
313 }
314
MatMultHermitianTranspose_LMVMBrdn_CompactDense(Mat B,Vec X,Vec Z)315 static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_CompactDense(Mat B, Vec X, Vec Z)
316 {
317 PetscFunctionBegin;
318 PetscCall(BroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_PRIMAL, X, Z));
319 PetscFunctionReturn(PETSC_SUCCESS);
320 }
321
MatMult_LMVMBrdn_Dense(Mat B,Vec X,Vec Z)322 PETSC_UNUSED static PetscErrorCode MatMult_LMVMBrdn_Dense(Mat B, Vec X, Vec Z)
323 {
324 PetscFunctionBegin;
325 PetscCall(BroydenKernel_Dense(B, MATLMVM_MODE_PRIMAL, X, Z));
326 PetscFunctionReturn(PETSC_SUCCESS);
327 }
328
MatMultHermitianTranspose_LMVMBrdn_Dense(Mat B,Vec X,Vec Z)329 PETSC_UNUSED static PetscErrorCode MatMultHermitianTranspose_LMVMBrdn_Dense(Mat B, Vec X, Vec Z)
330 {
331 PetscFunctionBegin;
332 PetscCall(BroydenKernelHermitianTranspose_Dense(B, MATLMVM_MODE_PRIMAL, X, Z));
333 PetscFunctionReturn(PETSC_SUCCESS);
334 }
335
MatUpdate_LMVMBrdn(Mat B,Vec X,Vec F)336 static PetscErrorCode MatUpdate_LMVMBrdn(Mat B, Vec X, Vec F)
337 {
338 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
339 Mat_Brdn *brdn = (Mat_Brdn *)lmvm->ctx;
340 PetscBool cache_YtFprev = (lmvm->mult_alg != MAT_LMVM_MULT_RECURSIVE) ? lmvm->cache_gradient_products : PETSC_FALSE;
341
342 PetscFunctionBegin;
343 if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
344 if (lmvm->prev_set) {
345 LMBasis Y;
346 PetscInt oldest, next;
347 Vec Fprev_old = NULL;
348 Vec YtFprev_old = NULL;
349 LMProducts YtY = NULL;
350
351 PetscCall(MatLMVMGetRange(B, &oldest, &next));
352 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
353
354 if (cache_YtFprev) {
355 if (!brdn->YtFprev) PetscCall(LMBasisCreateRow(Y, &brdn->YtFprev));
356 PetscCall(LMBasisGetWorkVec(Y, &Fprev_old));
357 PetscCall(MatLMVMGetUpdatedProducts(B, LMBASIS_Y, LMBASIS_Y, LMBLOCK_UPPER_TRIANGLE, &YtY));
358 PetscCall(LMProductsGetNextColumn(YtY, &YtFprev_old));
359 PetscCall(VecCopy(lmvm->Fprev, Fprev_old));
360 PetscCall(VecCopy(brdn->YtFprev, YtFprev_old));
361 }
362 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
363 PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
364 PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
365 /* Accept the update */
366 PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
367 if (cache_YtFprev) {
368 PetscInt oldest_new, next_new;
369
370 PetscCall(MatLMVMGetRange(B, &oldest_new, &next_new));
371 // Compute the one new Y_i^T Fprev_old value
372 PetscCall(LMBasisGEMVH(Y, next, next_new, 1.0, Fprev_old, 0.0, YtFprev_old));
373 PetscCall(LMBasisGEMVH(Y, oldest_new, next_new, 1.0, F, 0.0, brdn->YtFprev));
374 PetscCall(LMBasisSetCachedProduct(Y, F, brdn->YtFprev));
375 PetscCall(VecAXPBY(YtFprev_old, 1.0, -1.0, brdn->YtFprev));
376 PetscCall(LMProductsRestoreNextColumn(YtY, &YtFprev_old));
377 PetscCall(LMBasisRestoreWorkVec(Y, &Fprev_old));
378 }
379 } else if (cache_YtFprev) {
380 LMBasis Y;
381
382 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
383 if (!brdn->YtFprev) PetscCall(LMBasisCreateRow(Y, &brdn->YtFprev));
384 }
385 /* Save the solution and function to be used in the next update */
386 PetscCall(VecCopy(X, lmvm->Xprev));
387 PetscCall(VecCopy(F, lmvm->Fprev));
388 lmvm->prev_set = PETSC_TRUE;
389 PetscFunctionReturn(PETSC_SUCCESS);
390 }
391
MatReset_LMVMBrdn(Mat B,MatLMVMResetMode mode)392 static PetscErrorCode MatReset_LMVMBrdn(Mat B, MatLMVMResetMode mode)
393 {
394 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
395 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
396
397 PetscFunctionBegin;
398 if (MatLMVMResetClearsBases(mode)) {
399 for (PetscInt i = 0; i < BROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisDestroy(&lbrdn->basis[i]));
400 for (PetscInt i = 0; i < BROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsDestroy(&lbrdn->products[i]));
401 PetscCall(VecDestroy(&lbrdn->YtFprev));
402 } else {
403 for (PetscInt i = 0; i < BROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisReset(lbrdn->basis[i]));
404 for (PetscInt i = 0; i < BROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsReset(lbrdn->products[i]));
405 if (lbrdn->YtFprev) PetscCall(VecZeroEntries(lbrdn->YtFprev));
406 }
407 PetscFunctionReturn(PETSC_SUCCESS);
408 }
409
MatDestroy_LMVMBrdn(Mat B)410 static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
411 {
412 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
413
414 PetscFunctionBegin;
415 PetscCall(MatReset_LMVMBrdn(B, MAT_LMVM_RESET_ALL));
416 PetscCall(PetscFree(lmvm->ctx));
417 PetscCall(MatDestroy_LMVM(B));
418 PetscFunctionReturn(PETSC_SUCCESS);
419 }
420
MatLMVMSetMultAlgorithm_Brdn(Mat B)421 static PetscErrorCode MatLMVMSetMultAlgorithm_Brdn(Mat B)
422 {
423 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
424
425 PetscFunctionBegin;
426 switch (lmvm->mult_alg) {
427 case MAT_LMVM_MULT_RECURSIVE:
428 lmvm->ops->mult = MatMult_LMVMBrdn_Recursive;
429 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBrdn_Recursive;
430 lmvm->ops->solve = MatSolve_LMVMBrdn_Recursive;
431 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_Recursive;
432 break;
433 case MAT_LMVM_MULT_DENSE:
434 lmvm->ops->mult = MatMult_LMVMBrdn_Dense;
435 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBrdn_Dense;
436 lmvm->ops->solve = MatSolve_LMVMBrdn_CompactDense;
437 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_CompactDense;
438 break;
439 case MAT_LMVM_MULT_COMPACT_DENSE:
440 lmvm->ops->mult = MatMult_LMVMBrdn_CompactDense;
441 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBrdn_CompactDense;
442 lmvm->ops->solve = MatSolve_LMVMBrdn_CompactDense;
443 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBrdn_CompactDense;
444 break;
445 }
446 PetscFunctionReturn(PETSC_SUCCESS);
447 }
448
MatCreate_LMVMBrdn(Mat B)449 PetscErrorCode MatCreate_LMVMBrdn(Mat B)
450 {
451 Mat_LMVM *lmvm;
452 Mat_Brdn *lbrdn;
453
454 PetscFunctionBegin;
455 PetscCall(MatCreate_LMVM(B));
456 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBROYDEN));
457 B->ops->destroy = MatDestroy_LMVMBrdn;
458
459 lmvm = (Mat_LMVM *)B->data;
460 lmvm->ops->reset = MatReset_LMVMBrdn;
461 lmvm->ops->update = MatUpdate_LMVMBrdn;
462 lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_Brdn;
463
464 PetscCall(PetscNew(&lbrdn));
465 lmvm->ctx = (void *)lbrdn;
466
467 PetscCall(MatLMVMSetMultAlgorithm_Brdn(B));
468 PetscFunctionReturn(PETSC_SUCCESS);
469 }
470
471 /*@
472 MatCreateLMVMBroyden - Creates a limited-memory "good" Broyden-type approximation
473 matrix used for a Jacobian. L-Brdn is not guaranteed to be symmetric or
474 positive-definite.
475
476 To use the L-Brdn matrix with other vector types, the matrix must be
477 created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
478 This ensures that the internal storage and work vectors are duplicated from the
479 correct type of vector.
480
481 Collective
482
483 Input Parameters:
484 + comm - MPI communicator
485 . n - number of local rows for storage vectors
486 - N - global size of the storage vectors
487
488 Output Parameter:
489 . B - the matrix
490
491 Options Database Keys:
492 + -mat_lmvm_hist_size - the number of history vectors to keep
493 . -mat_lmvm_mult_algorithm - the algorithm to use for multiplication (recursive, dense, compact_dense)
494 . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed
495 - -mat_lmvm_debug - (developer) perform internal debugging checks
496
497 Level: intermediate
498
499 Note:
500 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
501 paradigm instead of this routine directly.
502
503 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
504 `MatCreateLMVMBFGS()`, `MatCreateLMVMBadBroyden()`, `MatCreateLMVMSymBroyden()`
505 @*/
MatCreateLMVMBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)506 PetscErrorCode MatCreateLMVMBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
507 {
508 PetscFunctionBegin;
509 PetscCall(KSPInitializePackage());
510 PetscCall(MatCreate(comm, B));
511 PetscCall(MatSetSizes(*B, n, n, N, N));
512 PetscCall(MatSetType(*B, MATLMVMBROYDEN));
513 PetscCall(MatSetUp(*B));
514 PetscFunctionReturn(PETSC_SUCCESS);
515 }
516