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 `MATLMVM` matrix.
5: The first time the function is called for an `MATLMVM` 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 `MATLMVM` matrix in place of J0, the J0
10: matrix is also updated recursively.
12: Input Parameters:
13: + B - A `MATLMVM` 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 - A `MATLMVM` 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 - A `MATLMVM` 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 - A `MATLMVM` 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.
143: Input Parameters:
144: + B - A `MATLMVM` matrix
145: - J0 - The initial Jacobian matrix
147: Level: advanced
149: Notes:
150: The inverse of this initial Jacobian is applied
151: using an internal `KSP` solver, which defaults to `KSPGMRES`.
153: This internal `KSP` solver has the "mat_lmvm_" option
154: prefix.
156: Another `MATLMVM` matrix can be used in place of
157: `J0`, in which case updating the outer `MATLMVM` matrix will
158: also trigger the update for the inner `MATLMVM` matrix. This
159: is useful in cases where a full-memory diagonal approximation
160: such as `MATLMVMDIAGBRDN` is used in place of `J0`.
162: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
163: @*/
164: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
165: {
166: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
167: PetscBool same;
169: PetscFunctionBegin;
172: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
173: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
174: PetscCall(MatLMVMClearJ0(B));
175: PetscCall(MatDestroy(&lmvm->J0));
176: PetscCall(PetscObjectReference((PetscObject)J0));
177: lmvm->J0 = J0;
178: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
179: if (!same && lmvm->square) PetscCall(KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: MatLMVMSetJ0PC - Allows the user to define a `PC` object that
185: acts as the initial inverse-Jacobian matrix.
187: Input Parameters:
188: + B - A `MATLMVM` matrix
189: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
191: Level: advanced
193: Note:
194: `J0pc` should already contain all the operators necessary for its application.
195: The `MATLMVM` matrix only calls `PCApply()` without changing any other
196: options.
198: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
199: @*/
200: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
201: {
202: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
203: PetscBool same;
204: MPI_Comm comm = PetscObjectComm((PetscObject)B);
206: PetscFunctionBegin;
209: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
210: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
211: PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
212: PetscCall(MatLMVMClearJ0(B));
213: PetscCall(PetscObjectReference((PetscObject)J0pc));
214: lmvm->J0pc = J0pc;
215: lmvm->user_pc = PETSC_TRUE;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
221: KSP solver for the initial inverse-Jacobian approximation.
223: Input Parameters:
224: + B - A `MATLMVM` matrix
225: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
227: Level: advanced
229: Note:
230: The `KSP` solver should already contain all the operators
231: necessary to perform the inversion. The `MATLMVM` matrix only
232: calls `KSPSolve()` without changing any other options.
234: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
235: @*/
236: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
237: {
238: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
239: PetscBool same;
240: MPI_Comm comm = PetscObjectComm((PetscObject)B);
242: PetscFunctionBegin;
245: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
246: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
247: PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
248: PetscCall(MatLMVMClearJ0(B));
249: PetscCall(KSPDestroy(&lmvm->J0ksp));
250: PetscCall(PetscObjectReference((PetscObject)J0ksp));
251: lmvm->J0ksp = J0ksp;
252: lmvm->user_ksp = PETSC_TRUE;
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: MatLMVMGetJ0 - Returns a pointer to the internal `J0` matrix.
259: Input Parameter:
260: . B - A `MATLMVM` matrix
262: Output Parameter:
263: . J0 - `Mat` object for defining the initial Jacobian
265: Level: advanced
267: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
268: @*/
269: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
270: {
271: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
272: PetscBool same;
274: PetscFunctionBegin;
276: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
277: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
278: *J0 = lmvm->J0;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@
283: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
284: associated with the initial Jacobian.
286: Input Parameter:
287: . B - A `MATLMVM` matrix
289: Output Parameter:
290: . J0pc - `PC` object for defining the initial inverse-Jacobian
292: Level: advanced
294: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
295: @*/
296: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
297: {
298: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
299: PetscBool same;
301: PetscFunctionBegin;
303: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
304: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
305: if (lmvm->J0pc) {
306: *J0pc = lmvm->J0pc;
307: } else {
308: PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
309: }
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
315: associated with the initial Jacobian.
317: Input Parameter:
318: . B - A `MATLMVM` matrix
320: Output Parameter:
321: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
323: Level: advanced
325: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
326: @*/
327: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
328: {
329: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
330: PetscBool same;
332: PetscFunctionBegin;
334: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
335: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
336: *J0ksp = lmvm->J0ksp;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
342: matrix-vector product with the initial Jacobian.
344: Input Parameters:
345: + B - A `MATLMVM` matrix
346: - X - vector to multiply with J0
348: Output Parameter:
349: . Y - resulting vector for the operation
351: Level: advanced
353: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
354: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
355: @*/
356: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
357: {
358: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
359: PetscBool same, hasMult;
360: MPI_Comm comm = PetscObjectComm((PetscObject)B);
361: Mat Amat, Pmat;
363: PetscFunctionBegin;
367: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
368: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
369: PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
370: VecCheckMatCompatible(B, X, 2, Y, 3);
371: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
372: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
373: if (lmvm->user_pc) {
374: PetscCall(PCGetOperators(lmvm->J0pc, &Amat, &Pmat));
375: } else if (lmvm->user_ksp) {
376: PetscCall(KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat));
377: } else {
378: Amat = lmvm->J0;
379: }
380: PetscCall(MatHasOperation(Amat, MATOP_MULT, &hasMult));
381: if (hasMult) {
382: /* product is available, use it */
383: PetscCall(MatMult(Amat, X, Y));
384: } else {
385: /* there's no product, so treat J0 as identity */
386: PetscCall(VecCopy(X, Y));
387: }
388: } else if (lmvm->user_scale) {
389: if (lmvm->J0diag) {
390: /* User has defined a diagonal vector for J0 */
391: PetscCall(VecPointwiseMult(X, lmvm->J0diag, Y));
392: } else {
393: /* User has defined a scalar value for J0 */
394: PetscCall(VecAXPBY(Y, lmvm->J0scalar, 0.0, X));
395: }
396: } else {
397: /* There is no J0 representation so just apply an identity matrix */
398: PetscCall(VecCopy(X, Y));
399: }
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@
404: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
405: inverse to the given vector.
407: Input Parameters:
408: + B - A `MATLMVM` matrix
409: - X - vector to "multiply" with J0^{-1}
411: Output Parameter:
412: . Y - resulting vector for the operation
414: Level: advanced
416: Note:
417: The specific form of the application
418: depends on whether the user provided a scaling factor, a J0 matrix,
419: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
420: provided, the function simply does an identity matrix application
421: (vector copy).
423: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
424: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
425: @*/
426: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
427: {
428: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
429: PetscBool same, hasSolve;
430: MPI_Comm comm = PetscObjectComm((PetscObject)B);
432: PetscFunctionBegin;
436: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
437: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
438: PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
439: VecCheckMatCompatible(B, X, 2, Y, 3);
441: /* Invert the initial Jacobian onto q (or apply scaling) */
442: if (lmvm->user_pc) {
443: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
444: PetscCall(PCApply(lmvm->J0pc, X, Y));
445: } else if (lmvm->user_ksp) {
446: /* User has defined a J0 or a custom KSP so just perform a solution */
447: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
448: } else if (lmvm->J0) {
449: PetscCall(MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve));
450: if (hasSolve) {
451: PetscCall(MatSolve(lmvm->J0, X, Y));
452: } else {
453: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
454: }
455: } else if (lmvm->user_scale) {
456: if (lmvm->J0diag) {
457: PetscCall(VecPointwiseDivide(X, Y, lmvm->J0diag));
458: } else {
459: PetscCall(VecAXPBY(Y, 1.0 / lmvm->J0scalar, 0.0, X));
460: }
461: } else {
462: /* There is no J0 representation so just apply an identity matrix */
463: PetscCall(VecCopy(X, Y));
464: }
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: MatLMVMIsAllocated - Returns a boolean flag that shows whether
470: the necessary data structures for the underlying matrix is allocated.
472: Input Parameter:
473: . B - A `MATLMVM` matrix
475: Output Parameter:
476: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
478: Level: intermediate
480: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
481: @*/
482: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
483: {
484: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
485: PetscBool same;
487: PetscFunctionBegin;
489: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
490: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
491: *flg = PETSC_FALSE;
492: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: /*@
497: MatLMVMAllocate - Produces all necessary common memory for
498: LMVM approximations based on the solution and function vectors
499: provided.
501: Input Parameters:
502: + B - A `MATLMVM` matrix
503: . X - Solution vector
504: - F - Function vector
506: Level: intermediate
508: Note:
509: If `MatSetSizes()` and `MatSetUp()` have not been called
510: before `MatLMVMAllocate()`, the allocation will read sizes from
511: the provided vectors and update the matrix.
513: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
514: @*/
515: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
516: {
517: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
518: PetscBool same;
519: VecType vtype;
521: PetscFunctionBegin;
525: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
526: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
527: PetscCall(VecGetType(X, &vtype));
528: PetscCall(MatSetVecType(B, vtype));
529: PetscCall((*lmvm->ops->allocate)(B, X, F));
530: if (lmvm->J0) {
531: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
532: if (same) PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
533: }
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.
540: Input Parameter:
541: . B - A `MATLMVM` matrix
543: Level: intermediate
545: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
546: @*/
547: PetscErrorCode MatLMVMResetShift(Mat B)
548: {
549: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
550: PetscBool same;
552: PetscFunctionBegin;
554: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
555: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
556: lmvm->shift = 0.0;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /*@
561: MatLMVMReset - Flushes all of the accumulated updates out of
562: the `MATLMVM` approximation.
564: Input Parameters:
565: + B - A `MATLMVM` matrix
566: - destructive - flag for enabling destruction of data structures
568: Level: intermediate
570: Note:
571: In practice, this will not actually
572: destroy the data associated with the updates. It simply resets
573: counters, which leads to existing data being overwritten, and
574: `MatSolve()` being applied as if there are no updates. A boolean
575: flag is available to force destruction of the update vectors.
577: If the user has provided another LMVM matrix as J0, the J0
578: matrix is also reset in this function.
580: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
581: @*/
582: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
583: {
584: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
585: PetscBool same;
587: PetscFunctionBegin;
589: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
590: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
591: PetscCall((*lmvm->ops->reset)(B, destructive));
592: if (lmvm->J0) {
593: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
594: if (same) PetscCall(MatLMVMReset(lmvm->J0, destructive));
595: }
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*@
600: MatLMVMSetHistorySize - Set the number of past iterates to be
601: stored for the construction of the limited-memory quasi-Newton update.
603: Input Parameters:
604: + B - A `MATLMVM` matrix
605: - hist_size - number of past iterates (default 5)
607: Options Database Key:
608: . -mat_lmvm_hist_size <m> - set number of past iterates
610: Level: beginner
612: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
613: @*/
614: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
615: {
616: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
617: PetscBool same;
618: Vec X, F;
620: PetscFunctionBegin;
622: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
623: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
624: if (hist_size > 0) {
625: lmvm->m = hist_size;
626: if (lmvm->allocated && lmvm->m != lmvm->m_old) {
627: PetscCall(VecDuplicate(lmvm->Xprev, &X));
628: PetscCall(VecDuplicate(lmvm->Fprev, &F));
629: PetscCall(MatLMVMReset(B, PETSC_TRUE));
630: PetscCall(MatLMVMAllocate(B, X, F));
631: PetscCall(VecDestroy(&X));
632: PetscCall(VecDestroy(&F));
633: }
634: } else PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@
639: MatLMVMGetUpdateCount - Returns the number of accepted updates.
641: Input Parameter:
642: . B - A `MATLMVM` matrix
644: Output Parameter:
645: . nupdates - number of accepted updates
647: Level: intermediate
649: Note:
650: This number may be greater than the total number of update vectors
651: stored in the matrix. The counters are reset when `MatLMVMReset()`
652: is called.
654: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
655: @*/
656: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
657: {
658: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
659: PetscBool same;
661: PetscFunctionBegin;
663: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
664: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
665: *nupdates = lmvm->nupdates;
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: MatLMVMGetRejectCount - Returns the number of rejected updates.
671: The counters are reset when `MatLMVMReset()` is called.
673: Input Parameter:
674: . B - A `MATLMVM` matrix
676: Output Parameter:
677: . nrejects - number of rejected updates
679: Level: intermediate
681: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
682: @*/
683: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
684: {
685: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
686: PetscBool same;
688: PetscFunctionBegin;
690: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
691: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
692: *nrejects = lmvm->nrejects;
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }