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: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
446: VecCheckMatCompatible(B, X, 2, Y, 3);
448: /* Invert the initial Jacobian onto q (or apply scaling) */
449: if (lmvm->user_pc) {
450: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
451: PetscCall(PCApply(lmvm->J0pc, X, Y));
452: } else if (lmvm->user_ksp) {
453: /* User has defined a J0 or a custom KSP so just perform a solution */
454: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
455: } else if (lmvm->J0) {
456: PetscCall(MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve));
457: if (hasSolve) {
458: PetscCall(MatSolve(lmvm->J0, X, Y));
459: } else {
460: PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
461: }
462: } else if (lmvm->user_scale) {
463: if (lmvm->J0diag) {
464: PetscCall(VecPointwiseDivide(X, Y, lmvm->J0diag));
465: } else {
466: PetscCall(VecAXPBY(Y, 1.0 / lmvm->J0scalar, 0.0, X));
467: }
468: } else {
469: /* There is no J0 representation so just apply an identity matrix */
470: PetscCall(VecCopy(X, Y));
471: }
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /*@
476: MatLMVMIsAllocated - Returns a boolean flag that shows whether
477: the necessary data structures for the underlying matrix is allocated.
479: Input Parameter:
480: . B - A `MATLMVM` matrix
482: Output Parameter:
483: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
485: Level: intermediate
487: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
488: @*/
489: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
490: {
491: Mat_LMVM *lmvm;
492: PetscBool same;
494: PetscFunctionBegin;
496: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
497: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
498: *flg = PETSC_FALSE;
499: lmvm = (Mat_LMVM *)B->data;
500: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: /*@
505: MatLMVMAllocate - Produces all necessary common memory for
506: LMVM approximations based on the solution and function vectors
507: provided.
509: Input Parameters:
510: + B - A `MATLMVM` matrix
511: . X - Solution vector
512: - F - Function vector
514: Level: intermediate
516: Note:
517: If `MatSetSizes()` and `MatSetUp()` have not been called
518: before `MatLMVMAllocate()`, the allocation will read sizes from
519: the provided vectors and update the matrix.
521: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
522: @*/
523: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
524: {
525: Mat_LMVM *lmvm;
526: PetscBool same;
527: VecType vtype;
529: PetscFunctionBegin;
533: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
534: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
535: lmvm = (Mat_LMVM *)B->data;
536: PetscCall(VecGetType(X, &vtype));
537: PetscCall(MatSetVecType(B, vtype));
538: PetscCall((*lmvm->ops->allocate)(B, X, F));
539: if (lmvm->J0) PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: /*@
544: MatLMVMResetShift - Zero the shift factor for a `MATLMVM`.
546: Input Parameter:
547: . B - A `MATLMVM` matrix
549: Level: intermediate
551: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
552: @*/
553: PetscErrorCode MatLMVMResetShift(Mat B)
554: {
555: Mat_LMVM *lmvm;
556: PetscBool same;
558: PetscFunctionBegin;
560: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
561: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
562: lmvm = (Mat_LMVM *)B->data;
563: lmvm->shift = 0.0;
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@
568: MatLMVMReset - Flushes all of the accumulated updates out of
569: the `MATLMVM` approximation.
571: Input Parameters:
572: + B - A `MATLMVM` matrix
573: - destructive - flag for enabling destruction of data structures
575: Level: intermediate
577: Note:
578: In practice, this will not actually
579: destroy the data associated with the updates. It simply resets
580: counters, which leads to existing data being overwritten, and
581: `MatSolve()` being applied as if there are no updates. A boolean
582: flag is available to force destruction of the update vectors.
584: If the user has provided another LMVM matrix as J0, the J0
585: matrix is also reset in this function.
587: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
588: @*/
589: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
590: {
591: Mat_LMVM *lmvm;
592: PetscBool same;
594: PetscFunctionBegin;
596: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
597: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
598: lmvm = (Mat_LMVM *)B->data;
599: PetscCall((*lmvm->ops->reset)(B, destructive));
600: if (lmvm->J0) PetscCall(MatLMVMReset(lmvm->J0, destructive));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@
605: MatLMVMSetHistorySize - Set the number of past iterates to be
606: stored for the construction of the limited-memory quasi-Newton update.
608: Input Parameters:
609: + B - A `MATLMVM` matrix
610: - hist_size - number of past iterates (default 5)
612: Options Database Key:
613: . -mat_lmvm_hist_size <m> - set number of past iterates
615: Level: beginner
617: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
618: @*/
619: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
620: {
621: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
622: PetscBool same;
624: PetscFunctionBegin;
626: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
627: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
628: PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
629: {
630: PetscBool reallocate = PETSC_FALSE;
631: Vec X = NULL, F = NULL;
632: //lmvm->m = hist_size;
633: if (lmvm->allocated && hist_size != lmvm->m) {
634: PetscCall(VecDuplicate(lmvm->Xprev, &X));
635: PetscCall(VecDuplicate(lmvm->Fprev, &F));
636: PetscCall(MatLMVMReset(B, PETSC_TRUE));
637: reallocate = PETSC_TRUE;
638: }
639: lmvm->m = hist_size;
640: if (reallocate) {
641: PetscCall(MatLMVMAllocate(B, X, F));
642: PetscCall(VecDestroy(&X));
643: PetscCall(VecDestroy(&F));
644: }
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: PetscErrorCode MatLMVMGetHistorySize(Mat B, PetscInt *hist_size)
650: {
651: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
652: PetscBool same;
654: PetscFunctionBegin;
656: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
657: if (!same) PetscFunctionReturn(PETSC_SUCCESS);
658: *hist_size = lmvm->m;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: MatLMVMGetUpdateCount - Returns the number of accepted updates.
665: Input Parameter:
666: . B - A `MATLMVM` matrix
668: Output Parameter:
669: . nupdates - number of accepted updates
671: Level: intermediate
673: Note:
674: This number may be greater than the total number of update vectors
675: stored in the matrix. The counters are reset when `MatLMVMReset()`
676: is called.
678: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
679: @*/
680: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
681: {
682: Mat_LMVM *lmvm;
683: PetscBool same;
685: PetscFunctionBegin;
687: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
688: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
689: lmvm = (Mat_LMVM *)B->data;
690: *nupdates = lmvm->nupdates;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: MatLMVMGetRejectCount - Returns the number of rejected updates.
696: The counters are reset when `MatLMVMReset()` is called.
698: Input Parameter:
699: . B - A `MATLMVM` matrix
701: Output Parameter:
702: . nrejects - number of rejected updates
704: Level: intermediate
706: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`
707: @*/
708: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
709: {
710: Mat_LMVM *lmvm;
711: PetscBool same;
713: PetscFunctionBegin;
715: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
716: PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
717: lmvm = (Mat_LMVM *)B->data;
718: *nrejects = lmvm->nrejects;
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }