Actual source code: lmvmutils.c
1: #include <petscdevice.h>
2: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: #include <petsc/private/deviceimpl.h>
5: /*@
6: MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to an `MATLMVM` matrix.
7: The first time the function is called for an `MATLMVM` matrix, no update is
8: applied, but the given X and F vectors are stored for use as Xprev and
9: Fprev in the next update.
11: If the user has provided another `MATLMVM` matrix in place of J0, the J0
12: matrix is also updated recursively.
14: Input Parameters:
15: + B - A `MATLMVM` matrix
16: . X - Solution vector
17: - F - Function vector
19: Level: intermediate
21: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMAllocate()`
22: @*/
23: PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
24: {
25: Mat_LMVM *lmvm;
26: PetscBool same;
28: PetscFunctionBegin;
32: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
33: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
34: lmvm = (Mat_LMVM *)B->data;
35: if (!lmvm->allocated) {
36: PetscCall(MatLMVMAllocate(B, X, F));
37: } else {
38: VecCheckMatCompatible(B, X, 2, F, 3);
39: }
40: if (lmvm->J0) {
41: /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
42: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
43: if (same) PetscCall(MatLMVMUpdate(lmvm->J0, X, F));
44: }
45: PetscCall((*lmvm->ops->update)(B, X, F));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: /*@
50: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
51: an identity matrix (scale = 1.0).
53: Input Parameter:
54: . B - A `MATLMVM` matrix
56: Level: advanced
58: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
59: @*/
60: PetscErrorCode MatLMVMClearJ0(Mat B)
61: {
62: Mat_LMVM *lmvm;
63: PetscBool same;
65: PetscFunctionBegin;
67: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
68: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
69: lmvm = (Mat_LMVM *)B->data;
70: lmvm->user_pc = PETSC_FALSE;
71: lmvm->user_ksp = PETSC_FALSE;
72: lmvm->user_scale = PETSC_FALSE;
73: lmvm->J0scalar = 1.0;
74: PetscCall(VecDestroy(&lmvm->J0diag));
75: PetscCall(MatDestroy(&lmvm->J0));
76: PetscCall(PCDestroy(&lmvm->J0pc));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /*@
81: MatLMVMSetJ0Scale - Allows the user to define a scalar value
82: mu such that J0 = mu*I.
84: Input Parameters:
85: + B - A `MATLMVM` matrix
86: - scale - Scalar value mu that defines the initial Jacobian
88: Level: advanced
90: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
91: @*/
92: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
93: {
94: Mat_LMVM *lmvm;
95: PetscBool same;
97: PetscFunctionBegin;
99: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
100: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
101: lmvm = (Mat_LMVM *)B->data;
102: PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Scaling is available only for square LMVM matrices");
103: PetscCall(MatLMVMClearJ0(B));
104: lmvm->J0scalar = scale;
105: lmvm->user_scale = PETSC_TRUE;
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*@
110: MatLMVMSetJ0Diag - Allows the user to define a vector
111: V such that J0 = diag(V).
113: Input Parameters:
114: + B - A `MATLMVM` matrix
115: - V - Vector that defines the diagonal of the initial Jacobian
117: Level: advanced
119: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
120: @*/
121: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
122: {
123: Mat_LMVM *lmvm;
124: PetscBool same;
126: PetscFunctionBegin;
129: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
130: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
131: lmvm = (Mat_LMVM *)B->data;
132: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "Matrix must be allocated before setting diagonal scaling");
133: PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Diagonal scaling is available only for square LMVM matrices");
134: VecCheckSameSize(V, 2, lmvm->Fprev, 3);
136: PetscCall(MatLMVMClearJ0(B));
137: if (!lmvm->J0diag) PetscCall(VecDuplicate(V, &lmvm->J0diag));
138: PetscCall(VecCopy(V, lmvm->J0diag));
139: lmvm->user_scale = PETSC_TRUE;
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@
144: MatLMVMSetJ0 - Allows the user to define the initial
145: Jacobian matrix from which the LMVM-type approximation is
146: built up.
148: Input Parameters:
149: + B - A `MATLMVM` matrix
150: - J0 - The initial Jacobian matrix
152: Level: advanced
154: Notes:
155: The inverse of this initial Jacobian is applied
156: using an internal `KSP` solver, which defaults to `KSPGMRES`.
158: This internal `KSP` solver has the "mat_lmvm_" option
159: prefix.
161: Another `MATLMVM` matrix can be used in place of
162: `J0`, in which case updating the outer `MATLMVM` matrix will
163: also trigger the update for the inner `MATLMVM` matrix. This
164: is useful in cases where a full-memory diagonal approximation
165: such as `MATLMVMDIAGBRDN` is used in place of `J0`.
167: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
168: @*/
169: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
170: {
171: Mat_LMVM *lmvm;
172: PetscBool same;
174: PetscFunctionBegin;
177: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
178: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
179: lmvm = (Mat_LMVM *)B->data;
180: PetscCall(MatLMVMClearJ0(B));
181: PetscCall(MatDestroy(&lmvm->J0));
182: PetscCall(PetscObjectReference((PetscObject)J0));
183: lmvm->J0 = J0;
184: PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
185: if (!same && lmvm->square) PetscCall(KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*@
190: MatLMVMSetJ0PC - Allows the user to define a `PC` object that
191: acts as the initial inverse-Jacobian matrix.
193: Input Parameters:
194: + B - A `MATLMVM` matrix
195: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
197: Level: advanced
199: Note:
200: `J0pc` should already contain all the operators necessary for its application.
201: The `MATLMVM` matrix only calls `PCApply()` without changing any other
202: options.
204: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
205: @*/
206: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
207: {
208: Mat_LMVM *lmvm;
209: PetscBool same;
211: PetscFunctionBegin;
214: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
215: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
216: lmvm = (Mat_LMVM *)B->data;
217: PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
218: PetscCall(MatLMVMClearJ0(B));
219: PetscCall(PetscObjectReference((PetscObject)J0pc));
220: lmvm->J0pc = J0pc;
221: lmvm->user_pc = PETSC_TRUE;
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*@
226: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
227: KSP solver for the initial inverse-Jacobian approximation.
229: Input Parameters:
230: + B - A `MATLMVM` matrix
231: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
233: Level: advanced
235: Note:
236: The `KSP` solver should already contain all the operators
237: necessary to perform the inversion. The `MATLMVM` matrix only
238: calls `KSPSolve()` without changing any other options.
240: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
241: @*/
242: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
243: {
244: Mat_LMVM *lmvm;
245: PetscBool same;
247: PetscFunctionBegin;
250: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
251: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
252: lmvm = (Mat_LMVM *)B->data;
253: PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
254: PetscCall(MatLMVMClearJ0(B));
255: PetscCall(KSPDestroy(&lmvm->J0ksp));
256: PetscCall(PetscObjectReference((PetscObject)J0ksp));
257: lmvm->J0ksp = J0ksp;
258: lmvm->user_ksp = PETSC_TRUE;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*@
263: MatLMVMGetJ0 - Returns a pointer to the internal `J0` matrix.
265: Input Parameter:
266: . B - A `MATLMVM` matrix
268: Output Parameter:
269: . J0 - `Mat` object for defining the initial Jacobian
271: Level: advanced
273: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
274: @*/
275: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
276: {
277: Mat_LMVM *lmvm;
278: PetscBool same;
280: PetscFunctionBegin;
282: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
283: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
284: lmvm = (Mat_LMVM *)B->data;
285: *J0 = lmvm->J0;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
291: associated with the initial Jacobian.
293: Input Parameter:
294: . B - A `MATLMVM` matrix
296: Output Parameter:
297: . J0pc - `PC` object for defining the initial inverse-Jacobian
299: Level: advanced
301: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
302: @*/
303: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
304: {
305: Mat_LMVM *lmvm;
306: PetscBool same;
308: PetscFunctionBegin;
310: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
311: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
312: lmvm = (Mat_LMVM *)B->data;
313: if (lmvm->J0pc) {
314: *J0pc = lmvm->J0pc;
315: } else {
316: PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
317: }
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@
322: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
323: associated with the initial Jacobian.
325: Input Parameter:
326: . B - A `MATLMVM` matrix
328: Output Parameter:
329: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
331: Level: advanced
333: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
334: @*/
335: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
336: {
337: Mat_LMVM *lmvm;
338: PetscBool same;
340: PetscFunctionBegin;
342: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
343: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
344: lmvm = (Mat_LMVM *)B->data;
345: *J0ksp = lmvm->J0ksp;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@
350: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
351: matrix-vector product with the initial Jacobian.
353: Input Parameters:
354: + B - A `MATLMVM` matrix
355: - X - vector to multiply with J0
357: Output Parameter:
358: . Y - resulting vector for the operation
360: Level: advanced
362: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
363: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
364: @*/
365: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
366: {
367: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
368: PetscBool same, hasMult;
369: Mat Amat, Pmat;
371: PetscFunctionBegin;
375: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
376: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
377: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
378: VecCheckMatCompatible(B, X, 2, Y, 3);
379: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
380: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
381: if (lmvm->user_pc) {
382: PetscCall(PCGetOperators(lmvm->J0pc, &Amat, &Pmat));
383: } else if (lmvm->user_ksp) {
384: PetscCall(KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat));
385: } else {
386: Amat = lmvm->J0;
387: }
388: PetscCall(MatHasOperation(Amat, MATOP_MULT, &hasMult));
389: if (hasMult) {
390: /* product is available, use it */
391: PetscCall(MatMult(Amat, X, Y));
392: } else {
393: /* there's no product, so treat J0 as identity */
394: PetscCall(VecCopy(X, Y));
395: }
396: } else if (lmvm->user_scale) {
397: if (lmvm->J0diag) {
398: /* User has defined a diagonal vector for J0 */
399: PetscCall(VecPointwiseMult(X, lmvm->J0diag, Y));
400: } else {
401: /* User has defined a scalar value for J0 */
402: PetscCall(VecAXPBY(Y, lmvm->J0scalar, 0.0, X));
403: }
404: } else {
405: /* There is no J0 representation so just apply an identity matrix */
406: PetscCall(VecCopy(X, Y));
407: }
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*@
412: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
413: inverse to the given vector.
415: Input Parameters:
416: + B - A `MATLMVM` matrix
417: - X - vector to "multiply" with J0^{-1}
419: Output Parameter:
420: . Y - resulting vector for the operation
422: Level: advanced
424: Note:
425: The specific form of the application
426: depends on whether the user provided a scaling factor, a J0 matrix,
427: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
428: provided, the function simply does an identity matrix application
429: (vector copy).
431: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
432: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
433: @*/
434: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
435: {
436: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
437: PetscBool same, hasSolve;
439: PetscFunctionBegin;
443: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
444: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
445: lmvm = (Mat_LMVM *)B->data;
446: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
447: VecCheckMatCompatible(B, X, 2, Y, 3);
449: /* Invert the initial Jacobian onto q (or apply scaling) */
450: if (lmvm->user_pc) {
451: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
452: PetscCall(PCApply(lmvm->J0pc, X, Y));
453: } else if (lmvm->user_ksp) {
454: /* User has defined a J0 or a custom KSP so just perform a solution */
455: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
456: } else if (lmvm->J0) {
457: PetscCall(MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve));
458: if (hasSolve) {
459: PetscCall(MatSolve(lmvm->J0, X, Y));
460: } else {
461: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
462: }
463: } else if (lmvm->user_scale) {
464: if (lmvm->J0diag) {
465: PetscCall(VecPointwiseDivide(X, Y, lmvm->J0diag));
466: } else {
467: PetscCall(VecAXPBY(Y, 1.0 / lmvm->J0scalar, 0.0, X));
468: }
469: } else {
470: /* There is no J0 representation so just apply an identity matrix */
471: PetscCall(VecCopy(X, Y));
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: MatLMVMIsAllocated - Returns a boolean flag that shows whether
478: the necessary data structures for the underlying matrix is allocated.
480: Input Parameter:
481: . B - A `MATLMVM` matrix
483: Output Parameter:
484: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
486: Level: intermediate
488: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
489: @*/
490: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
491: {
492: Mat_LMVM *lmvm;
493: PetscBool same;
495: PetscFunctionBegin;
497: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
498: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
499: *flg = PETSC_FALSE;
500: lmvm = (Mat_LMVM *)B->data;
501: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@
506: MatLMVMAllocate - Produces all necessary common memory for
507: LMVM approximations based on the solution and function vectors
508: provided.
510: Input Parameters:
511: + B - A `MATLMVM` matrix
512: . X - Solution vector
513: - F - Function vector
515: Level: intermediate
517: Note:
518: If `MatSetSizes()` and `MatSetUp()` have not been called
519: before `MatLMVMAllocate()`, the allocation will read sizes from
520: the provided vectors and update the matrix.
522: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
523: @*/
524: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
525: {
526: Mat_LMVM *lmvm;
527: PetscBool same;
528: VecType vtype;
530: PetscFunctionBegin;
534: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
535: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
536: lmvm = (Mat_LMVM *)B->data;
537: PetscCall(VecGetType(X, &vtype));
538: PetscCall(MatSetVecType(B, vtype));
539: PetscCall((*lmvm->ops->allocate)(B, X, F));
540: if (lmvm->J0) PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.
547: Input Parameter:
548: . B - A `MATLMVM` matrix
550: Level: intermediate
552: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
553: @*/
554: PetscErrorCode MatLMVMResetShift(Mat B)
555: {
556: Mat_LMVM *lmvm;
557: PetscBool same;
559: PetscFunctionBegin;
561: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
562: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
563: lmvm = (Mat_LMVM *)B->data;
564: lmvm->shift = 0.0;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@
569: MatLMVMReset - Flushes all of the accumulated updates out of
570: the `MATLMVM` approximation.
572: Input Parameters:
573: + B - A `MATLMVM` matrix
574: - destructive - flag for enabling destruction of data structures
576: Level: intermediate
578: Note:
579: In practice, this will not actually
580: destroy the data associated with the updates. It simply resets
581: counters, which leads to existing data being overwritten, and
582: `MatSolve()` being applied as if there are no updates. A boolean
583: flag is available to force destruction of the update vectors.
585: If the user has provided another LMVM matrix as J0, the J0
586: matrix is also reset in this function.
588: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
589: @*/
590: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
591: {
592: Mat_LMVM *lmvm;
593: PetscBool same;
595: PetscFunctionBegin;
597: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
598: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
599: lmvm = (Mat_LMVM *)B->data;
600: PetscCall((*lmvm->ops->reset)(B, destructive));
601: if (lmvm->J0) PetscCall(MatLMVMReset(lmvm->J0, destructive));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@
606: MatLMVMSetHistorySize - Set the number of past iterates to be
607: stored for the construction of the limited-memory quasi-Newton update.
609: Input Parameters:
610: + B - A `MATLMVM` matrix
611: - hist_size - number of past iterates (default 5)
613: Options Database Key:
614: . -mat_lmvm_hist_size <m> - set number of past iterates
616: Level: beginner
618: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
619: @*/
620: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
621: {
622: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
623: PetscBool same;
625: PetscFunctionBegin;
627: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
628: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
629: PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
630: {
631: PetscBool reallocate = PETSC_FALSE;
632: Vec X = NULL, F = NULL;
633: //lmvm->m = hist_size;
634: if (lmvm->allocated && hist_size != lmvm->m) {
635: PetscCall(VecDuplicate(lmvm->Xprev, &X));
636: PetscCall(VecDuplicate(lmvm->Fprev, &F));
637: PetscCall(MatLMVMReset(B, PETSC_TRUE));
638: reallocate = PETSC_TRUE;
639: }
640: lmvm->m = hist_size;
641: if (reallocate) {
642: PetscCall(MatLMVMAllocate(B, X, F));
643: PetscCall(VecDestroy(&X));
644: PetscCall(VecDestroy(&F));
645: }
646: }
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode MatLMVMGetHistorySize(Mat B, PetscInt *hist_size)
651: {
652: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
653: PetscBool same;
655: PetscFunctionBegin;
657: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
658: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
659: *hist_size = lmvm->m;
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@
664: MatLMVMGetUpdateCount - Returns the number of accepted updates.
666: Input Parameter:
667: . B - A `MATLMVM` matrix
669: Output Parameter:
670: . nupdates - number of accepted updates
672: Level: intermediate
674: Note:
675: This number may be greater than the total number of update vectors
676: stored in the matrix. The counters are reset when `MatLMVMReset()`
677: is called.
679: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
680: @*/
681: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
682: {
683: Mat_LMVM *lmvm;
684: PetscBool same;
686: PetscFunctionBegin;
688: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
689: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
690: lmvm = (Mat_LMVM *)B->data;
691: *nupdates = lmvm->nupdates;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: MatLMVMGetRejectCount - Returns the number of rejected updates.
697: The counters are reset when `MatLMVMReset()` is called.
699: Input Parameter:
700: . B - A `MATLMVM` matrix
702: Output Parameter:
703: . nrejects - number of rejected updates
705: Level: intermediate
707: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
708: @*/
709: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
710: {
711: Mat_LMVM *lmvm;
712: PetscBool same;
714: PetscFunctionBegin;
716: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
717: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
718: lmvm = (Mat_LMVM *)B->data;
719: *nrejects = lmvm->nrejects;
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }