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: PetscCall(KSPSetOperators(lmvm->J0ksp, J0, J0));
400: }
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: MatLMVMSetJ0PC - Allows the user to define a `PC` object that acts as the initial inverse-Jacobian matrix.
407: Input Parameters:
408: + B - A `MATLMVM` matrix
409: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
411: Level: advanced
413: Notes:
414: `J0pc` should already contain all the operators necessary for its application. The `MATLMVM` matrix only calls
415: `PCApply()` without changing any other options.
417: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0PC()` is
418: called, then `B` will adopt the sizes and layouts of the operators of `J0pc`, and these sizes will be final.
420: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
421: @*/
422: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
423: {
424: Mat_LMVM *lmvm;
425: PetscBool same, mat_set, pmat_set;
426: PC current_pc;
427: Mat J0 = NULL;
429: PetscFunctionBegin;
432: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
433: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
434: lmvm = (Mat_LMVM *)B->data;
435: PetscCall(PCGetOperatorsSet(J0pc, &mat_set, &pmat_set));
436: PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
437: if (mat_set) PetscCall(PCGetOperators(J0pc, &J0, NULL));
438: else PetscCall(PCGetOperators(J0pc, NULL, &J0));
439: PetscCall(KSPGetPC(lmvm->J0ksp, ¤t_pc));
440: if (J0pc == current_pc && J0 == lmvm->J0) PetscFunctionReturn(PETSC_SUCCESS);
441: PetscCall(MatLMVMSetJ0(B, J0));
442: PetscCall(KSPSetPC(lmvm->J0ksp, J0pc));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured KSP solver for the initial inverse-Jacobian
448: approximation.
450: Input Parameters:
451: + B - A `MATLMVM` matrix
452: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
454: Level: advanced
456: Note:
457: The `KSP` solver should already contain all the operators necessary to perform the inversion. The `MATLMVM` matrix
458: only calls `KSPSolve()` without changing any other options.
460: If the sizes of `B` have not been specified (using `MatSetSizes()` or `MatSetLayouts()`) before `MatLMVMSetJ0KSP()` is
461: called, then `B` will adopt the sizes and layouts of the operators of `J0ksp`, and these sizes will be final.
463: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
464: @*/
465: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
466: {
467: Mat_LMVM *lmvm;
468: PetscBool same, mat_set, pmat_set;
469: Mat J0;
471: PetscFunctionBegin;
474: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
475: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
476: lmvm = (Mat_LMVM *)B->data;
477: PetscCall(KSPGetOperatorsSet(J0ksp, &mat_set, &pmat_set));
478: PetscCheck(mat_set || pmat_set, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "PC has not operators, call PCSetOperators() before MatLMVMSetJ0PC()");
479: if (mat_set) PetscCall(KSPGetOperators(J0ksp, &J0, NULL));
480: else PetscCall(KSPGetOperators(J0ksp, NULL, &J0));
481: if (J0ksp == lmvm->J0ksp && lmvm->J0 == J0) PetscFunctionReturn(PETSC_SUCCESS);
482: PetscCall(MatLMVMSetJ0(B, J0));
483: if (J0ksp != lmvm->J0ksp) {
484: lmvm->created_J0ksp = PETSC_FALSE;
485: lmvm->disable_ksp_viewers = PETSC_FALSE; // if the user supplies a more complicated KSP, don't turn off viewers
486: }
487: PetscCall(PetscObjectReference((PetscObject)J0ksp));
488: PetscCall(KSPDestroy(&lmvm->J0ksp));
489: lmvm->J0ksp = J0ksp;
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@
494: MatLMVMGetJ0 - Returns a pointer to the internal `J0` matrix.
496: Input Parameter:
497: . B - A `MATLMVM` matrix
499: Output Parameter:
500: . J0 - `Mat` object for defining the initial Jacobian
502: Level: advanced
504: Note:
506: If `J0` was created by `B` it will have the options prefix `-mat_lmvm_J0_`.
508: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
509: @*/
510: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
511: {
512: Mat_LMVM *lmvm;
513: PetscBool same;
515: PetscFunctionBegin;
517: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
518: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
519: lmvm = (Mat_LMVM *)B->data;
520: *J0 = lmvm->J0;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
526: associated with the initial Jacobian.
528: Input Parameter:
529: . B - A `MATLMVM` matrix
531: Output Parameter:
532: . J0pc - `PC` object for defining the initial inverse-Jacobian
534: Level: advanced
536: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
537: @*/
538: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
539: {
540: Mat_LMVM *lmvm;
541: PetscBool same;
543: PetscFunctionBegin;
545: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
546: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
547: lmvm = (Mat_LMVM *)B->data;
548: PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /*@
553: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
554: associated with the initial Jacobian.
556: Input Parameter:
557: . B - A `MATLMVM` matrix
559: Output Parameter:
560: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
562: Level: advanced
564: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
565: @*/
566: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
567: {
568: Mat_LMVM *lmvm;
569: PetscBool same;
571: PetscFunctionBegin;
573: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
574: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
575: lmvm = (Mat_LMVM *)B->data;
576: *J0ksp = lmvm->J0ksp;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
582: matrix-vector product with the initial Jacobian.
584: Input Parameters:
585: + B - A `MATLMVM` matrix
586: - X - vector to multiply with J0
588: Output Parameter:
589: . Y - resulting vector for the operation
591: Level: advanced
593: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
594: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
595: @*/
596: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
597: {
598: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
600: PetscFunctionBegin;
601: PetscCall(MatMult(lmvm->J0, X, Y));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0HermitianTranspose(Mat B, Vec X, Vec Y)
606: {
607: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
609: PetscFunctionBegin;
610: PetscCall(MatMultHermitianTranspose(lmvm->J0, X, Y));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*@
615: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
616: inverse to the given vector.
618: Input Parameters:
619: + B - A `MATLMVM` matrix
620: - X - vector to "multiply" with J0^{-1}
622: Output Parameter:
623: . Y - resulting vector for the operation
625: Level: advanced
627: Note:
628: The specific form of the application
629: depends on whether the user provided a scaling factor, a J0 matrix,
630: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
631: provided, the function simply does an identity matrix application
632: (vector copy).
634: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
635: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
636: @*/
637: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
638: {
639: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
641: PetscFunctionBegin;
642: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
643: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
644: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvTranspose(Mat B, Vec X, Vec Y)
649: {
650: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
652: PetscFunctionBegin;
653: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPushCreateViewerOff(PETSC_TRUE));
654: PetscCall(KSPSolveTranspose(lmvm->J0ksp, X, Y));
655: if (lmvm->disable_ksp_viewers) PetscCall(PetscOptionsPopCreateViewerOff());
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: PETSC_INTERN PetscErrorCode MatLMVMApplyJ0InvHermitianTranspose(Mat B, Vec X, Vec Y)
660: {
661: PetscFunctionBegin;
662: if (!PetscDefined(USE_COMPLEX)) {
663: PetscCall(MatLMVMApplyJ0InvTranspose(B, X, Y));
664: } else {
665: Vec X_conj;
667: PetscCall(VecDuplicate(X, &X_conj));
668: PetscCall(VecCopy(X, X_conj));
669: PetscCall(VecConjugate(X_conj));
670: PetscCall(MatLMVMApplyJ0InvTranspose(B, X_conj, Y));
671: PetscCall(VecConjugate(Y));
672: PetscCall(VecDestroy(&X_conj));
673: }
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: /*@
678: MatLMVMIsAllocated - Returns a boolean flag that shows whether
679: the necessary data structures for the underlying matrix is allocated.
681: Input Parameter:
682: . B - A `MATLMVM` matrix
684: Output Parameter:
685: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
687: Level: intermediate
689: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
690: @*/
691: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
692: {
693: PetscBool same;
695: PetscFunctionBegin;
697: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
698: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
699: *flg = B->preallocated;
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: /*@
704: MatLMVMAllocate - Produces all necessary common memory for
705: LMVM approximations based on the solution and function vectors
706: provided.
708: Input Parameters:
709: + B - A `MATLMVM` matrix
710: . X - Solution vector
711: - F - Function vector
713: Level: intermediate
715: Note:
716: If `MatSetSizes()` and `MatSetUp()` have not been called
717: before `MatLMVMAllocate()`, the row layout of `B` will be set to match `F`
718: and the column layout of `B` will be set to match `X`.
720: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
721: @*/
722: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
723: {
724: Mat_LMVM *lmvm;
725: PetscBool same;
727: PetscFunctionBegin;
731: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
732: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
733: lmvm = (Mat_LMVM *)B->data;
734: PetscCall(MatAllocate_LMVM(B, X, F));
735: PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: /*@
740: MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.
742: Input Parameter:
743: . B - A `MATLMVM` matrix
745: Level: intermediate
747: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
748: @*/
749: PetscErrorCode MatLMVMResetShift(Mat B)
750: {
751: Mat_LMVM *lmvm;
752: PetscBool same;
754: PetscFunctionBegin;
756: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
757: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
758: lmvm = (Mat_LMVM *)B->data;
759: lmvm->shift = 0.0;
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: PETSC_INTERN PetscErrorCode MatLMVMReset_Internal(Mat B, MatLMVMResetMode mode)
764: {
765: Mat_LMVM *lmvm;
766: PetscBool same;
768: PetscFunctionBegin;
770: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
771: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
772: lmvm = (Mat_LMVM *)B->data;
773: PetscCall(MatLMVMReset_Internal(lmvm->J0, mode));
774: if (lmvm->ops->reset) PetscCall((*lmvm->ops->reset)(B, mode));
775: PetscCall(MatReset_LMVM(B, mode));
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: /*@
780: MatLMVMReset - Flushes all of the accumulated updates out of
781: the `MATLMVM` approximation.
783: Input Parameters:
784: + B - A `MATLMVM` matrix
785: - destructive - flag for enabling destruction of data structures
787: Level: intermediate
789: Note:
790: In practice, this will not actually
791: destroy the data associated with the updates. It simply resets
792: counters, which leads to existing data being overwritten, and
793: `MatSolve()` being applied as if there are no updates. A boolean
794: flag is available to force destruction of the update vectors.
796: If the user has provided another LMVM matrix as J0, the J0
797: matrix is also reset to the identity matrix in this function.
799: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
800: @*/
801: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
802: {
803: Mat_LMVM *lmvm;
804: PetscBool same;
806: PetscFunctionBegin;
808: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
809: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
810: lmvm = (Mat_LMVM *)B->data;
811: PetscCall(PetscInfo(B, "Reseting %s after %" PetscInt_FMT " iterations\n", ((PetscObject)B)->type_name, lmvm->k));
812: PetscCall(MatLMVMReset_Internal(B, destructive ? MAT_LMVM_RESET_ALL : MAT_LMVM_RESET_HISTORY));
813: ++lmvm->nresets;
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /*@
818: MatLMVMSetHistorySize - Set the number of past iterates to be
819: stored for the construction of the limited-memory quasi-Newton update.
821: Input Parameters:
822: + B - A `MATLMVM` matrix
823: - hist_size - number of past iterates (default 5)
825: Options Database Key:
826: . -mat_lmvm_hist_size <m> - set number of past iterates
828: Level: beginner
830: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
831: @*/
832: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
833: {
834: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
835: PetscBool same;
837: PetscFunctionBegin;
839: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
840: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
841: PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
842: if (lmvm->m != hist_size) PetscCall(MatLMVMReset_Internal(B, MAT_LMVM_RESET_BASES));
843: lmvm->m = hist_size;
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: PetscErrorCode MatLMVMGetHistorySize(Mat B, PetscInt *hist_size)
848: {
849: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
850: PetscBool same;
852: PetscFunctionBegin;
854: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
855: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
856: *hist_size = lmvm->m;
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /*@
861: MatLMVMGetUpdateCount - Returns the number of accepted updates.
863: Input Parameter:
864: . B - A `MATLMVM` matrix
866: Output Parameter:
867: . nupdates - number of accepted updates
869: Level: intermediate
871: Note:
872: This number may be greater than the total number of update vectors
873: stored in the matrix (`MatLMVMGetHistorySize()`). The counters are reset when `MatLMVMReset()`
874: is called.
876: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
877: @*/
878: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
879: {
880: Mat_LMVM *lmvm;
881: PetscBool same;
883: PetscFunctionBegin;
885: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
886: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
887: lmvm = (Mat_LMVM *)B->data;
888: *nupdates = lmvm->nupdates;
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: /*@
893: MatLMVMGetRejectCount - Returns the number of rejected updates.
894: The counters are reset when `MatLMVMReset()` is called.
896: Input Parameter:
897: . B - A `MATLMVM` matrix
899: Output Parameter:
900: . nrejects - number of rejected updates
902: Level: intermediate
904: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
905: @*/
906: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
907: {
908: Mat_LMVM *lmvm;
909: PetscBool same;
911: PetscFunctionBegin;
913: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
914: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
915: lmvm = (Mat_LMVM *)B->data;
916: *nrejects = lmvm->nrejects;
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: PETSC_INTERN PetscErrorCode MatLMVMGetJ0Scalar(Mat B, PetscBool *is_scalar, PetscScalar *scale)
921: {
922: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
924: PetscFunctionBegin;
925: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->J0, MATCONSTANTDIAGONAL, is_scalar));
926: if (*is_scalar) PetscCall(MatConstantDiagonalGetConstant(lmvm->J0, scale));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode MatLMVMUpdateOpVecs(Mat B, LMBasis X, LMBasis OpX, PetscErrorCode (*op)(Mat, Vec, Vec))
931: {
932: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
933: PetscObjectId J0_id;
934: PetscObjectState J0_state;
935: PetscInt oldest, next;
937: PetscFunctionBegin;
938: PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
939: PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
940: PetscCall(LMBasisGetRange(X, &oldest, &next));
941: if (OpX->operator_id != J0_id || OpX->operator_state != J0_state) {
942: // invalidate OpX
943: OpX->k = oldest;
944: OpX->operator_id = J0_id;
945: OpX->operator_state = J0_state;
946: PetscCall(LMBasisSetCachedProduct(OpX, NULL, NULL));
947: }
948: OpX->k = PetscMax(OpX->k, oldest);
949: for (PetscInt i = OpX->k; i < next; i++) {
950: Vec x_i, op_x_i;
952: PetscCall(LMBasisGetVecRead(X, i, &x_i));
953: PetscCall(LMBasisGetNextVec(OpX, &op_x_i));
954: PetscCall(op(B, x_i, op_x_i));
955: PetscCall(LMBasisRestoreNextVec(OpX, &op_x_i));
956: PetscCall(LMBasisRestoreVecRead(X, i, &x_i));
957: }
958: 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");
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: static PetscErrorCode MatLMVMUpdateOpDiffVecs(Mat B, LMBasis Y, PetscScalar alpha, LMBasis OpX, LMBasis YmalphaOpX)
963: {
964: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
965: PetscInt start;
966: PetscObjectId J0_id;
967: PetscObjectState J0_state;
968: PetscInt oldest, next;
970: PetscFunctionBegin;
971: PetscCall(PetscObjectGetId((PetscObject)lmvm->J0, &J0_id));
972: PetscCall(PetscObjectStateGet((PetscObject)lmvm->J0, &J0_state));
973: PetscAssert(Y->m == OpX->m, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Incompatible Y and OpX in MatLMVMUpdateOpDiffVecs()");
974: PetscAssert(Y->k == OpX->k, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "Stale OpX in MatLMVMUpdateOpDiffVecs()");
975: PetscCall(LMBasisGetRange(Y, &oldest, &next));
976: if (YmalphaOpX->operator_id != J0_id || YmalphaOpX->operator_state != J0_state) {
977: // invalidate OpX
978: YmalphaOpX->k = oldest;
979: YmalphaOpX->operator_id = J0_id;
980: YmalphaOpX->operator_state = J0_state;
981: PetscCall(LMBasisSetCachedProduct(YmalphaOpX, NULL, NULL));
982: }
983: YmalphaOpX->k = PetscMax(YmalphaOpX->k, oldest);
984: start = YmalphaOpX->k;
985: if (next - start == Y->m) { // full matrix AXPY
986: PetscCall(MatCopy(Y->vecs, YmalphaOpX->vecs, SAME_NONZERO_PATTERN));
987: PetscCall(MatAXPY(YmalphaOpX->vecs, -alpha, OpX->vecs, SAME_NONZERO_PATTERN));
988: YmalphaOpX->k = Y->k;
989: } else {
990: for (PetscInt i = start; i < next; i++) {
991: Vec y_i, op_x_i, y_m_op_x_i;
993: PetscCall(LMBasisGetVecRead(Y, i, &y_i));
994: PetscCall(LMBasisGetVecRead(OpX, i, &op_x_i));
995: PetscCall(LMBasisGetNextVec(YmalphaOpX, &y_m_op_x_i));
996: PetscCall(VecAXPBYPCZ(y_m_op_x_i, 1.0, -alpha, 0.0, y_i, op_x_i));
997: PetscCall(LMBasisRestoreNextVec(YmalphaOpX, &y_m_op_x_i));
998: PetscCall(LMBasisRestoreVecRead(OpX, i, &op_x_i));
999: PetscCall(LMBasisRestoreVecRead(Y, i, &y_i));
1000: }
1001: }
1002: 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");
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedBasis(Mat B, MatLMVMBasisType type, LMBasis *basis_p, MatLMVMBasisType *returned_type, PetscScalar *scale)
1007: {
1008: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1009: LMBasis basis;
1010: PetscBool is_scalar;
1011: PetscScalar scale_;
1013: PetscFunctionBegin;
1014: switch (type) {
1015: case LMBASIS_S:
1016: case LMBASIS_Y:
1017: *basis_p = lmvm->basis[type];
1018: if (returned_type) *returned_type = type;
1019: if (scale) *scale = 1.0;
1020: break;
1021: case LMBASIS_B0S:
1022: case LMBASIS_H0Y:
1023: // if B_0 = gamma * I, do not actually compute these bases
1024: PetscAssertPointer(returned_type, 4);
1025: PetscAssertPointer(scale, 5);
1026: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1027: if (is_scalar) {
1028: *basis_p = lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y];
1029: *returned_type = (type == LMBASIS_B0S) ? LMBASIS_S : LMBASIS_Y;
1030: *scale = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
1031: } else {
1032: LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];
1034: *returned_type = type;
1035: *scale = 1.0;
1036: if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
1037: basis = lmvm->basis[type];
1038: PetscCall(MatLMVMUpdateOpVecs(B, orig_basis, basis, (type == LMBASIS_B0S) ? MatLMVMApplyJ0Fwd : MatLMVMApplyJ0Inv));
1039: *basis_p = basis;
1040: }
1041: break;
1042: case LMBASIS_S_MINUS_H0Y:
1043: case LMBASIS_Y_MINUS_B0S: {
1044: MatLMVMBasisType op_basis_t = (type == LMBASIS_S_MINUS_H0Y) ? LMBASIS_H0Y : LMBASIS_B0S;
1045: LMBasis op_basis;
1047: if (returned_type) *returned_type = type;
1048: if (scale) *scale = 1.0;
1049: if (!lmvm->basis[type]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(type) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lmvm->basis[type]));
1050: basis = lmvm->basis[type];
1051: PetscCall(MatLMVMGetUpdatedBasis(B, op_basis_t, &op_basis, &op_basis_t, &scale_));
1052: PetscCall(MatLMVMUpdateOpDiffVecs(B, lmvm->basis[MatLMVMBasisSizeOf(type)], scale_, op_basis, basis));
1053: *basis_p = basis;
1054: } break;
1055: default:
1056: PetscUnreachable();
1057: }
1058: basis = *basis_p;
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: PETSC_INTERN PetscErrorCode MatLMVMBasisGetVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
1063: {
1064: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1065: PetscBool is_scalar;
1066: PetscScalar scale_;
1068: PetscFunctionBegin;
1069: switch (type) {
1070: case LMBASIS_B0S:
1071: case LMBASIS_H0Y:
1072: // if B_0 = gamma * I, do not actually compute these bases
1073: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1074: if (is_scalar) {
1075: *scale = (type == LMBASIS_B0S) ? scale_ : (1.0 / scale_);
1076: PetscCall(LMBasisGetVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
1077: } else if (lmvm->do_not_cache_J0_products) {
1078: Vec tmp;
1079: Vec w;
1080: LMBasis orig_basis = (type == LMBASIS_B0S) ? lmvm->basis[LMBASIS_S] : lmvm->basis[LMBASIS_Y];
1081: LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];
1083: PetscCall(LMBasisGetVecRead(orig_basis, i, &w));
1084: PetscCall(LMBasisGetWorkVec(size_basis, &tmp));
1085: if (type == LMBASIS_B0S) {
1086: PetscCall(MatLMVMApplyJ0Fwd(B, w, tmp));
1087: } else {
1088: PetscCall(MatLMVMApplyJ0Inv(B, w, tmp));
1089: }
1090: PetscCall(LMBasisRestoreVecRead(orig_basis, i, &w));
1091: *scale = 1.0;
1092: *y = tmp;
1093: } else {
1094: LMBasis basis;
1095: PetscScalar dummy;
1097: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &type, &dummy));
1098: PetscCall(LMBasisGetVecRead(basis, i, y));
1099: *scale = 1.0;
1100: }
1101: break;
1102: default:
1103: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisGetVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y. Use MatLMVMGetUpdatedBasis() and LMBasisGetVecRead().");
1104: }
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: PETSC_INTERN PetscErrorCode MatLMVMBasisRestoreVecRead(Mat B, MatLMVMBasisType type, PetscInt i, Vec *y, PetscScalar *scale)
1109: {
1110: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1111: PetscBool is_scalar;
1112: PetscScalar scale_;
1114: PetscFunctionBegin;
1115: switch (type) {
1116: case LMBASIS_B0S:
1117: case LMBASIS_H0Y:
1118: // if B_0 = gamma * I, do not actually compute these bases
1119: PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &scale_));
1120: if (is_scalar) {
1121: PetscCall(LMBasisRestoreVecRead(lmvm->basis[type == LMBASIS_B0S ? LMBASIS_S : LMBASIS_Y], i, y));
1122: } else if (lmvm->do_not_cache_J0_products) {
1123: LMBasis size_basis = lmvm->basis[MatLMVMBasisSizeOf(type)];
1125: PetscCall(LMBasisRestoreWorkVec(size_basis, y));
1126: } else {
1127: PetscCall(LMBasisRestoreVecRead(lmvm->basis[type], i, y));
1128: }
1129: break;
1130: default:
1131: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "MatLMVMBasisRestoreVecRead() is only for LMBASIS_B0S and LMBASIS_H0Y. Use MatLMVMGetUpdatedBasis() and LMBasisRestoreVecRead().");
1132: }
1133: PetscFunctionReturn(PETSC_SUCCESS);
1134: }
1136: PETSC_INTERN PetscErrorCode MatLMVMGetRange(Mat B, PetscInt *oldest, PetscInt *next)
1137: {
1138: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1140: PetscFunctionBegin;
1141: PetscCall(LMBasisGetRange(lmvm->basis[LMBASIS_S], oldest, next));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: PETSC_INTERN PetscErrorCode MatLMVMGetWorkRow(Mat B, Vec *array_p)
1146: {
1147: LMBasis basis;
1149: PetscFunctionBegin;
1150: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
1151: PetscCall(LMBasisGetWorkRow(basis, array_p));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: PETSC_INTERN PetscErrorCode MatLMVMRestoreWorkRow(Mat B, Vec *array_p)
1156: {
1157: LMBasis basis;
1159: PetscFunctionBegin;
1160: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &basis, NULL, NULL));
1161: PetscCall(LMBasisRestoreWorkRow(basis, array_p));
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: 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)
1166: {
1167: LMBasis S;
1168: Vec B0x;
1170: PetscFunctionBegin;
1171: PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
1172: PetscCall(LMBasisGetWorkVec(S, &B0x));
1173: PetscCall(op(B, x, B0x));
1174: PetscCall(LMBasisGEMVH(S, oldest, next, alpha, B0x, beta, y));
1175: PetscCall(LMBasisRestoreWorkVec(S, &B0x));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: 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)
1180: {
1181: LMBasis S, Y;
1182: Vec S_x;
1183: Vec B0S_x;
1185: PetscFunctionBegin;
1186: PetscCall(MatLMVMGetUpdatedBasis(B, type_S, &S, NULL, NULL));
1187: PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, NULL, NULL));
1188: PetscCall(LMBasisGetWorkVec(S, &S_x));
1189: PetscCall(LMBasisGEMV(S, oldest, next, alpha, x, 0.0, S_x));
1190: PetscCall(LMBasisGetWorkVec(Y, &B0S_x));
1191: PetscCall(op(B, S_x, B0S_x));
1192: PetscCall(VecAYPX(y, beta, B0S_x));
1193: PetscCall(LMBasisRestoreWorkVec(Y, &B0S_x));
1194: PetscCall(LMBasisRestoreWorkVec(S, &S_x));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMVH(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec v, PetscScalar beta, Vec array)
1199: {
1200: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1201: PetscBool cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
1202: LMBasis basis;
1203: MatLMVMBasisType basis_t;
1204: PetscScalar gamma;
1206: PetscFunctionBegin;
1207: if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
1208: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &basis_t, &gamma));
1209: PetscCall(LMBasisGEMVH(basis, oldest, next, alpha * gamma, v, beta, array));
1210: } else {
1211: switch (type) {
1212: case LMBASIS_B0S:
1213: PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_S, MatLMVMApplyJ0HermitianTranspose, v, beta, array));
1214: break;
1215: case LMBASIS_H0Y:
1216: PetscCall(MatLMVMApplyOpThenVecs(alpha, B, oldest, next, LMBASIS_Y, MatLMVMApplyJ0InvHermitianTranspose, v, beta, array));
1217: break;
1218: case LMBASIS_Y_MINUS_B0S:
1219: PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_Y], oldest, next, alpha, v, beta, array));
1220: PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_B0S, oldest, next, -alpha, v, 1.0, array));
1221: break;
1222: case LMBASIS_S_MINUS_H0Y:
1223: PetscCall(LMBasisGEMVH(lmvm->basis[LMBASIS_S], oldest, next, alpha, v, beta, array));
1224: PetscCall(MatLMVMBasisGEMVH(B, LMBASIS_H0Y, oldest, next, -alpha, v, 1.0, array));
1225: break;
1226: default:
1227: PetscUnreachable();
1228: }
1229: }
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: // x must come from MatLMVMGetRowWork()
1234: PETSC_INTERN PetscErrorCode MatLMVMBasisGEMV(Mat B, MatLMVMBasisType type, PetscInt oldest, PetscInt next, PetscScalar alpha, Vec x, PetscScalar beta, Vec y)
1235: {
1236: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1237: PetscBool cache_J0_products = lmvm->do_not_cache_J0_products ? PETSC_FALSE : PETSC_TRUE;
1238: LMBasis basis;
1240: PetscFunctionBegin;
1241: if (cache_J0_products || type == LMBASIS_S || type == LMBASIS_Y) {
1242: PetscScalar gamma;
1243: MatLMVMBasisType base_type;
1245: PetscCall(MatLMVMGetUpdatedBasis(B, type, &basis, &base_type, &gamma));
1246: PetscCall(LMBasisGEMV(basis, oldest, next, alpha * gamma, x, beta, y));
1247: } else {
1248: switch (type) {
1249: case LMBASIS_B0S:
1250: PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_S, LMBASIS_Y, MatLMVMApplyJ0Fwd, x, beta, y));
1251: break;
1252: case LMBASIS_H0Y:
1253: PetscCall(MatLMVMApplyVecsThenOp(alpha, B, oldest, next, LMBASIS_Y, LMBASIS_S, MatLMVMApplyJ0Inv, x, beta, y));
1254: break;
1255: case LMBASIS_Y_MINUS_B0S:
1256: PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_Y], oldest, next, alpha, x, beta, y));
1257: PetscCall(MatLMVMBasisGEMV(B, LMBASIS_B0S, oldest, next, -alpha, x, 1.0, y));
1258: break;
1259: case LMBASIS_S_MINUS_H0Y:
1260: PetscCall(LMBasisGEMV(lmvm->basis[LMBASIS_S], oldest, next, alpha, x, beta, y));
1261: PetscCall(MatLMVMBasisGEMV(B, LMBASIS_H0Y, oldest, next, -alpha, x, 1.0, y));
1262: break;
1263: default:
1264: PetscUnreachable();
1265: }
1266: }
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: PETSC_INTERN PetscErrorCode MatLMVMCreateProducts(Mat B, LMBlockType block_type, LMProducts *products)
1271: {
1272: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1274: PetscFunctionBegin;
1275: PetscCall(LMProductsCreate(lmvm->basis[LMBASIS_S], block_type, products));
1276: (*products)->debug = lmvm->debug;
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: static PetscErrorCode MatLMVMProductsUpdate(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type)
1281: {
1282: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1283: LMBasis X, Y;
1284: MatLMVMBasisType true_type_X, true_type_Y;
1285: PetscScalar alpha_X, alpha_Y;
1286: PetscInt oldest, next;
1287: LMProducts G;
1289: PetscFunctionBegin;
1290: PetscCall(MatLMVMGetUpdatedBasis(B, type_X, &X, &true_type_X, &alpha_X));
1291: PetscCall(MatLMVMGetUpdatedBasis(B, type_Y, &Y, &true_type_Y, &alpha_Y));
1292: 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]));
1293: PetscCall(LMProductsUpdate(lmvm->products[block_type][true_type_X][true_type_Y], X, Y));
1294: if (true_type_X == type_X && true_type_Y == type_Y) PetscFunctionReturn(PETSC_SUCCESS);
1295: if (!lmvm->products[block_type][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, block_type, &lmvm->products[block_type][type_X][type_Y]));
1296: G = lmvm->products[block_type][type_X][type_Y];
1297: PetscCall(MatLMVMGetRange(B, &oldest, &next));
1298: PetscCall(LMProductsPrepare(G, lmvm->J0, oldest, next));
1299: if (G->k < lmvm->k) {
1300: PetscCall(LMProductsCopy(lmvm->products[block_type][true_type_X][true_type_Y], lmvm->products[block_type][type_X][type_Y]));
1301: if (alpha_X * alpha_Y != 1.0) PetscCall(LMProductsScale(lmvm->products[block_type][type_X][type_Y], alpha_X * alpha_Y));
1302: }
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: PETSC_INTERN PetscErrorCode MatLMVMGetUpdatedProducts(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, LMBlockType block_type, LMProducts *lmwd)
1307: {
1308: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1310: PetscFunctionBegin;
1311: PetscCall(MatLMVMProductsUpdate(B, type_X, type_Y, block_type));
1312: *lmwd = lmvm->products[block_type][type_X][type_Y];
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: PETSC_INTERN PetscErrorCode MatLMVMProductsInsertDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar v)
1317: {
1318: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1319: LMProducts products;
1321: PetscFunctionBegin;
1322: if (!lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y]));
1323: products = lmvm->products[LMBLOCK_DIAGONAL][type_X][type_Y];
1324: PetscCall(LMProductsInsertNextDiagonalValue(products, idx, v));
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: PETSC_INTERN PetscErrorCode MatLMVMProductsGetDiagonalValue(Mat B, MatLMVMBasisType type_X, MatLMVMBasisType type_Y, PetscInt idx, PetscScalar *v)
1329: {
1330: LMProducts products = NULL;
1332: PetscFunctionBegin;
1333: PetscCall(MatLMVMGetUpdatedProducts(B, type_X, type_Y, LMBLOCK_DIAGONAL, &products));
1334: PetscCall(LMProductsGetDiagonalValue(products, idx, v));
1335: PetscFunctionReturn(PETSC_SUCCESS);
1336: }