Actual source code: lmvmutils.c
1: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: /*@
4: MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to an LMVM-type matrix.
5: The first time the function is called for an LMVM-type matrix, no update is
6: applied, but the given X and F vectors are stored for use as Xprev and
7: Fprev in the next update.
9: If the user has provided another LMVM-type matrix in place of J0, the J0
10: matrix is also updated recursively.
12: Input Parameters:
13: + B - An LMVM-type matrix
14: . X - Solution vector
15: - F - Function vector
17: Level: intermediate
19: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMAllocate()`
20: @*/
21: PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
22: {
23: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
24: PetscBool same;
26: PetscFunctionBegin;
30: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
31: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
32: if (!lmvm->allocated) {
33: PetscCall(MatLMVMAllocate(B, X, F));
34: } else {
35: VecCheckMatCompatible(B, X, 2, F, 3);
36: }
37: if (lmvm->J0) {
38: /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
39: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
40: if (same) PetscCall(MatLMVMUpdate(lmvm->J0, X, F));
41: }
42: PetscCall((*lmvm->ops->update)(B, X, F));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
48: an identity matrix (scale = 1.0).
50: Input Parameter:
51: . B - An LMVM-type matrix
53: Level: advanced
55: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
56: @*/
57: PetscErrorCode MatLMVMClearJ0(Mat B)
58: {
59: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
60: PetscBool same;
62: PetscFunctionBegin;
64: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
65: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
66: lmvm->user_pc = PETSC_FALSE;
67: lmvm->user_ksp = PETSC_FALSE;
68: lmvm->user_scale = PETSC_FALSE;
69: lmvm->J0scalar = 1.0;
70: PetscCall(VecDestroy(&lmvm->J0diag));
71: PetscCall(MatDestroy(&lmvm->J0));
72: PetscCall(PCDestroy(&lmvm->J0pc));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: MatLMVMSetJ0Scale - Allows the user to define a scalar value
78: mu such that J0 = mu*I.
80: Input Parameters:
81: + B - An LMVM-type matrix
82: - scale - Scalar value mu that defines the initial Jacobian
84: Level: advanced
86: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
87: @*/
88: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
89: {
90: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
91: PetscBool same;
93: PetscFunctionBegin;
95: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
96: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
97: PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Scaling is available only for square LMVM matrices");
98: PetscCall(MatLMVMClearJ0(B));
99: lmvm->J0scalar = scale;
100: lmvm->user_scale = PETSC_TRUE;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: MatLMVMSetJ0Diag - Allows the user to define a vector
106: V such that J0 = diag(V).
108: Input Parameters:
109: + B - An LMVM-type matrix
110: - V - Vector that defines the diagonal of the initial Jacobian
112: Level: advanced
114: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
115: @*/
116: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
117: {
118: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
119: PetscBool same;
120: MPI_Comm comm = PetscObjectComm((PetscObject)B);
122: PetscFunctionBegin;
125: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
126: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
127: PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "Matrix must be allocated before setting diagonal scaling");
128: PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Diagonal scaling is available only for square LMVM matrices");
129: VecCheckSameSize(V, 2, lmvm->Fprev, 3);
131: PetscCall(MatLMVMClearJ0(B));
132: if (!lmvm->J0diag) PetscCall(VecDuplicate(V, &lmvm->J0diag));
133: PetscCall(VecCopy(V, lmvm->J0diag));
134: lmvm->user_scale = PETSC_TRUE;
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: /*@
139: MatLMVMSetJ0 - Allows the user to define the initial
140: Jacobian matrix from which the LMVM-type approximation is
141: built up. Inverse of this initial Jacobian is applied
142: using an internal `KSP` solver, which defaults to `KSPGMRES`.
143: This internal `KSP` solver has the "mat_lmvm_" option
144: prefix.
146: Note that another LMVM-type matrix can be used in place of
147: J0, in which case updating the outer LMVM-type matrix will
148: also trigger the update for the inner LMVM-type matrix. This
149: is useful in cases where a full-memory diagonal approximation
150: such as `MATLMVMDIAGBRDN` is used in place of J0.
152: Input Parameters:
153: + B - An LMVM-type matrix
154: - J0 - The initial Jacobian matrix
156: Level: advanced
158: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
159: @*/
160: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
161: {
162: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
163: PetscBool same;
165: PetscFunctionBegin;
168: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
169: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
170: PetscCall(MatLMVMClearJ0(B));
171: PetscCall(MatDestroy(&lmvm->J0));
172: PetscCall(PetscObjectReference((PetscObject)J0));
173: lmvm->J0 = J0;
174: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
175: if (!same && lmvm->square) PetscCall(KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /*@
180: MatLMVMSetJ0PC - Allows the user to define a `PC` object that
181: acts as the initial inverse-Jacobian matrix. This `PC` should
182: already contain all the operators necessary for its application.
183: The LMVM-type matrix only calls `PCApply()` without changing any other
184: options.
186: Input Parameters:
187: + B - An LMVM-type matrix
188: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
190: Level: advanced
192: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
193: @*/
194: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
195: {
196: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
197: PetscBool same;
198: MPI_Comm comm = PetscObjectComm((PetscObject)B);
200: PetscFunctionBegin;
203: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
204: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
205: PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
206: PetscCall(MatLMVMClearJ0(B));
207: PetscCall(PetscObjectReference((PetscObject)J0pc));
208: lmvm->J0pc = J0pc;
209: lmvm->user_pc = PETSC_TRUE;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@
214: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
215: KSP solver for the initial inverse-Jacobian approximation.
216: This `KSP` solver should already contain all the operators
217: necessary to perform the inversion. The LMVM-type matrix only
218: calls `KSPSolve()` without changing any other options.
220: Input Parameters:
221: + B - An LMVM-type matrix
222: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
224: Level: advanced
226: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
227: @*/
228: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
229: {
230: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
231: PetscBool same;
232: MPI_Comm comm = PetscObjectComm((PetscObject)B);
234: PetscFunctionBegin;
237: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
238: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
239: PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
240: PetscCall(MatLMVMClearJ0(B));
241: PetscCall(KSPDestroy(&lmvm->J0ksp));
242: PetscCall(PetscObjectReference((PetscObject)J0ksp));
243: lmvm->J0ksp = J0ksp;
244: lmvm->user_ksp = PETSC_TRUE;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.
251: Input Parameters:
252: . B - An LMVM-type matrix
254: Output Parameter:
255: . J0 - `Mat` object for defining the initial Jacobian
257: Level: advanced
259: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
260: @*/
261: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
262: {
263: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
264: PetscBool same;
266: PetscFunctionBegin;
268: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
269: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
270: *J0 = lmvm->J0;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
276: associated with the initial Jacobian.
278: Input Parameter:
279: . B - An LMVM-type matrix
281: Output Parameter:
282: . J0pc - `PC` object for defining the initial inverse-Jacobian
284: Level: advanced
286: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
287: @*/
288: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
289: {
290: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
291: PetscBool same;
293: PetscFunctionBegin;
295: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
296: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
297: if (lmvm->J0pc) {
298: *J0pc = lmvm->J0pc;
299: } else {
300: PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
301: }
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
307: associated with the initial Jacobian.
309: Input Parameter:
310: . B - An LMVM-type matrix
312: Output Parameter:
313: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
315: Level: advanced
317: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
318: @*/
319: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
320: {
321: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
322: PetscBool same;
324: PetscFunctionBegin;
326: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
327: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
328: *J0ksp = lmvm->J0ksp;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*@
333: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
334: matrix-vector product with the initial Jacobian.
336: Input Parameters:
337: + B - An LMVM-type matrix
338: - X - vector to multiply with J0
340: Output Parameter:
341: . Y - resulting vector for the operation
343: Level: advanced
345: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
346: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
347: @*/
348: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
349: {
350: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
351: PetscBool same, hasMult;
352: MPI_Comm comm = PetscObjectComm((PetscObject)B);
353: Mat Amat, Pmat;
355: PetscFunctionBegin;
359: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
360: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
361: PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
362: VecCheckMatCompatible(B, X, 2, Y, 3);
363: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
364: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
365: if (lmvm->user_pc) {
366: PetscCall(PCGetOperators(lmvm->J0pc, &Amat, &Pmat));
367: } else if (lmvm->user_ksp) {
368: PetscCall(KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat));
369: } else {
370: Amat = lmvm->J0;
371: }
372: PetscCall(MatHasOperation(Amat, MATOP_MULT, &hasMult));
373: if (hasMult) {
374: /* product is available, use it */
375: PetscCall(MatMult(Amat, X, Y));
376: } else {
377: /* there's no product, so treat J0 as identity */
378: PetscCall(VecCopy(X, Y));
379: }
380: } else if (lmvm->user_scale) {
381: if (lmvm->J0diag) {
382: /* User has defined a diagonal vector for J0 */
383: PetscCall(VecPointwiseMult(X, lmvm->J0diag, Y));
384: } else {
385: /* User has defined a scalar value for J0 */
386: PetscCall(VecAXPBY(Y, lmvm->J0scalar, 0.0, X));
387: }
388: } else {
389: /* There is no J0 representation so just apply an identity matrix */
390: PetscCall(VecCopy(X, Y));
391: }
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: /*@
396: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
397: inverse to the given vector. The specific form of the application
398: depends on whether the user provided a scaling factor, a J0 matrix,
399: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
400: provided, the function simply does an identity matrix application
401: (vector copy).
403: Input Parameters:
404: + B - An LMVM-type matrix
405: - X - vector to "multiply" with J0^{-1}
407: Output Parameter:
408: . Y - resulting vector for the operation
410: Level: advanced
412: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
413: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
414: @*/
415: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
416: {
417: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
418: PetscBool same, hasSolve;
419: MPI_Comm comm = PetscObjectComm((PetscObject)B);
421: PetscFunctionBegin;
425: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
426: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
427: PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
428: VecCheckMatCompatible(B, X, 2, Y, 3);
430: /* Invert the initial Jacobian onto q (or apply scaling) */
431: if (lmvm->user_pc) {
432: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
433: PetscCall(PCApply(lmvm->J0pc, X, Y));
434: } else if (lmvm->user_ksp) {
435: /* User has defined a J0 or a custom KSP so just perform a solution */
436: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
437: } else if (lmvm->J0) {
438: PetscCall(MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve));
439: if (hasSolve) {
440: PetscCall(MatSolve(lmvm->J0, X, Y));
441: } else {
442: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
443: }
444: } else if (lmvm->user_scale) {
445: if (lmvm->J0diag) {
446: PetscCall(VecPointwiseDivide(X, Y, lmvm->J0diag));
447: } else {
448: PetscCall(VecAXPBY(Y, 1.0 / lmvm->J0scalar, 0.0, X));
449: }
450: } else {
451: /* There is no J0 representation so just apply an identity matrix */
452: PetscCall(VecCopy(X, Y));
453: }
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: MatLMVMIsAllocated - Returns a boolean flag that shows whether
459: the necessary data structures for the underlying matrix is allocated.
461: Input Parameter:
462: . B - An LMVM-type matrix
464: Output Parameter:
465: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
467: Level: intermediate
469: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
470: @*/
472: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
473: {
474: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
475: PetscBool same;
477: PetscFunctionBegin;
479: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
480: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
481: *flg = PETSC_FALSE;
482: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /*@
487: MatLMVMAllocate - Produces all necessary common memory for
488: LMVM approximations based on the solution and function vectors
489: provided. If `MatSetSizes()` and `MatSetUp()` have not been called
490: before `MatLMVMAllocate()`, the allocation will read sizes from
491: the provided vectors and update the matrix.
493: Input Parameters:
494: + B - An LMVM-type matrix
495: . X - Solution vector
496: - F - Function vector
498: Level: intermediate
500: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
501: @*/
502: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
503: {
504: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
505: PetscBool same;
506: VecType vtype;
508: PetscFunctionBegin;
512: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
513: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
514: PetscCall(VecGetType(X, &vtype));
515: PetscCall(MatSetVecType(B, vtype));
516: PetscCall((*lmvm->ops->allocate)(B, X, F));
517: if (lmvm->J0) {
518: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
519: if (same) PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
520: }
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: MatLMVMResetShift - Zero the shift factor.
527: Input Parameter:
528: . B - An LMVM-type matrix
530: Level: intermediate
532: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
533: @*/
534: PetscErrorCode MatLMVMResetShift(Mat B)
535: {
536: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
537: PetscBool same;
539: PetscFunctionBegin;
541: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
542: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
543: lmvm->shift = 0.0;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: MatLMVMReset - Flushes all of the accumulated updates out of
549: the LMVM approximation. In practice, this will not actually
550: destroy the data associated with the updates. It simply resets
551: counters, which leads to existing data being overwritten, and
552: `MatSolve()` being applied as if there are no updates. A boolean
553: flag is available to force destruction of the update vectors.
555: If the user has provided another LMVM matrix as J0, the J0
556: matrix is also reset in this function.
558: Input Parameters:
559: + B - An LMVM-type matrix
560: - destructive - flag for enabling destruction of data structures
562: Level: intermediate
564: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
565: @*/
566: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
567: {
568: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
569: PetscBool same;
571: PetscFunctionBegin;
573: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
574: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
575: PetscCall((*lmvm->ops->reset)(B, destructive));
576: if (lmvm->J0) {
577: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
578: if (same) PetscCall(MatLMVMReset(lmvm->J0, destructive));
579: }
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: /*@
584: MatLMVMSetHistorySize - Set the number of past iterates to be
585: stored for the construction of the limited-memory QN update.
587: Input Parameters:
588: + B - An LMVM-type matrix
589: - hist_size - number of past iterates (default 5)
591: Options Database Key:
592: . -mat_lmvm_hist_size <m> - set number of past iterates
594: Level: beginner
596: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
597: @*/
598: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
599: {
600: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
601: PetscBool same;
602: Vec X, F;
604: PetscFunctionBegin;
606: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
607: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
608: if (hist_size > 0) {
609: lmvm->m = hist_size;
610: if (lmvm->allocated && lmvm->m != lmvm->m_old) {
611: PetscCall(VecDuplicate(lmvm->Xprev, &X));
612: PetscCall(VecDuplicate(lmvm->Fprev, &F));
613: PetscCall(MatLMVMReset(B, PETSC_TRUE));
614: PetscCall(MatLMVMAllocate(B, X, F));
615: PetscCall(VecDestroy(&X));
616: PetscCall(VecDestroy(&F));
617: }
618: } else PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /*@
623: MatLMVMGetUpdateCount - Returns the number of accepted updates.
624: This number may be greater than the total number of update vectors
625: stored in the matrix. The counters are reset when `MatLMVMReset()`
626: is called.
628: Input Parameter:
629: . B - An LMVM-type matrix
631: Output Parameter:
632: . nupdates - number of accepted updates
634: Level: intermediate
636: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
637: @*/
638: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
639: {
640: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
641: PetscBool same;
643: PetscFunctionBegin;
645: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
646: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
647: *nupdates = lmvm->nupdates;
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
651: /*@
652: MatLMVMGetRejectCount - Returns the number of rejected updates.
653: The counters are reset when `MatLMVMReset()` is called.
655: Input Parameter:
656: . B - An LMVM-type matrix
658: Output Parameter:
659: . nrejects - number of rejected updates
661: Level: intermediate
663: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
664: @*/
665: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
666: {
667: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
668: PetscBool same;
670: PetscFunctionBegin;
672: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
673: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
674: *nrejects = lmvm->nrejects;
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }