Actual source code: brdn.c
1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>
2: #include <petscblaslapack.h>
4: // Broyden is dual to bad Broyden: MatSolve routines are dual to bad Broyden MatMult routines
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: }
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: }
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: }
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: }
34: /* The Broyden kernel can be written as
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
38: This means that we can write B_k x using an intermediate variable alpha_k as
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: */
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;
56: PetscFunctionBegin;
57: PetscCall(MatLMVMGetRange(B, &oldest, &next));
58: if (next > oldest) {
59: LMProducts StS;
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));
65: PetscCall(LMBasisGetWorkVec(S, &G));
66: PetscCall(VecCopy(X, G));
67: PetscCall(VecZeroEntries(BX));
69: for (PetscInt i = next - 1; i >= oldest; i--) {
70: Vec s_i, y_i;
71: PetscScalar sitg, sitsi, alphai;
73: PetscCall(LMBasisGetVecRead(S, i, &s_i));
74: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
75: PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));
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));
82: PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
83: PetscCall(LMBasisRestoreVecRead(S, i, &s_i));
84: }
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: }
97: /* The adjoint of the kernel is
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
101: This means that we can write B_k^H x using an intermediate variable w_k = B_{k-1}^H x
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: */
107: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
108: {
109: PetscInt oldest, next;
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;
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));
124: for (PetscInt i = oldest; i < next; i++) {
125: Vec s_i, y_i;
126: PetscScalar sitBHX, sitsi, yitx;
128: PetscCall(LMBasisGetVecRead(S, i, &s_i));
129: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
130: PetscCall(LMProductsGetDiagonalValue(StS, i, &sitsi));
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: }
144: /*
145: The Broyden kernel can be written as
147: B_k = B_0 + (Y - B_0 S) (triu(S^H S))^{-1} S^H
149: where triu is the upper triangular component. We solve by back substitution each time
150: we apply
151: */
153: PETSC_INTERN PetscErrorCode BroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX)
154: {
155: PetscInt oldest, next;
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;
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: }
178: /*
179: The adjoint of the above formula is
181: B_k^H = B_0^H + S^H (triu(S^H S))^{-H} (Y - B_0 S)^H
182: */
184: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
185: {
186: PetscInt oldest, next;
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;
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: }
209: /*
210: The dense Broyden formula use for CompactDense can be written in a non-update form as
212: B_k = [B_0 | Y] [ I | -S ] [ I ]
213: [---+----] [----------------------]
214: [ 0 | I ] [ triu(S^H S)^{-1} S^H ]
216: The advantage of this form is B_0 appears only once and the (Y - B_0 S) vectors are not needed
217: */
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;
228: PetscFunctionBegin;
229: PetscCall(MatLMVMGetRange(B, &oldest, &next));
230: if (next > oldest) {
231: LMProducts StS;
232: Vec StX;
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: }
256: /*
257: The adoint of the above formula is
259: B_k^h = [I | S triu(S^H S)^{-H} ] [ I | 0 ] [ B_0^H ]
260: [------+---] [-------]
261: [ -S^H | I ] [ Y^H ]
262: */
264: PETSC_INTERN PetscErrorCode BroydenKernelHermitianTranspose_Dense(Mat B, MatLMVMMode mode, Vec X, Vec BHX)
265: {
266: PetscInt oldest, next;
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;
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: }
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: }
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: }
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: }
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: }
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: }
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: }
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;
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;
351: PetscCall(MatLMVMGetRange(B, &oldest, &next));
352: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
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;
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;
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: }
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;
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: }
410: static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
411: {
412: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
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: }
421: static PetscErrorCode MatLMVMSetMultAlgorithm_Brdn(Mat B)
422: {
423: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
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: }
449: PetscErrorCode MatCreate_LMVMBrdn(Mat B)
450: {
451: Mat_LMVM *lmvm;
452: Mat_Brdn *lbrdn;
454: PetscFunctionBegin;
455: PetscCall(MatCreate_LMVM(B));
456: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBROYDEN));
457: B->ops->destroy = MatDestroy_LMVMBrdn;
459: lmvm = (Mat_LMVM *)B->data;
460: lmvm->ops->reset = MatReset_LMVMBrdn;
461: lmvm->ops->update = MatUpdate_LMVMBrdn;
462: lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_Brdn;
464: PetscCall(PetscNew(&lbrdn));
465: lmvm->ctx = (void *)lbrdn;
467: PetscCall(MatLMVMSetMultAlgorithm_Brdn(B));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
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.
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.
481: Collective
483: Input Parameters:
484: + comm - MPI communicator
485: . n - number of local rows for storage vectors
486: - N - global size of the storage vectors
488: Output Parameter:
489: . B - the matrix
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
497: Level: intermediate
499: Note:
500: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
501: paradigm instead of this routine directly.
503: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
504: `MatCreateLMVMBFGS()`, `MatCreateLMVMBadBroyden()`, `MatCreateLMVMSymBroyden()`
505: @*/
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: }