Actual source code: badbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>
3: // Bad Broyden is dual to Broyden: MatSolve routines are dual to Broyden MatMult routines
5: static PetscErrorCode MatSolve_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX)
6: {
7: PetscFunctionBegin;
8: PetscCall(BroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
9: PetscFunctionReturn(PETSC_SUCCESS);
10: }
12: static PetscErrorCode MatSolve_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX)
13: {
14: PetscFunctionBegin;
15: PetscCall(BroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: PETSC_UNUSED static PetscErrorCode MatSolve_LMVMBadBrdn_Dense(Mat B, Vec F, Vec dX)
20: {
21: PetscFunctionBegin;
22: PetscCall(BroydenKernel_Dense(B, MATLMVM_MODE_DUAL, F, dX));
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX)
27: {
28: PetscFunctionBegin;
29: PetscCall(BroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_DUAL, F, dX));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX)
34: {
35: PetscFunctionBegin;
36: PetscCall(BroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_DUAL, F, dX));
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: PETSC_UNUSED static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_Dense(Mat B, Vec F, Vec dX)
41: {
42: PetscFunctionBegin;
43: PetscCall(BroydenKernelHermitianTranspose_Dense(B, MATLMVM_MODE_DUAL, F, dX));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /*
48: The bad Broyden kernel can be written as
50: B_{k+1} x = B_k x + (y_k - B_k s_k)^T * (y_k^T B_k s_k)^{-1} y_k^T B_k x
51: = (I + (y_k - B_k s_k)^T (y_k^T B_k s_k)^{-1} y_k^T) (B_k x)
52: \______________________ _________________________/
53: V
54: recursive rank-1 update
56: When the basis (y_k - B_k s_k) and the products (y_k^T B_k s_k) have been computed, the product can be computed by
57: application of rank-1 updates from oldest to newest
58: */
60: static PetscErrorCode BadBroydenKernel_Recursive_Inner(Mat B, MatLMVMMode mode, PetscInt oldest, PetscInt next, Vec B0X)
61: {
62: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
63: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
64: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
65: LMBasis Y_minus_BkS = lbrdn->basis[LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode)];
66: LMBasis Y;
67: LMProducts YtBkS = lbrdn->products[LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode)];
69: PetscFunctionBegin;
70: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
71: // These cannot be parallelized, notice the data dependence
72: for (PetscInt i = oldest; i < next; i++) {
73: Vec y_i, yimbisi;
74: PetscScalar yitbix;
75: PetscScalar yitbisi;
77: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
78: PetscCall(VecDot(B0X, y_i, &yitbix));
79: PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
80: PetscCall(LMProductsGetDiagonalValue(YtBkS, i, &yitbisi));
81: PetscCall(LMBasisGetVecRead(Y_minus_BkS, i, &yimbisi));
82: PetscCall(VecAXPY(B0X, yitbix / yitbisi, yimbisi));
83: PetscCall(LMBasisRestoreVecRead(Y_minus_BkS, i, &yimbisi));
84: }
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: /*
89: Compute the basis vectors (y_k - B_k s_k) and dot products (y_k^T B_k s_k) recursively
90: */
92: static PetscErrorCode BadBroydenRecursiveBasisUpdate(Mat B, MatLMVMMode mode)
93: {
94: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
95: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
96: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
97: MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
98: LMBasis Y_minus_BkS;
99: LMProducts YtBkS;
100: BroydenBasisType Y_minus_BkS_t = LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode);
101: BroydenProductsType YtBkS_t = LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode);
102: PetscInt oldest, next;
103: PetscInt products_oldest;
104: LMBasis Y;
105: PetscInt start;
107: PetscFunctionBegin;
108: if (!lbrdn->basis[Y_minus_BkS_t]) PetscCall(LMBasisCreate(Y_t == LMBASIS_Y ? lmvm->Fprev : lmvm->Xprev, lmvm->m, &lbrdn->basis[Y_minus_BkS_t]));
109: Y_minus_BkS = lbrdn->basis[Y_minus_BkS_t];
110: if (!lbrdn->products[YtBkS_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lbrdn->products[YtBkS_t]));
111: YtBkS = lbrdn->products[YtBkS_t];
112: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
113: PetscCall(MatLMVMGetRange(B, &oldest, &next));
114: // invalidate computed values if J0 has changed
115: PetscCall(LMProductsPrepare(YtBkS, lmvm->J0, oldest, next));
116: products_oldest = PetscMax(0, YtBkS->k - lmvm->m);
117: if (oldest > products_oldest) {
118: // recursion is starting from a different starting index, it must be recomputed
119: YtBkS->k = oldest;
120: }
121: Y_minus_BkS->k = start = YtBkS->k;
122: // recompute each column vector in Y_minus_BkS, and the product y_k^T B_k s_k, in order
123: for (PetscInt j = start; j < next; j++) {
124: Vec p_j, y_j, B0s_j;
125: PetscScalar yjtbjsj, alpha;
127: PetscCall(LMBasisGetWorkVec(Y_minus_BkS, &p_j));
129: // p_j starts as B_0 * s_j
130: PetscCall(MatLMVMBasisGetVecRead(B, B0S_t, j, &B0s_j, &alpha));
131: PetscCall(VecAXPBY(p_j, alpha, 0.0, B0s_j));
132: PetscCall(MatLMVMBasisRestoreVecRead(B, B0S_t, j, &B0s_j, &alpha));
134: /* Use the matmult kernel to compute p_j = B_j * p_j
135: (if j == oldest, p_j is already correct) */
136: if (j > oldest) PetscCall(BadBroydenKernel_Recursive_Inner(B, mode, oldest, j, p_j));
137: PetscCall(LMBasisGetVecRead(Y, j, &y_j));
138: PetscCall(VecDot(p_j, y_j, &yjtbjsj));
139: PetscCall(LMProductsInsertNextDiagonalValue(YtBkS, j, yjtbjsj));
140: PetscCall(VecAYPX(p_j, -1.0, y_j));
141: PetscCall(LMBasisRestoreVecRead(Y, j, &y_j));
142: PetscCall(LMBasisSetNextVec(Y_minus_BkS, p_j));
143: PetscCall(LMBasisRestoreWorkVec(Y_minus_BkS, &p_j));
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: PETSC_INTERN PetscErrorCode BadBroydenKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec Y)
149: {
150: PetscInt oldest, next;
152: PetscFunctionBegin;
153: PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, Y));
154: PetscCall(MatLMVMGetRange(B, &oldest, &next));
155: if (next > oldest) {
156: PetscCall(BadBroydenRecursiveBasisUpdate(B, mode));
157: PetscCall(BadBroydenKernel_Recursive_Inner(B, mode, oldest, next, Y));
158: }
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*
163: The adjoint of the recursive bad Broyden kernel is
165: B_{k+1}^T x = B_k^T (I + y_k (s_k^T B_k^T y_k)^{-1} (y_k - B_k s_k)^T) x
166: \____________________ ___________________________/
167: V
168: recursive rank-1 update
170: which can be computed by application of rank-1 updates from newest to oldest
171: */
173: static PetscErrorCode BadBroydenKernelHermitianTranspose_Recursive_Inner(Mat B, MatLMVMMode mode, PetscInt oldest, PetscInt next, Vec X)
174: {
175: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
176: BroydenBasisType Y_minus_BkS_t = LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode);
177: BroydenProductsType YtBkS_t = LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode);
178: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
179: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
180: LMBasis Y_minus_BkS = lbrdn->basis[Y_minus_BkS_t];
181: LMBasis Y;
182: LMProducts YtBkS = lbrdn->products[YtBkS_t];
184: PetscFunctionBegin;
185: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
186: // These cannot be parallelized, notice the data dependence
187: for (PetscInt i = next - 1; i >= oldest; i--) {
188: Vec yimbisi, y_i;
189: PetscScalar yimbisitx;
190: PetscScalar yitbisi;
192: PetscCall(LMBasisGetVecRead(Y_minus_BkS, i, &yimbisi));
193: PetscCall(VecDot(X, yimbisi, &yimbisitx));
194: PetscCall(LMBasisRestoreVecRead(Y_minus_BkS, i, &yimbisi));
195: PetscCall(LMProductsGetDiagonalValue(YtBkS, i, &yitbisi));
196: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
197: PetscCall(VecAXPY(X, yimbisitx / PetscConj(yitbisi), y_i));
198: PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
199: }
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: PETSC_INTERN PetscErrorCode BadBroydenKernelHermitianTranspose_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BX)
204: {
205: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
206: PetscInt oldest, next;
207: Vec G = X;
208: LMBasis Y;
210: PetscFunctionBegin;
211: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
212: PetscCall(MatLMVMGetRange(B, &oldest, &next));
213: if (next > oldest) {
214: PetscCall(LMBasisGetWorkVec(Y, &G));
215: PetscCall(VecCopy(X, G));
216: PetscCall(BadBroydenRecursiveBasisUpdate(B, mode));
217: PetscCall(BadBroydenKernelHermitianTranspose_Recursive_Inner(B, mode, oldest, next, G));
218: }
219: PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, G, BX));
220: if (next > oldest) PetscCall(LMBasisRestoreWorkVec(Y, &G));
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*
225: The bad Broyden kernel can be written as
227: B_k = B_0 + (Y - B_0 S) (Y^H B_0 S - stril(Y^H Y))^{-1} Y^H B_0
228: = (I + (Y - B_0 S) (Y^H B_0 S - stril(Y^H Y))^{-1} Y^H) B_0
230: where stril is the strictly lower triangular component. We compute and factorize
231: the small matrix in order to apply a single rank m update
232: */
234: static PetscErrorCode BadBroydenCompactProductsUpdate(Mat B, MatLMVMMode mode)
235: {
236: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
237: MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
238: BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode);
239: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
240: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
241: LMProducts YtB0S, YtY, YtB0S_minus_YtY;
242: LMBasis Y, B0S;
243: PetscScalar alpha;
244: PetscInt oldest, k, next;
245: PetscBool local_is_nonempty;
246: Mat ytbsmyty_local;
248: PetscFunctionBegin;
249: if (!lbrdn->products[YtB0S_minus_YtY_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_FULL, &lbrdn->products[YtB0S_minus_YtY_t]));
250: YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t];
251: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
252: PetscCall(MatLMVMGetUpdatedBasis(B, B0S_t, &B0S, &B0S_t, &alpha));
253: PetscCall(MatLMVMGetRange(B, &oldest, &next));
254: // invalidate computed values if J0 has changed
255: PetscCall(LMProductsPrepare(YtB0S_minus_YtY, lmvm->J0, oldest, next));
256: PetscCall(LMProductsGetLocalMatrix(YtB0S_minus_YtY, &ytbsmyty_local, &k, &local_is_nonempty));
257: if (k < next) {
258: Mat yty_local, ytbs_local;
260: PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, B0S_t, LMBLOCK_FULL, &YtB0S));
261: PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, Y_t, LMBLOCK_STRICT_UPPER_TRIANGLE, &YtY));
262: PetscCall(LMProductsGetLocalMatrix(YtB0S, &ytbs_local, NULL, NULL));
263: PetscCall(LMProductsGetLocalMatrix(YtY, &yty_local, NULL, NULL));
264: if (local_is_nonempty) {
265: PetscCall(MatSetUnfactored(ytbsmyty_local));
266: PetscCall(MatCopy(yty_local, ytbsmyty_local, SAME_NONZERO_PATTERN));
267: PetscCall(MatTranspose(ytbsmyty_local, MAT_INPLACE_MATRIX, &ytbsmyty_local));
268: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ytbsmyty_local));
269: PetscCall(MatScale(ytbsmyty_local, -1.0));
270: PetscCall(MatAXPY(ytbsmyty_local, alpha, ytbs_local, SAME_NONZERO_PATTERN));
271: PetscCall(LMProductsOnesOnUnusedDiagonal(ytbsmyty_local, oldest, next));
272: PetscCall(MatLUFactor(ytbsmyty_local, NULL, NULL, NULL));
273: }
274: PetscCall(LMProductsRestoreLocalMatrix(YtY, &yty_local, NULL));
275: PetscCall(LMProductsRestoreLocalMatrix(YtB0S, &ytbs_local, NULL));
276: }
277: PetscCall(LMProductsRestoreLocalMatrix(YtB0S_minus_YtY, &ytbsmyty_local, &next));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: PETSC_INTERN PetscErrorCode BadBroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
282: {
283: PetscInt oldest, next;
285: PetscFunctionBegin;
286: PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX));
287: PetscCall(MatLMVMGetRange(B, &oldest, &next));
288: if (next > oldest) {
289: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
290: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
291: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
292: MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
293: BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode);
294: LMProducts YtB0S_minus_YtY;
295: LMBasis Y;
296: Vec YtB0X, v;
298: PetscCall(BadBroydenCompactProductsUpdate(B, mode));
299: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
300: YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t];
301: PetscCall(MatLMVMGetWorkRow(B, &YtB0X));
302: PetscCall(MatLMVMGetWorkRow(B, &v));
303: PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, BX, 0.0, YtB0X));
304: PetscCall(LMProductsSolve(YtB0S_minus_YtY, oldest, next, YtB0X, v, /* ^H */ PETSC_FALSE));
305: PetscCall(MatLMVMBasisGEMV(B, Y_minus_B0S_t, oldest, next, 1.0, v, 1.0, BX));
306: PetscCall(MatLMVMRestoreWorkRow(B, &v));
307: PetscCall(MatLMVMRestoreWorkRow(B, &YtB0X));
308: }
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*
313: The adjoint of the above formula for the bad Broyden kernel is
315: B_k^H = B_0^H + B_0^H Y (Y^H B_0 S - stril(Y^H Y))^{-H} (Y - B_0 S)^H
316: = B_0^H (I + Y (Y^H B_0 S - stril(Y^H Y))^{-H} (Y - B_0 S)^H)
317: */
319: PETSC_INTERN PetscErrorCode BadBroydenKernelHermitianTranspose_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
320: {
321: MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
322: PetscInt oldest, next;
323: Vec G = X;
324: LMBasis Y = NULL;
326: PetscFunctionBegin;
327: PetscCall(MatLMVMGetRange(B, &oldest, &next));
328: if (next > oldest) {
329: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
330: Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx;
331: MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode);
332: BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode);
333: LMProducts YtB0S_minus_YtY;
334: Vec YmB0StG, v;
336: PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
337: PetscCall(LMBasisGetWorkVec(Y, &G));
338: PetscCall(VecCopy(X, G));
339: PetscCall(BadBroydenCompactProductsUpdate(B, mode));
340: YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t];
341: PetscCall(MatLMVMGetWorkRow(B, &YmB0StG));
342: PetscCall(MatLMVMGetWorkRow(B, &v));
343: PetscCall(MatLMVMBasisGEMVH(B, Y_minus_B0S_t, oldest, next, 1.0, G, 0.0, YmB0StG));
344: PetscCall(LMProductsSolve(YtB0S_minus_YtY, oldest, next, YmB0StG, v, /* ^H */ PETSC_TRUE));
345: PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, v, 1.0, G));
346: PetscCall(MatLMVMRestoreWorkRow(B, &v));
347: PetscCall(MatLMVMRestoreWorkRow(B, &YmB0StG));
348: }
349: PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, G, BHX));
350: if (next > oldest) PetscCall(LMBasisRestoreWorkVec(Y, &G));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode MatMult_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX)
355: {
356: PetscFunctionBegin;
357: PetscCall(BadBroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, F, dX));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: static PetscErrorCode MatMult_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX)
362: {
363: PetscFunctionBegin;
364: PetscCall(BadBroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, F, dX));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode MatMultHermitianTranspose_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX)
369: {
370: PetscFunctionBegin;
371: PetscCall(BadBroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_PRIMAL, F, dX));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode MatMultHermitianTranspose_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX)
376: {
377: PetscFunctionBegin;
378: PetscCall(BadBroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_PRIMAL, F, dX));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode MatLMVMSetMultAlgorithm_BadBrdn(Mat B)
383: {
384: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
386: PetscFunctionBegin;
387: switch (lmvm->mult_alg) {
388: case MAT_LMVM_MULT_RECURSIVE:
389: lmvm->ops->mult = MatMult_LMVMBadBrdn_Recursive;
390: lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_Recursive;
391: lmvm->ops->solve = MatSolve_LMVMBadBrdn_Recursive;
392: lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_Recursive;
393: break;
394: case MAT_LMVM_MULT_DENSE:
395: lmvm->ops->mult = MatMult_LMVMBadBrdn_CompactDense;
396: lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_CompactDense;
397: lmvm->ops->solve = MatSolve_LMVMBadBrdn_Dense;
398: lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_Dense;
399: break;
400: case MAT_LMVM_MULT_COMPACT_DENSE:
401: lmvm->ops->mult = MatMult_LMVMBadBrdn_CompactDense;
402: lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_CompactDense;
403: lmvm->ops->solve = MatSolve_LMVMBadBrdn_CompactDense;
404: lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_CompactDense;
405: break;
406: }
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: PetscErrorCode MatCreate_LMVMBadBrdn(Mat B)
411: {
412: Mat_LMVM *lmvm;
414: PetscFunctionBegin;
415: PetscCall(MatCreate_LMVMBrdn(B));
416: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBROYDEN));
417: lmvm = (Mat_LMVM *)B->data;
418: lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_BadBrdn;
419: lmvm->cache_gradient_products = PETSC_TRUE;
420: PetscCall(MatLMVMSetMultAlgorithm_BadBrdn(B));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type
426: approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be
427: symmetric or positive-definite.
429: To use the L-BadBrdn matrix with other vector types, the matrix must be
430: created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
431: This ensures that the internal storage and work vectors are duplicated from the
432: correct type of vector.
434: Collective
436: Input Parameters:
437: + comm - MPI communicator
438: . n - number of local rows for storage vectors
439: - N - global size of the storage vectors
441: Output Parameter:
442: . B - the matrix
444: Options Database Keys:
445: + -mat_lmvm_hist_size - the number of history vectors to keep
446: . -mat_lmvm_mult_algorithm - the algorithm to use for multiplication (recursive, dense, compact_dense)
447: . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed
448: - -mat_lmvm_debug - (developer) perform internal debugging checks
450: Level: intermediate
452: Note:
453: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
454: paradigm instead of this routine directly.
456: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBADBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
457: `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMSymBroyden()`
458: @*/
459: PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
460: {
461: PetscFunctionBegin;
462: PetscCall(KSPInitializePackage());
463: PetscCall(MatCreate(comm, B));
464: PetscCall(MatSetSizes(*B, n, n, N, N));
465: PetscCall(MatSetType(*B, MATLMVMBADBROYDEN));
466: PetscCall(MatSetUp(*B));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }