Actual source code: lmvmutils.c
1: #include <petscdevice.h>
2: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: #include <petsc/private/deviceimpl.h>
4: #include <petscblaslapack.h>
6: /*@
7: MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to a `MATLMVM` matrix.
9: Input Parameters:
10: + B - A `MATLMVM` matrix
11: . X - Solution vector
12: - F - Function vector
14: Level: intermediate
16: Notes:
18: The first time this function is called for a `MATLMVM` matrix, no update is applied, but the given X and F vectors
19: are stored for use as Xprev and Fprev in the next update.
21: If the user has provided another `MATLMVM` matrix for the reference Jacobian (using `MatLMVMSetJ0()`, for example),
22: that matrix is also updated recursively.
24: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMUpdate()` is
25: called, the row size and layout of `B` will be set to match `F` and the column size and layout of `B` will be set to
26: match `X`, and these sizes will be final.
28: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMAllocate()`
29: @*/
30: PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
31: {
32: Mat_LMVM *lmvm;
33: PetscBool same;
35: PetscFunctionBegin;
39: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
40: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
41: /* If B has specified layouts, this will check X and F are compatible;
42: if B does not have specified layouts, this will adopt them, so that
43: this pattern is okay
45: MatCreate(comm, &B);
46: MatLMVMSetType(B, MATLMVMBFGS);
47: MatLMVMUpdate(B, X, F);
48: */
49: PetscCall(MatLMVMUseVecLayoutsIfCompatible(B, X, F));
50: MatCheckPreallocated(B, 1);
51: PetscCall(PetscLogEventBegin(MATLMVM_Update, NULL, NULL, NULL, NULL));
52: lmvm = (Mat_LMVM *)B->data;
53: PetscCall(MatLMVMUpdate(lmvm->J0, X, F));
54: PetscCall((*lmvm->ops->update)(B, X, F));
55: PetscCall(PetscLogEventEnd(MATLMVM_Update, NULL, NULL, NULL, NULL));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode MatLMVMCreateJ0(Mat B, Mat *J0)
60: {
61: PetscLayout rmap, cmap;
62: VecType vec_type;
63: const char *prefix;
65: PetscFunctionBegin;
66: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), J0));
67: PetscCall(MatGetLayouts(B, &rmap, &cmap));
68: PetscCall(MatSetLayouts(*J0, rmap, cmap));
69: PetscCall(MatGetVecType(B, &vec_type));
70: PetscCall(MatSetVecType(*J0, vec_type));
71: PetscCall(MatGetOptionsPrefix(B, &prefix));
72: PetscCall(MatSetOptionsPrefix(*J0, prefix));
73: PetscCall(MatAppendOptionsPrefix(*J0, "mat_lmvm_J0_"));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode MatLMVMCreateJ0KSP(Mat B, KSP *ksp)
78: {
79: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
80: const char *prefix;
82: PetscFunctionBegin;
83: PetscCall(KSPCreate(PetscObjectComm((PetscObject)B), ksp));
84: PetscCall(KSPSetOperators(*ksp, lmvm->J0, lmvm->J0));
85: PetscCall(PetscObjectIncrementTabLevel((PetscObject)B, (PetscObject)*ksp, 1));
86: PetscCall(MatGetOptionsPrefix(B, &prefix));
87: PetscCall(KSPSetOptionsPrefix(*ksp, prefix));
88: PetscCall(KSPAppendOptionsPrefix(*ksp, "mat_lmvm_J0_"));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode MatLMVMCreateJ0KSP_ExactInverse(Mat B, KSP *ksp)
93: {
94: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
95: PC pc;
97: PetscFunctionBegin;
98: PetscCall(MatLMVMCreateJ0KSP(B, ksp));
99: PetscCall(KSPSetType(*ksp, KSPPREONLY));
100: PetscCall(KSPGetPC(*ksp, &pc));
101: PetscCall(PCSetType(pc, PCMAT));
102: PetscCall(PCMatSetApplyOperation(pc, MATOP_SOLVE));
103: lmvm->disable_ksp_viewers = PETSC_TRUE;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*@
108: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
109: an identity matrix (scale = 1.0).
111: Input Parameter:
112: . B - A `MATLMVM` matrix
114: Level: advanced
116: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
117: @*/
118: PetscErrorCode MatLMVMClearJ0(Mat B)
119: {
120: Mat_LMVM *lmvm;
121: PetscBool same;
123: PetscFunctionBegin;
125: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
126: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
127: lmvm = (Mat_LMVM *)B->data;
128: PetscCall(MatDestroy(&lmvm->J0));
129: PetscCall(KSPDestroy(&lmvm->J0ksp));
130: PetscCall(MatLMVMCreateJ0(B, &lmvm->J0));
131: PetscCall(MatSetType(lmvm->J0, MATCONSTANTDIAGONAL));
132: PetscCall(MatZeroEntries(lmvm->J0));
133: PetscCall(MatShift(lmvm->J0, 1.0));
134: PetscCall(MatLMVMCreateJ0KSP_ExactInverse(B, &lmvm->J0ksp));
135: lmvm->created_J0 = PETSC_TRUE;
136: lmvm->created_J0ksp = PETSC_TRUE;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*@
141: MatLMVMSetJ0Scale - Allows the user to define a scalar value
142: mu such that J0 = mu*I.
144: Input Parameters:
145: + B - A `MATLMVM` matrix
146: - scale - Scalar value mu that defines the initial Jacobian
148: Level: advanced
150: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
151: @*/
152: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
153: {
154: Mat_LMVM *lmvm;
155: PetscBool same;
156: PetscBool isconstant;
158: PetscFunctionBegin;
160: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
161: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
162: lmvm = (Mat_LMVM *)B->data;
163: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATCONSTANTDIAGONAL, &isconstant));
164: if (!isconstant) PetscCall(MatLMVMClearJ0(B));
165: PetscCall(MatZeroEntries(lmvm->J0));
166: PetscCall(MatShift(lmvm->J0, scale));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: #if PetscDefined(USE_DEBUG)
172: do { \
173: if (!(layout)->setupcalled) { \
174: PetscMPIInt global[2]; \
175: global[0] = (PetscMPIInt)(v); \
176: global[1] = -global[0]; \
177: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &global[0], 2, MPI_INT, MPI_MIN, ((layout)->comm))); \
178: PetscCheck(global[1] == -global[0], ((layout)->comm), PETSC_ERR_ARG_WRONGSTATE, "PetscLayout has size == PETSC_DECIDE and local size == PETSC_DETERMINE on only some processes"); \
179: } \
180: } while (0)
181: #else
183: do { \
184: (void)(comm); \
185: (void)(v); \
186: } while (0)
187: #endif
189: static PetscErrorCode MatLMVMCheckArgumentLayout(PetscLayout b, PetscLayout a)
190: {
191: PetscBool b_is_unspecified, a_is_specified, are_compatible;
192: PetscLayout b_setup = NULL, a_setup = NULL;
194: PetscFunctionBegin;
195: if (b == a) PetscFunctionReturn(PETSC_SUCCESS); // a layout is compatible with itself
196: if (b->setupcalled && a->setupcalled) {
197: // run the standard checks that are guaranteed to error on at least one process if the layouts are incompatible
198: PetscCheck(b->N == a->N, b->comm, PETSC_ERR_ARG_SIZ, "argument layout (size %" PetscInt_FMT ") is incompatible with MatLMVM layout (size %" PetscInt_FMT ")", a->N, b->N);
199: PetscCheck(b->n == a->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "argument layout (local size %" PetscInt_FMT ") is incompatible with MatLMVM layout (local size %" PetscInt_FMT ")", a->n, b->n);
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
202: a_is_specified = (a->n >= 0) || (a->N >= 0) ? PETSC_TRUE : PETSC_FALSE;
204: PetscCheck(a_is_specified, a->comm, PETSC_ERR_ARG_WRONGSTATE, "argument layout has n == PETSC_DETERMINE and N == PETSC_DECIDE, size must be specified first");
205: b_is_unspecified = (b->n < 0) && (b->N < 0) ? PETSC_TRUE : PETSC_FALSE;
207: if (b_is_unspecified) PetscFunctionReturn(PETSC_SUCCESS); // any layout can replace an unspecified layout
208: // we don't want to change the setup states in this check, so make duplicates if they have not been setup
209: if (!b->setupcalled) {
210: PetscCall(PetscLayoutDuplicate(b, &b_setup));
211: PetscCall(PetscLayoutSetUp(b_setup));
212: } else PetscCall(PetscLayoutReference(b, &b_setup));
213: if (!a->setupcalled) {
214: PetscCall(PetscLayoutDuplicate(a, &a_setup));
215: PetscCall(PetscLayoutSetUp(a_setup));
216: } else PetscCall(PetscLayoutReference(a, &a_setup));
217: PetscCall(PetscLayoutCompare(b_setup, a_setup, &are_compatible));
218: PetscCall(PetscLayoutDestroy(&a_setup));
219: PetscCall(PetscLayoutDestroy(&b_setup));
220: PetscCheck(are_compatible, b->comm, PETSC_ERR_ARG_SIZ, "argument layout (size %" PetscInt_FMT ") is incompatible with MatLMVM layout (size %" PetscInt_FMT ")", a->N, b->N);
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: static PetscErrorCode MatLMVMUseJ0LayoutsIfCompatible(Mat B, Mat J0)
225: {
226: PetscFunctionBegin;
227: PetscCall(MatLMVMCheckArgumentLayout(B->rmap, J0->rmap));
228: PetscCall(MatLMVMCheckArgumentLayout(B->cmap, J0->cmap));
229: PetscCall(PetscLayoutSetUp(J0->rmap));
230: PetscCall(PetscLayoutSetUp(J0->cmap));
231: PetscCall(PetscLayoutReference(J0->rmap, &B->rmap));
232: PetscCall(PetscLayoutReference(J0->cmap, &B->cmap));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: static PetscErrorCode MatLMVMUseJ0DiagLayoutsIfCompatible(Mat B, Vec J0_diag)
237: {
238: PetscFunctionBegin;
239: PetscCall(MatLMVMCheckArgumentLayout(B->rmap, J0_diag->map));
240: PetscCall(MatLMVMCheckArgumentLayout(B->cmap, J0_diag->map));
241: PetscCall(PetscLayoutSetUp(J0_diag->map));
242: PetscCall(PetscLayoutReference(J0_diag->map, &B->rmap));
243: PetscCall(PetscLayoutReference(J0_diag->map, &B->cmap));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: PETSC_INTERN PetscErrorCode MatLMVMUseVecLayoutsIfCompatible(Mat B, Vec X, Vec F)
248: {
249: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
251: PetscFunctionBegin;
252: PetscCall(MatLMVMCheckArgumentLayout(B->rmap, F->map));
253: PetscCall(MatLMVMCheckArgumentLayout(B->cmap, X->map));
254: PetscCall(PetscLayoutSetUp(F->map));
255: PetscCall(PetscLayoutSetUp(X->map));
256: PetscCall(PetscLayoutReference(F->map, &B->rmap));
257: PetscCall(PetscLayoutReference(X->map, &B->cmap));
258: if (lmvm->created_J0) {
259: PetscCall(PetscLayoutReference(B->rmap, &lmvm->J0->rmap));
260: PetscCall(PetscLayoutReference(B->cmap, &lmvm->J0->cmap));
261: }
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: MatLMVMSetJ0Diag - Allows the user to define a vector
267: V such that J0 = diag(V).
269: Input Parameters:
270: + B - An LMVM-type matrix
271: - V - Vector that defines the diagonal of the initial Jacobian: values are copied, V is not referenced
273: Level: advanced
275: Note:
276: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0Diag()` is
277: called, the rows and columns of `B` will each have the size and layout of `V`, and these sizes will be final.
279: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
280: @*/
281: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
282: {
283: Mat J0diag;
284: PetscBool same;
285: VecType vec_type;
287: PetscFunctionBegin;
290: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
291: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
292: PetscCheckSameComm(B, 1, V, 2);
293: PetscCall(MatLMVMUseJ0DiagLayoutsIfCompatible(B, V));
294: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &J0diag));
295: PetscCall(MatSetLayouts(J0diag, V->map, V->map));
296: PetscCall(VecGetType(V, &vec_type));
297: PetscCall(MatSetVecType(J0diag, vec_type));
298: PetscCall(MatSetType(J0diag, MATDIAGONAL));
299: PetscCall(MatDiagonalSet(J0diag, V, INSERT_VALUES));
300: PetscCall(MatLMVMSetJ0(B, J0diag));
301: PetscCall(MatDestroy(&J0diag));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: PETSC_INTERN PetscErrorCode MatLMVMGetJ0InvDiag(Mat B, Vec *V)
306: {
307: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
308: PetscBool isvdiag;
310: PetscFunctionBegin;
311: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATDIAGONAL, &isvdiag));
312: if (!isvdiag) {
313: PetscCall(MatLMVMClearJ0(B));
314: PetscCall(MatSetType(lmvm->J0, MATDIAGONAL));
315: PetscCall(MatZeroEntries(lmvm->J0));
316: PetscCall(MatShift(lmvm->J0, 1.0));
317: }
318: PetscCall(MatDiagonalGetInverseDiagonal(lmvm->J0, V));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: PETSC_INTERN PetscErrorCode MatLMVMRestoreJ0InvDiag(Mat B, Vec *V)
323: {
324: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
326: PetscFunctionBegin;
327: PetscCall(MatDiagonalRestoreInverseDiagonal(lmvm->J0, V));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: PETSC_INTERN PetscErrorCode MatLMVMJ0KSPIsExact(Mat B, PetscBool *is_exact)
332: {
333: PetscBool is_preonly, is_pcmat, has_pmat;
334: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
335: Mat pc_pmat;
336: PC pc;
337: MatOperation matop;
339: PetscFunctionBegin;
340: *is_exact = PETSC_FALSE;
341: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0ksp, KSPPREONLY, &is_preonly));
342: if (!is_preonly) PetscFunctionReturn(PETSC_SUCCESS);
343: PetscCall(KSPGetPC(lmvm->J0ksp, &pc));
344: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMAT, &is_pcmat));
345: if (!is_pcmat) PetscFunctionReturn(PETSC_SUCCESS);
346: PetscCall(PCGetOperatorsSet(pc, NULL, &has_pmat));
347: if (!has_pmat) PetscFunctionReturn(PETSC_SUCCESS);
348: PetscCall(PCGetOperators(pc, NULL, &pc_pmat));
349: if (pc_pmat != lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
350: PetscCall(PCMatGetApplyOperation(pc, &matop));
351: *is_exact = (matop == MATOP_SOLVE) ? PETSC_TRUE : PETSC_FALSE;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: MatLMVMSetJ0 - Allows the user to define the initial Jacobian matrix from which the LMVM-type approximation is built
357: up.
359: Input Parameters:
360: + B - An LMVM-type matrix
361: - J0 - The initial Jacobian matrix, will be referenced by B.
363: Level: advanced
365: Notes:
366: A KSP is created for inverting J0 with prefix `-mat_lmvm_J0_` and J0 is set to both operators in `KSPSetOperators()`.
367: If you want to use a separate preconditioning matrix, use `MatLMVMSetJ0KSP()` directly.
369: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0()` is
370: called, then `B` will adopt the sizes and layouts of `J0`, and these sizes will be final.
372: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
373: @*/
374: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
375: {
376: Mat_LMVM *lmvm;
377: PetscBool same;
378: PetscBool J0_has_solve;
380: PetscFunctionBegin;
383: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
384: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
385: lmvm = (Mat_LMVM *)B->data;
386: if (J0 == lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
387: PetscCheckSameComm(B, 1, J0, 2);
388: PetscCall(MatLMVMUseJ0LayoutsIfCompatible(B, J0));
389: PetscCall(PetscObjectReference((PetscObject)J0));
390: PetscCall(MatDestroy(&lmvm->J0));
391: lmvm->J0 = J0;
392: lmvm->created_J0 = PETSC_FALSE;
393: PetscCall(MatHasOperation(J0, MATOP_SOLVE, &J0_has_solve));
394: if (J0_has_solve) {
395: PetscCall(KSPDestroy(&lmvm->J0ksp));
396: PetscCall(MatLMVMCreateJ0KSP_ExactInverse(B, &lmvm->J0ksp));
397: lmvm->created_J0ksp = PETSC_TRUE;
398: } else {
399: if (lmvm->created_J0ksp) {
400: PetscBool is_preonly, is_pcmat = PETSC_FALSE, is_pcmat_solve = PETSC_FALSE;
401: PC pc;
403: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0ksp, KSPPREONLY, &is_preonly));
404: PetscCall(KSPGetPC(lmvm->J0ksp, &pc));
405: if (pc) {
406: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMAT, &is_pcmat));
407: if (is_pcmat) {
408: MatOperation matop;
410: PetscCall(PCMatGetApplyOperation(pc, &matop));
411: if (matop == MATOP_SOLVE) is_pcmat_solve = PETSC_TRUE;
412: }
413: }
414: if (is_preonly && is_pcmat_solve) {
415: /* The KSP is one created by LMVM for a mat that has a MatSolve() implementation. Because this new J0 doesn't, change it to
416: a default KSP */
417: PetscCall(KSPDestroy(&lmvm->J0ksp));
418: PetscCall(MatLMVMCreateJ0KSP(B, &lmvm->J0ksp));
419: }
420: }
421: PetscCall(KSPSetOperators(lmvm->J0ksp, J0, J0));
422: }
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@
427: MatLMVMSetJ0PC - Allows the user to define a `PC` object that acts as the initial inverse-Jacobian matrix.
429: Input Parameters:
430: + B - A `MATLMVM` matrix
431: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
433: Level: advanced
435: Notes:
436: `J0pc` should already contain all the operators necessary for its application. The `MATLMVM` matrix only calls
437: `PCApply()` without changing any other options.
439: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0PC()` is
440: called, then `B` will adopt the sizes and layouts of the operators of `J0pc`, and these sizes will be final.
442: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
443: @*/
444: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
445: {
446: Mat_LMVM *lmvm;
447: PetscBool same, mat_set, pmat_set;
448: PC current_pc;
449: Mat J0 = NULL;
451: PetscFunctionBegin;
454: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
455: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
456: lmvm = (Mat_LMVM *)B->data;
457: PetscCall(PCGetOperatorsSet(J0pc, &mat_set, &pmat_set));
458: PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
459: if (mat_set) PetscCall(PCGetOperators(J0pc, &J0, NULL));
460: else PetscCall(PCGetOperators(J0pc, NULL, &J0));
461: PetscCall(KSPGetPC(lmvm->J0ksp, ¤t_pc));
462: if (J0pc == current_pc && J0 == lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
463: PetscCall(MatLMVMSetJ0(B, J0));
464: PetscCall(KSPSetPC(lmvm->J0ksp, J0pc));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured KSP solver for the initial inverse-Jacobian
470: approximation.
472: Input Parameters:
473: + B - A `MATLMVM` matrix
474: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
476: Level: advanced
478: Note:
479: The `KSP` solver should already contain all the operators necessary to perform the inversion. The `MATLMVM` matrix
480: only calls `KSPSolve()` without changing any other options.
482: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0KSP()` is
483: called, then `B` will adopt the sizes and layouts of the operators of `J0ksp`, and these sizes will be final.
485: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
486: @*/
487: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
488: {
489: Mat_LMVM *lmvm;
490: PetscBool same, mat_set, pmat_set;
491: Mat J0;
493: PetscFunctionBegin;
496: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
497: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
498: lmvm = (Mat_LMVM *)B->data;
499: PetscCall(KSPGetOperatorsSet(J0ksp, &mat_set, &pmat_set));
500: PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
501: if (mat_set) PetscCall(KSPGetOperators(J0ksp, &J0, NULL));
502: else PetscCall(KSPGetOperators(J0ksp, NULL, &J0));
503: if (J0ksp == lmvm->J0ksp && lmvm->J0 == J0) PetscFunctionReturn(PETSC_SUCCESS);
504: PetscCall(MatLMVMSetJ0(B, J0));
505: if (J0ksp != lmvm->J0ksp) {
506: lmvm->created_J0ksp = PETSC_FALSE;
507: lmvm->disable_ksp_viewers = PETSC_FALSE; // if the user supplies a more complicated KSP, don't turn off viewers
508: }
509: PetscCall(PetscObjectReference((PetscObject)J0ksp));
510: PetscCall(KSPDestroy(&lmvm->J0ksp));
511: lmvm->J0ksp = J0ksp;
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: /*@
516: MatLMVMGetJ0 - Returns a pointer to the internal `J0` matrix.
518: Input Parameter:
519: . B - A `MATLMVM` matrix
521: Output Parameter:
522: . J0 - `Mat` object for defining the initial Jacobian
524: Level: advanced
526: Note:
528: If `J0` was created by `B` it will have the options prefix `-mat_lmvm_J0_`.
530: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
531: @*/
532: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
533: {
534: Mat_LMVM *lmvm;
535: PetscBool same;
537: PetscFunctionBegin;
539: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
540: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
541: lmvm = (Mat_LMVM *)B->data;
542: *J0 = lmvm->J0;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
548: associated with the initial Jacobian.
550: Input Parameter:
551: . B - A `MATLMVM` matrix
553: Output Parameter:
554: . J0pc - `PC` object for defining the initial inverse-Jacobian
556: Level: advanced
558: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
559: @*/
560: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
561: {
562: Mat_LMVM *lmvm;
563: PetscBool same;
565: PetscFunctionBegin;
567: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
568: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
569: lmvm = (Mat_LMVM *)B->data;
570: PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: /*@
575: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
576: associated with the initial Jacobian.
578: Input Parameter:
579: . B - A `MATLMVM` matrix
581: Output Parameter:
582: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
584: Level: advanced
586: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
587: @*/
588: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
589: {
590: Mat_LMVM *lmvm;
591: PetscBool same;
593: PetscFunctionBegin;
595: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
596: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
597: lmvm = (Mat_LMVM *)B->data;
598: *J0ksp = lmvm->J0ksp;
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*@
603: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
604: matrix-vector product with the initial Jacobian.
606: Input Parameters:
607: + B - A `MATLMVM` matrix
608: - X - vector to multiply with J0
610: Output Parameter:
611: . Y - resulting vector for the operation
613: Level: advanced
615: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
616: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
617: @*/
618: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
619: {
620: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
622: PetscFunctionBegin;
623: PetscCall(MatMult(lmvm->J0, X, Y));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0HermitianTranspose(Mat B, Vec X, Vec Y)
628: {
629: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
631: PetscFunctionBegin;
632: PetscCall(MatMultHermitianTranspose(lmvm->J0, X, Y));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
638: inverse to the given vector.
640: Input Parameters:
641: + B - A `MATLMVM` matrix
642: - X - vector to "multiply" with J0^{-1}
644: Output Parameter:
645: . Y - resulting vector for the operation
647: Level: advanced
649: Note:
650: The specific form of the application
651: depends on whether the user provided a scaling factor, a J0 matrix,
652: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
653: provided, the function simply does an identity matrix application
654: (vector copy).
656: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
657: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
658: @*/
659: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
660: {
661: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
663: PetscFunctionBegin;
664: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
665: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
666: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvTranspose(Mat B, Vec X, Vec Y)
671: {
672: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
674: PetscFunctionBegin;
675: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
676: PetscCall(KSPSolveTranspose(lmvm->J0ksp, X, Y));
677: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvHermitianTranspose(Mat B, Vec X, Vec Y)
682: {
683: PetscFunctionBegin;
684: if (!PetscDefined(USE_COMPLEX)) {
685: PetscCall(MatLMVMApplyJ0InvTranspose(B, X, Y));
686: } else {
687: Vec X_conj;
689: PetscCall(VecDuplicate(X, &X_conj));
690: PetscCall(VecCopy(X, X_conj));
691: PetscCall(VecConjugate(X_conj));
692: PetscCall(MatLMVMApplyJ0InvTranspose(B, X_conj, Y));
693: PetscCall(VecConjugate(Y));
694: PetscCall(VecDestroy(&X_conj));
695: }
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@
700: MatLMVMIsAllocated - Returns a boolean flag that shows whether
701: the necessary data structures for the underlying matrix is allocated.
703: Input Parameter:
704: . B - A `MATLMVM` matrix
706: Output Parameter:
707: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
709: Level: intermediate
711: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
712: @*/
713: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
714: {
715: PetscBool same;
717: PetscFunctionBegin;
719: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
720: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
721: *flg = B->preallocated;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /*@
726: MatLMVMAllocate - Produces all necessary common memory for
727: LMVM approximations based on the solution and function vectors
728: provided.
730: Input Parameters:
731: + B - A `MATLMVM` matrix
732: . X - Solution vector
733: - F - Function vector
735: Level: intermediate
737: Note:
738: If `MatSetSizes()` and `MatSetUp()` have not been called
739: before `MatLMVMAllocate()`, the row layout of `B` will be set to match `F`
740: and the column layout of `B` will be set to match `X`.
742: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
743: @*/
744: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
745: {
746: Mat_LMVM *lmvm;
747: PetscBool same;
749: PetscFunctionBegin;
753: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
754: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
755: lmvm = (Mat_LMVM *)B->data;
756: PetscCall(MatAllocate_LMVM(B, X, F));
757: PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@
762: MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.
764: Input Parameter:
765: . B - A `MATLMVM` matrix
767: Level: intermediate
769: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
770: @*/
771: PetscErrorCode MatLMVMResetShift(Mat B)
772: {
773: Mat_LMVM *lmvm;
774: PetscBool same;
776: PetscFunctionBegin;
778: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
779: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
780: lmvm = (Mat_LMVM *)B->data;
781: lmvm->shift = 0.0;
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: PETSC_INTERN PetscErrorCode MatLMVMReset_Internal(Mat B, MatLMVMResetMode mode)
786: {
787: Mat_LMVM *lmvm;
788: PetscBool same;
790: PetscFunctionBegin;
792: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
793: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
794: lmvm = (Mat_LMVM *)B->data;
795: PetscCall(MatLMVMReset_Internal(lmvm->J0, mode));
796: if (lmvm->ops->reset) PetscCall((*lmvm->ops->reset)(B, mode));
797: PetscCall(MatReset_LMVM(B, mode));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: /*@
802: MatLMVMReset - Flushes all of the accumulated updates out of
803: the `MATLMVM` approximation.
805: Input Parameters:
806: + B - A `MATLMVM` matrix
807: - destructive - flag for enabling destruction of data structures
809: Level: intermediate
811: Note:
812: In practice, this will not actually
813: destroy the data associated with the updates. It simply resets
814: counters, which leads to existing data being overwritten, and
815: `MatSolve()` being applied as if there are no updates. A boolean
816: flag is available to force destruction of the update vectors.
818: If the user has provided another LMVM matrix as J0, the J0
819: matrix is also reset to the identity matrix in this function.
821: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
822: @*/
823: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
824: {
825: Mat_LMVM *lmvm;
826: PetscBool same;
828: PetscFunctionBegin;
830: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
831: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
832: lmvm = (Mat_LMVM *)B->data;
833: PetscCall(PetscInfo(B, "Resetting %s after %" PetscInt_FMT " iterations\n", ((PetscObject)B)->type_name, lmvm->k));
834: PetscCall(MatLMVMReset_Internal(B, destructive ? MAT_LMVM_RESET_ALL : MAT_LMVM_RESET_HISTORY));
835: ++lmvm->nresets;
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: /*@
840: MatLMVMSetHistorySize - Set the number of past iterates to be
841: stored for the construction of the limited-memory quasi-Newton update.
843: Input Parameters:
844: + B - A `MATLMVM` matrix
845: - hist_size - number of past iterates (default 5)
847: Options Database Key:
848: . -mat_lmvm_hist_size m - set number of past iterates
850: Level: beginner
852: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
853: @*/
854: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
855: {
856: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
857: PetscBool same;
859: PetscFunctionBegin;
861: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
862: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
863: PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
864: if (lmvm->m != hist_size) PetscCall(MatLMVMReset_Internal(B, MAT_LMVM_RESET_BASES));
865: lmvm->m = hist_size;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@
870: MatLMVMGetHistorySize - Get the number of past iterates stored for the construction of the limited-memory quasi-Newton update
872: Not Collective
874: Input Parameter:
875: . B - A `MATLMVM` matrix
877: Output Parameter:
878: . hist_size - number of past iterates stored
880: Level: intermediate
882: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetHistorySize()`, `MatLMVMGetUpdateCount()`
883: @*/
884: PetscErrorCode MatLMVMGetHistorySize(Mat B, PetscInt *hist_size)
885: {
886: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
887: PetscBool same;
889: PetscFunctionBegin;
891: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
892: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
893: *hist_size = lmvm->m;
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: /*@
898: MatLMVMGetUpdateCount - Returns the number of accepted updates.
900: Input Parameter:
901: . B - A `MATLMVM` matrix
903: Output Parameter:
904: . nupdates - number of accepted updates
906: Level: intermediate
908: Note:
909: This number may be greater than the total number of update vectors
910: stored in the matrix (`MatLMVMGetHistorySize()`). The counters are reset when `MatLMVMReset()`
911: is called.
913: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
914: @*/
915: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
916: {
917: Mat_LMVM *lmvm;
918: PetscBool same;
920: PetscFunctionBegin;
922: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
923: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
924: lmvm = (Mat_LMVM *)B->data;
925: *nupdates = lmvm->nupdates;
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@
930: MatLMVMGetRejectCount - Returns the number of rejected updates.
931: The counters are reset when `MatLMVMReset()` is called.
933: Input Parameter:
934: . B - A `MATLMVM` matrix
936: Output Parameter:
937: . nrejects - number of rejected updates
939: Level: intermediate
941: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
942: @*/
943: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
944: {
945: Mat_LMVM *lmvm;
946: PetscBool same;
948: PetscFunctionBegin;
950: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
951: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
952: lmvm = (Mat_LMVM *)B->data;
953: *nrejects = lmvm->nrejects;
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: PETSC_INTERN PetscErrorCode MatLMVMGetJ0Scalar(Mat B, PetscBool *is_scalar, PetscScalar *scale)
958: {
959: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
961: PetscFunctionBegin;
962: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATCONSTANTDIAGONAL, is_scalar));
963: if (*is_scalar) PetscCall(MatConstantDiagonalGetConstant(lmvm->J0, scale));
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: static PetscErrorCode MatLMVMUpdateOpVecs(Mat B, LMBasis X, LMBasis OpX, PetscErrorCode (*op)(Mat, Vec, Vec))
968: {
969: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
970: PetscObjectId J0_id;
971: PetscObjectState J0_state;
972: PetscInt oldest, next;
974: PetscFunctionBegin;
975: PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
976: PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
977: PetscCall(LMBasisGetRange(X, &oldest, &next));
978: if (OpX->operator_id != J0_id || OpX->operator_state != J0_state) {
979: // invalidate OpX
980: OpX->k = oldest;
981: OpX->operator_id = J0_id;
982: OpX->operator_state = J0_state;
983: PetscCall(LMBasisSetCachedProduct(OpX, NULL, NULL));
984: }
985: OpX->k = PetscMax(OpX->k, oldest);
986: for (PetscInt i = OpX->k; i < next; i++) {
987: Vec x_i, op_x_i;
989: PetscCall(LMBasisGetVecRead(X, i, &x_i));
990: PetscCall(LMBasisGetNextVec(OpX, &op_x_i));
991: PetscCall(op(B, x_i, op_x_i));
992: PetscCall(LMBasisRestoreNextVec(OpX, &op_x_i));
993: PetscCall(LMBasisRestoreVecRead(X, i, &x_i));
994: }
995: PetscAssert(OpX->k == X->k && OpX->operator_id == J0_id && OpX->operator_state == J0_state, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Invalid state for operator-updated LMBasis");
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: static PetscErrorCode MatLMVMUpdateOpDiffVecs(Mat B, LMBasis Y, PetscScalar alpha, LMBasis OpX, LMBasis YmalphaOpX)
1000: {
1001: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1002: PetscInt start;
1003: PetscObjectId J0_id;
1004: PetscObjectState J0_state;
1005: PetscInt oldest, next;
1007: PetscFunctionBegin;
1008: PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
1009: PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
1010: PetscAssert(Y->m == OpX->m, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Incompatible Y and OpX in MatLMVMUpdateOpDiffVecs()");
1011: PetscAssert(Y->k == OpX->k, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Stale OpX in MatLMVMUpdateOpDiffVecs()");
1012: PetscCall(LMBasisGetRange(Y, &oldest, &next));
1013: if (YmalphaOpX->operator_id != J0_id || YmalphaOpX->operator_state != J0_state) {
1014: // invalidate OpX
1015: YmalphaOpX->k = oldest;
1016: YmalphaOpX->operator_id = J0_id;
1017: YmalphaOpX->operator_state = J0_state;
1018: PetscCall(LMBasisSetCachedProduct(YmalphaOpX, NULL, NULL));
1019: }
1020: YmalphaOpX->k = PetscMax(YmalphaOpX->k, oldest);
1021: start = YmalphaOpX->k;
1022: if (next - start == Y->m) { // full matrix AXPY
1023: PetscCall(MatCopy(Y->vecs, YmalphaOpX->vecs, SAME_NONZERO_PATTERN));
1024: PetscCall(MatAXPY(YmalphaOpX->vecs, -alpha, OpX->vecs, SAME_NONZERO_PATTERN));
1025: YmalphaOpX->k = Y->k;
1026: } else {
1027: for (PetscInt i = start; i < next; i++) {
1028: Vec y_i, op_x_i, y_m_op_x_i;
1030: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
1031: PetscCall(LMBasisGetVecRead(OpX, i, &op_x_i));
1032: PetscCall(LMBasisGetNextVec(YmalphaOpX, &y_m_op_x_i));
1033: PetscCall(VecAXPBYPCZ(y_m_op_x_i, 1.0, -alpha, 0.0, y_i, op_x_i));
1034: PetscCall(LMBasisRestoreNextVec(YmalphaOpX, &y_m_op_x_i));
1035: PetscCall(LMBasisRestoreVecRead(OpX, i, &op_x_i));
1036: PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
1037: }
1038: }
1039: PetscAssert(YmalphaOpX->k == Y->k && YmalphaOpX->operator_id == J0_id && YmalphaOpX->operator_state == J0_state, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Invalid state for operator-updated LMBasis");
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedBasis(Mat B, MatLMVMBasisType type, LMBasis *basis_p, MatLMVMBasisType *returned_type, PetscScalar *scale)
1044: {
1045: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1046: LMBasis basis;
1047: PetscBool is_scalar;
1048: PetscScalar scale_;
1050: PetscFunctionBegin;
1051: switch (type) {
1052: case LMBASIS_S:
1053: case LMBASIS_Y:
1054: *basis_p = lmvm->basis[type];
1055: if (returned_type) *returned_type = type;
1056: if (scale) *scale = 1.0;
1057: break;
1058: case LMBASIS_B0S:
1059: case LMBASIS_H0Y:
1060: // if B_0 = gamma * I, do not actually compute these bases
1061: PetscAssertPointer(returned_type, 4);
1062: PetscAssertPointer(scale, 5);
1063: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1064: if (is_scalar) {
1065: *basis_p = lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y];
1066: *returned_type = (type == LMBASIS_B0S) ? LMBASIS_S : LMBASIS_Y;
1067: *scale = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
1068: } else {
1069: LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];
1071: *returned_type = type;
1072: *scale = 1.0;
1073: if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
1074: basis = lmvm->basis[type];
1075: PetscCall(MatLMVMUpdateOpVecs(B, orig_basis, basis, (type == LMBASIS_B0S) ? MatLMVMApplyJ0Fwd : MatLMVMApplyJ0Inv));
1076: *basis_p = basis;
1077: }
1078: break;
1079: case LMBASIS_S_MINUS_H0Y:
1080: case LMBASIS_Y_MINUS_B0S: {
1081: MatLMVMBasisType op_basis_t = (type == LMBASIS_S_MINUS_H0Y) ? LMBASIS_H0Y : LMBASIS_B0S;
1082: LMBasis op_basis;
1084: if (returned_type) *returned_type = type;
1085: if (scale) *scale = 1.0;
1086: if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
1087: basis = lmvm->basis[type];
1088: PetscCall(MatLMVMGetUpdatedBasis(B, op_basis_t, &op_basis, &op_basis_t, &scale_));
1089: PetscCall(MatLMVMUpdateOpDiffVecs(B, lmvm->basis[MatLMVMBasisSizeOf(type)], scale_, op_basis, basis));
1090: *basis_p = basis;
1091: } break;
1092: default:
1093: PetscUnreachable();
1094: }
1095: basis = *basis_p;
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: PETSC_INTERN PetscErrorCode MatLMVMBasisGetVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
1100: {
1101: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1102: PetscBool is_scalar;
1103: PetscScalar scale_;
1105: PetscFunctionBegin;
1106: switch (type) {
1107: case LMBASIS_B0S:
1108: case LMBASIS_H0Y:
1109: // if B_0 = gamma * I, do not actually compute these bases
1110: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1111: if (is_scalar) {
1112: *scale = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
1113: PetscCall(LMBasisGetVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
1114: } else if (lmvm->do_not_cache_J0_products) {
1115: Vec tmp;
1116: Vec w;
1117: LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];
1118: LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];
1120: PetscCall(LMBasisGetVecRead(orig_basis, i, &w));
1121: PetscCall(LMBasisGetWorkVec(size_basis, &tmp));
1122: if (type == LMBASIS_B0S) PetscCall(MatLMVMApplyJ0Fwd(B, w, tmp));
1123: else PetscCall(MatLMVMApplyJ0Inv(B, w, tmp));
1124: PetscCall(LMBasisRestoreVecRead(orig_basis, i, &w));
1125: *scale = 1.0;
1126: *y = tmp;
1127: } else {
1128: LMBasis basis;
1129: PetscScalar dummy;
1131: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &type, &dummy));
1132: PetscCall(LMBasisGetVecRead(basis, i, y));
1133: *scale = 1.0;
1134: }
1135: break;
1136: default:
1137: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisGetVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y. Use MatLMVMGetUpdatedBasis() and LMBasisGetVecRead().");
1138: }
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1142: PETSC_INTERN PetscErrorCode MatLMVMBasisRestoreVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
1143: {
1144: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1145: PetscBool is_scalar;
1146: PetscScalar scale_;
1148: PetscFunctionBegin;
1149: switch (type) {
1150: case LMBASIS_B0S:
1151: case LMBASIS_H0Y:
1152: // if B_0 = gamma * I, do not actually compute these bases
1153: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1154: if (is_scalar) {
1155: PetscCall(LMBasisRestoreVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
1156: } else if (lmvm->do_not_cache_J0_products) {
1157: LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];
1159: PetscCall(LMBasisRestoreWorkVec(size_basis, y));
1160: } else {
1161: PetscCall(LMBasisRestoreVecRead(lmvm->basis[type], i, y));
1162: }
1163: break;
1164: default:
1165: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisRestoreVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y. Use MatLMVMGetUpdatedBasis() and LMBasisRestoreVecRead().");
1166: }
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: PETSC_INTERN PetscErrorCode MatLMVMGetRange(Mat B, PetscInt *oldest, PetscInt *next)
1171: {
1172: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1174: PetscFunctionBegin;
1175: PetscCall(LMBasisGetRange(lmvm->basis[LMBASIS_S], oldest, next));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: PETSC_INTERN PetscErrorCode MatLMVMGetWorkRow(Mat B, Vec *array_p)
1180: {
1181: LMBasis basis;
1183: PetscFunctionBegin;
1184: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
1185: PetscCall(LMBasisGetWorkRow(basis, array_p));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: PETSC_INTERN PetscErrorCode MatLMVMRestoreWorkRow(Mat B, Vec *array_p)
1190: {
1191: LMBasis basis;
1193: PetscFunctionBegin;
1194: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
1195: PetscCall(LMBasisRestoreWorkRow(basis, array_p));
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscErrorCode MatLMVMApplyOpThenVecs(PetscScalar alpha, Mat B, PetscInt oldest, PetscInt next, MatLMVMBasisType type_S, PetscErrorCode (*op)(Mat, Vec, Vec), Vec x, PetscScalar beta, Vec y)
1200: {
1201: LMBasis S;
1202: Vec B0x;
1204: PetscFunctionBegin;
1205: PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
1206: PetscCall(LMBasisGetWorkVec(S, &B0x));
1207: PetscCall(op(B, x, B0x));
1208: PetscCall(LMBasisGEMVH(S, oldest, next, alpha, B0x, beta, y));
1209: PetscCall(LMBasisRestoreWorkVec(S, &B0x));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: static PetscErrorCode MatLMVMApplyVecsThenOp(PetscScalar alpha, Mat B, PetscInt oldest, PetscInt next, MatLMVMBasisType type_S, MatLMVMBasisType type_Y, PetscErrorCode (*op)(Mat, Vec, Vec), Vec x, PetscScalar beta, Vec y)
1214: {
1215: LMBasis S, Y;
1216: Vec S_x;
1217: Vec B0S_x;
1219: PetscFunctionBegin;
1220: PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
1221: PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, NULL, NULL));
1222: PetscCall(LMBasisGetWorkVec(S, &S_x));
1223: PetscCall(LMBasisGEMV(S, oldest, next, alpha, x, 0.0, S_x));
1224: PetscCall(LMBasisGetWorkVec(Y, &B0S_x));
1225: PetscCall(op(B, S_x, B0S_x));
1226: PetscCall(VecAYPX(y, beta, B0S_x));
1227: PetscCall(LMBasisRestoreWorkVec(Y, &B0S_x));
1228: PetscCall(LMBasisRestoreWorkVec(S, &S_x));
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMVH(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec v, PetscScalar beta, Vec array)
1233: {
1234: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1235: PetscBool cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
1236: LMBasis basis;
1237: MatLMVMBasisType basis_t;
1238: PetscScalar gamma;
1240: PetscFunctionBegin;
1241: if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
1242: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &basis_t, &gamma));
1243: PetscCall(LMBasisGEMVH(basis, oldest, next, alpha * gamma, v, beta, array));
1244: } else {
1245: switch (type) {
1246: case LMBASIS_B0S:
1247: PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_S, MatLMVMApplyJ0HermitianTranspose, v, beta, array));
1248: break;
1249: case LMBASIS_H0Y:
1250: PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_Y, MatLMVMApplyJ0InvHermitianTranspose, v, beta, array));
1251: break;
1252: case LMBASIS_Y_MINUS_B0S:
1253: PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_Y], oldest, next, alpha, v, beta, array));
1254: PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_B0S, oldest, next, -alpha, v, 1.0, array));
1255: break;
1256: case LMBASIS_S_MINUS_H0Y:
1257: PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_S], oldest, next, alpha, v, beta, array));
1258: PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_H0Y, oldest, next, -alpha, v, 1.0, array));
1259: break;
1260: default:
1261: PetscUnreachable();
1262: }
1263: }
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: // x must come from MatLMVMGetRowWork()
1268: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMV(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec x, PetscScalar beta, Vec y)
1269: {
1270: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1271: PetscBool cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
1272: LMBasis basis;
1274: PetscFunctionBegin;
1275: if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
1276: PetscScalar gamma;
1277: MatLMVMBasisType base_type;
1279: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &base_type, &gamma));
1280: PetscCall(LMBasisGEMV(basis, oldest, next, alpha * gamma, x, beta, y));
1281: } else {
1282: switch (type) {
1283: case LMBASIS_B0S:
1284: PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_S, LMBASIS_Y, MatLMVMApplyJ0Fwd, x, beta, y));
1285: break;
1286: case LMBASIS_H0Y:
1287: PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_Y, LMBASIS_S, MatLMVMApplyJ0Inv, x, beta, y));
1288: break;
1289: case LMBASIS_Y_MINUS_B0S:
1290: PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_Y], oldest, next, alpha, x, beta, y));
1291: PetscCall(MatLMVMBasisGEMV(B, LMBASIS_B0S, oldest, next, -alpha, x, 1.0, y));
1292: break;
1293: case LMBASIS_S_MINUS_H0Y:
1294: PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_S], oldest, next, alpha, x, beta, y));
1295: PetscCall(MatLMVMBasisGEMV(B, LMBASIS_H0Y, oldest, next, -alpha, x, 1.0, y));
1296: break;
1297: default:
1298: PetscUnreachable();
1299: }
1300: }
1301: PetscFunctionReturn(PETSC_SUCCESS);
1302: }
1304: PETSC_INTERN PetscErrorCode MatLMVMCreateProducts(Mat B, LMBlockType block_type, LMProducts *products)
1305: {
1306: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1308: PetscFunctionBegin;
1309: PetscCall(LMProductsCreate(lmvm->basis[LMBASIS_S], block_type, products));
1310: (*products)->debug = lmvm->debug;
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: static PetscErrorCode MatLMVMProductsUpdate(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type)
1315: {
1316: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1317: LMBasis X, Y;
1318: MatLMVMBasisType true_type_X, true_type_Y;
1319: PetscScalar alpha_X, alpha_Y;
1320: PetscInt oldest, next;
1321: LMProducts G;
1323: PetscFunctionBegin;
1324: PetscCall(MatLMVMGetUpdatedBasis(B, type_X, &X, &true_type_X, &alpha_X));
1325: PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, &true_type_Y, &alpha_Y));
1326: if (!lmvm->products[block_type][true_type_X][true_type_Y]) PetscCall(MatLMVMCreateProducts(B, block_type, &lmvm->products[block_type][true_type_X][true_type_Y]));
1327: PetscCall(LMProductsUpdate(lmvm->products[block_type][true_type_X][true_type_Y], X, Y));
1328: if (true_type_X == type_X && true_type_Y == type_Y) PetscFunctionReturn(PETSC_SUCCESS);
1329: if (!lmvm->products[block_type][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, block_type, &lmvm->products[block_type][type_X][type_Y]));
1330: G = lmvm->products[block_type][type_X][type_Y];
1331: PetscCall(MatLMVMGetRange(B, &oldest, &next));
1332: PetscCall(LMProductsPrepare(G, lmvm->J0, oldest, next));
1333: if (G->k < lmvm->k) {
1334: PetscCall(LMProductsCopy(lmvm->products[block_type][true_type_X][true_type_Y], lmvm->products[block_type][type_X][type_Y]));
1335: if (alpha_X * alpha_Y != 1.0) PetscCall(LMProductsScale(lmvm->products[block_type][type_X][type_Y], alpha_X * alpha_Y));
1336: }
1337: PetscFunctionReturn(PETSC_SUCCESS);
1338: }
1340: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedProducts(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type, LMProducts *lmwd)
1341: {
1342: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1344: PetscFunctionBegin;
1345: PetscCall(MatLMVMProductsUpdate(B, type_X, type_Y, block_type));
1346: *lmwd = lmvm->products[block_type][type_X][type_Y];
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }
1350: PETSC_INTERN PetscErrorCode MatLMVMProductsInsertDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar v)
1351: {
1352: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1353: LMProducts products;
1355: PetscFunctionBegin;
1356: if (!lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]));
1357: products = lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y];
1358: PetscCall(LMProductsInsertNextDiagonalValue(products, idx, v));
1359: PetscFunctionReturn(PETSC_SUCCESS);
1360: }
1362: PETSC_INTERN PetscErrorCode MatLMVMProductsGetDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar *v)
1363: {
1364: LMProducts products = NULL;
1366: PetscFunctionBegin;
1367: PetscCall(MatLMVMGetUpdatedProducts(B, type_X, type_Y, LMBLOCK_DIAGONAL, &products));
1368: PetscCall(LMProductsGetDiagonalValue(products, idx, v));
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }