Actual source code: lmvmutils.c
petsc-3.14.6 2021-03-30
1: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: /*@
4: MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to an LMVM matrix.
5: The first time the function is called for an LMVM 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 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: MatLMVMReset(), MatLMVMAllocate()
20: @*/
21: PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
22: {
23: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
24: PetscErrorCode ierr;
25: PetscBool same;
31: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
32: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
33: if (!lmvm->allocated) {
34: MatLMVMAllocate(B, X, F);
35: } else {
36: VecCheckMatCompatible(B, X, 2, F, 3);
37: }
38: if (lmvm->J0) {
39: /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
40: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
41: if (same) {
42: MatLMVMUpdate(lmvm->J0, X, F);
43: }
44: }
45: (*lmvm->ops->update)(B, X, F);
46: return(0);
47: }
49: /*------------------------------------------------------------*/
51: /*@
52: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
53: an identity matrix (scale = 1.0).
55: Input Parameters:
56: . B - An LMVM-type matrix
58: Level: advanced
60: .seealso: MatLMVMSetJ0()
61: @*/
62: PetscErrorCode MatLMVMClearJ0(Mat B)
63: {
64: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
65: PetscErrorCode ierr;
66: PetscBool same;
67: MPI_Comm comm = PetscObjectComm((PetscObject)B);
71: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
72: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
73: lmvm->user_pc = PETSC_FALSE;
74: lmvm->user_ksp = PETSC_FALSE;
75: lmvm->user_scale = PETSC_FALSE;
76: lmvm->J0scalar = 1.0;
77: VecDestroy(&lmvm->J0diag);
78: MatDestroy(&lmvm->J0);
79: PCDestroy(&lmvm->J0pc);
80: return(0);
81: }
83: /*------------------------------------------------------------*/
85: /*@
86: MatLMVMSetJ0Scale - Allows the user to define a scalar value
87: mu such that J0 = mu*I.
89: Input Parameters:
90: + B - An LMVM-type matrix
91: - scale - Scalar value mu that defines the initial Jacobian
93: Level: advanced
95: .seealso: MatLMVMSetDiagScale(), MatLMVMSetJ0()
96: @*/
97: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
98: {
99: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
100: PetscErrorCode ierr;
101: PetscBool same;
102: MPI_Comm comm = PetscObjectComm((PetscObject)B);
106: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
107: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
108: if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Scaling is available only for square LMVM matrices");
109: MatLMVMClearJ0(B);
110: lmvm->J0scalar = scale;
111: lmvm->user_scale = PETSC_TRUE;
112: return(0);
113: }
115: /*------------------------------------------------------------*/
117: /*@
118: MatLMVMSetJ0Diag - Allows the user to define a vector
119: V such that J0 = diag(V).
121: Input Parameters:
122: + B - An LMVM-type matrix
123: - V - Vector that defines the diagonal of the initial Jacobian
125: Level: advanced
127: .seealso: MatLMVMSetScale(), MatLMVMSetJ0()
128: @*/
129: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
130: {
131: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
132: PetscErrorCode ierr;
133: PetscBool same;
134: MPI_Comm comm = PetscObjectComm((PetscObject)B);
139: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
140: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
141: if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "Matrix must be allocated before setting diagonal scaling");
142: if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Diagonal scaling is available only for square LMVM matrices");
143: VecCheckSameSize(V, 2, lmvm->Fprev, 3);
144: MatLMVMClearJ0(B);
145: if (!lmvm->J0diag) {
146: VecDuplicate(V, &lmvm->J0diag);
147: }
148: VecCopy(V, lmvm->J0diag);
149: lmvm->user_scale = PETSC_TRUE;
150: return(0);
151: }
153: /*------------------------------------------------------------*/
155: /*@
156: MatLMVMSetJ0 - Allows the user to define the initial
157: Jacobian matrix from which the LMVM approximation is
158: built up. Inverse of this initial Jacobian is applied
159: using an internal KSP solver, which defaults to GMRES.
160: This internal KSP solver has the "mat_lmvm_" option
161: prefix.
163: Note that another LMVM matrix can be used in place of
164: J0, in which case updating the outer LMVM matrix will
165: also trigger the update for the inner LMVM matrix. This
166: is useful in cases where a full-memory diagonal approximation
167: such as MATLMVMDIAGBRDN is used in place of J0.
169: Input Parameters:
170: + B - An LMVM-type matrix
171: - J0 - The initial Jacobian matrix
173: Level: advanced
175: .seealso: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP()
176: @*/
177: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
178: {
179: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
180: PetscErrorCode ierr;
181: PetscBool same;
182: MPI_Comm comm = PetscObjectComm((PetscObject)B);
187: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
188: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
189: MatLMVMClearJ0(B);
190: MatDestroy(&lmvm->J0);
191: PetscObjectReference((PetscObject)J0);
192: lmvm->J0 = J0;
193: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
194: if (!same && lmvm->square) {
195: KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0);
196: }
197: return(0);
198: }
200: /*------------------------------------------------------------*/
202: /*@
203: MatLMVMSetJ0PC - Allows the user to define a PC object that
204: acts as the initial inverse-Jacobian matrix. This PC should
205: already contain all the operators necessary for its application.
206: The LMVM matrix only calls PCApply() without changing any other
207: options.
209: Input Parameters:
210: + B - An LMVM-type matrix
211: - J0pc - PC object where PCApply defines an inverse application for J0
213: Level: advanced
215: .seealso: MatLMVMGetJ0PC()
216: @*/
217: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
218: {
219: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
220: PetscErrorCode ierr;
221: PetscBool same;
222: MPI_Comm comm = PetscObjectComm((PetscObject)B);
227: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
228: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
229: if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
230: MatLMVMClearJ0(B);
231: PetscObjectReference((PetscObject)J0pc);
232: lmvm->J0pc = J0pc;
233: lmvm->user_pc = PETSC_TRUE;
234: return(0);
235: }
237: /*------------------------------------------------------------*/
239: /*@
240: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
241: KSP solver for the initial inverse-Jacobian approximation.
242: This KSP solver should already contain all the operators
243: necessary to perform the inversion. The LMVM matrix only
244: calls KSPSolve() without changing any other options.
246: Input Parameters:
247: + B - An LMVM-type matrix
248: - J0ksp - KSP solver for the initial inverse-Jacobian application
250: Level: advanced
252: .seealso: MatLMVMGetJ0KSP()
253: @*/
254: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
255: {
256: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
257: PetscErrorCode ierr;
258: PetscBool same;
259: MPI_Comm comm = PetscObjectComm((PetscObject)B);
264: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
265: if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
266: if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
267: MatLMVMClearJ0(B);
268: KSPDestroy(&lmvm->J0ksp);
269: PetscObjectReference((PetscObject)J0ksp);
270: lmvm->J0ksp = J0ksp;
271: lmvm->user_ksp = PETSC_TRUE;
272: return(0);
273: }
275: /*------------------------------------------------------------*/
277: /*@
278: MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.
280: Input Parameters:
281: . B - An LMVM-type matrix
283: Output Parameter:
284: . J0 - Mat object for defining the initial Jacobian
286: Level: advanced
288: .seealso: MatLMVMSetJ0()
289: @*/
290: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
291: {
292: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
293: PetscErrorCode ierr;
294: PetscBool same;
298: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
299: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
300: *J0 = lmvm->J0;
301: return(0);
302: }
304: /*------------------------------------------------------------*/
306: /*@
307: MatLMVMGetJ0PC - Returns a pointer to the internal PC object
308: associated with the initial Jacobian.
310: Input Parameters:
311: . B - An LMVM-type matrix
313: Output Parameter:
314: . J0pc - PC object for defining the initial inverse-Jacobian
316: Level: advanced
318: .seealso: MatLMVMSetJ0PC()
319: @*/
320: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
321: {
322: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
323: PetscErrorCode ierr;
324: PetscBool same;
328: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
329: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
330: if (lmvm->J0pc) {
331: *J0pc = lmvm->J0pc;
332: } else {
333: KSPGetPC(lmvm->J0ksp, J0pc);
334: }
335: return(0);
336: }
338: /*------------------------------------------------------------*/
340: /*@
341: MatLMVMGetJ0KSP - Returns a pointer to the internal KSP solver
342: associated with the initial Jacobian.
344: Input Parameters:
345: . B - An LMVM-type matrix
347: Output Parameter:
348: . J0ksp - KSP solver for defining the initial inverse-Jacobian
350: Level: advanced
352: .seealso: MatLMVMSetJ0KSP()
353: @*/
354: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
355: {
356: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
357: PetscErrorCode ierr;
358: PetscBool same;
362: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
363: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
364: *J0ksp = lmvm->J0ksp;
365: return(0);
366: }
368: /*------------------------------------------------------------*/
370: /*@
371: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
372: matrix-vector product with the initial Jacobian.
374: Input Parameters:
375: + B - An LMVM-type matrix
376: - X - vector to multiply with J0
378: Output Parameter:
379: . Y - resulting vector for the operation
381: Level: advanced
383: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
384: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Inv()
385: @*/
386: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
387: {
388: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
389: PetscErrorCode ierr;
390: PetscBool same, hasMult;
391: MPI_Comm comm = PetscObjectComm((PetscObject)B);
392: Mat Amat, Pmat;
398: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
399: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
400: if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
401: VecCheckMatCompatible(B, X, 2, Y, 3);
402: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
403: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
404: if (lmvm->user_pc){
405: PCGetOperators(lmvm->J0pc, &Amat, &Pmat);
406: } else if (lmvm->user_ksp) {
407: KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat);
408: } else {
409: Amat = lmvm->J0;
410: }
411: MatHasOperation(Amat, MATOP_MULT, &hasMult);
412: if (hasMult) {
413: /* product is available, use it */
414: MatMult(Amat, X, Y);
415: } else {
416: /* there's no product, so treat J0 as identity */
417: VecCopy(X, Y);
418: }
419: } else if (lmvm->user_scale) {
420: if (lmvm->J0diag) {
421: /* User has defined a diagonal vector for J0 */
422: VecPointwiseMult(X, lmvm->J0diag, Y);
423: } else {
424: /* User has defined a scalar value for J0 */
425: VecCopy(X, Y);
426: VecScale(Y, lmvm->J0scalar);
427: }
428: } else {
429: /* There is no J0 representation so just apply an identity matrix */
430: VecCopy(X, Y);
431: }
432: return(0);
433: }
435: /*------------------------------------------------------------*/
437: /*@
438: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
439: inverse to the given vector. The specific form of the application
440: depends on whether the user provided a scaling factor, a J0 matrix,
441: a J0 PC, or a J0 KSP object. If no form of the initial Jacobian is
442: provided, the function simply does an identity matrix application
443: (vector copy).
445: Input Parameters:
446: + B - An LMVM-type matrix
447: - X - vector to "multiply" with J0^{-1}
449: Output Parameter:
450: . Y - resulting vector for the operation
452: Level: advanced
454: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
455: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Fwd()
456: @*/
457: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
458: {
459: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
460: PetscErrorCode ierr;
461: PetscBool same, hasSolve;
462: MPI_Comm comm = PetscObjectComm((PetscObject)B);
468: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
469: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
470: if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
471: VecCheckMatCompatible(B, X, 2, Y, 3);
472: /* Invert the initial Jacobian onto q (or apply scaling) */
473: if (lmvm->user_pc) {
474: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
475: PCApply(lmvm->J0pc, X, Y);
476: } else if (lmvm->user_ksp) {
477: /* User has defined a J0 or a custom KSP so just perform a solution */
478: KSPSolve(lmvm->J0ksp, X, Y);
479: } else if (lmvm->J0) {
480: MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve);
481: if (hasSolve) {
482: MatSolve(lmvm->J0, X, Y);
483: } else {
484: KSPSolve(lmvm->J0ksp, X, Y);
485: }
486: } else if (lmvm->user_scale) {
487: if (lmvm->J0diag) {
488: VecPointwiseDivide(X, Y, lmvm->J0diag);
489: } else {
490: VecCopy(X, Y);
491: VecScale(Y, 1.0/lmvm->J0scalar);
492: }
493: } else {
494: /* There is no J0 representation so just apply an identity matrix */
495: VecCopy(X, Y);
496: }
497: return(0);
498: }
500: /*------------------------------------------------------------*/
502: /*@
503: MatLMVMIsAllocated - Returns a boolean flag that shows whether
504: the necessary data structures for the underlying matrix is allocated.
506: Input Parameters:
507: . B - An LMVM-type matrix
509: Output Parameter:
510: . flg - PETSC_TRUE if allocated, PETSC_FALSE otherwise
512: Level: intermediate
514: .seealso: MatLMVMAllocate(), MatLMVMReset()
515: @*/
517: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
518: {
519: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
520: PetscErrorCode ierr;
521: PetscBool same;
525: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
526: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
527: *flg = PETSC_FALSE;
528: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
529: return(0);
530: }
532: /*------------------------------------------------------------*/
534: /*@
535: MatLMVMAllocate - Produces all necessary common memory for
536: LMVM approximations based on the solution and function vectors
537: provided. If MatSetSizes() and MatSetUp() have not been called
538: before MatLMVMAllocate(), the allocation will read sizes from
539: the provided vectors and update the matrix.
541: Input Parameters:
542: + B - An LMVM-type matrix
543: . X - Solution vector
544: - F - Function vector
546: Level: intermediate
548: .seealso: MatLMVMReset(), MatLMVMUpdate()
549: @*/
550: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
551: {
552: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
553: PetscErrorCode ierr;
554: PetscBool same;
560: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
561: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
562: (*lmvm->ops->allocate)(B, X, F);
563: if (lmvm->J0) {
564: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
565: if (same) {
566: MatLMVMAllocate(lmvm->J0, X, F);
567: }
568: }
569: return(0);
570: }
572: /*------------------------------------------------------------*/
574: /*@
575: MatLMVMResetShift - Zero the shift factor.
577: Input Parameters:
578: . B - An LMVM-type matrix
580: Level: intermediate
582: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
583: @*/
584: PetscErrorCode MatLMVMResetShift(Mat B)
585: {
586: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
587: PetscErrorCode ierr;
588: PetscBool same;
592: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
593: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
594: lmvm->shift = 0.0;
595: return(0);
596: }
598: /*------------------------------------------------------------*/
600: /*@
601: MatLMVMReset - Flushes all of the accumulated updates out of
602: the LMVM approximation. In practice, this will not actually
603: destroy the data associated with the updates. It simply resets
604: counters, which leads to existing data being overwritten, and
605: MatSolve() being applied as if there are no updates. A boolean
606: flag is available to force destruction of the update vectors.
608: If the user has provided another LMVM matrix as J0, the J0
609: matrix is also reset in this function.
611: Input Parameters:
612: + B - An LMVM-type matrix
613: - destructive - flag for enabling destruction of data structures
615: Level: intermediate
617: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
618: @*/
619: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
620: {
621: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
622: PetscErrorCode ierr;
623: PetscBool same;
627: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
628: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
629: (*lmvm->ops->reset)(B, destructive);
630: if (lmvm->J0) {
631: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
632: if (same) {
633: MatLMVMReset(lmvm->J0, destructive);
634: }
635: }
636: return(0);
637: }
639: /*------------------------------------------------------------*/
641: /*@
642: MatLMVMSetHistorySize - Set the number of past iterates to be
643: stored for the construction of the limited-memory QN update.
645: Input Parameters:
646: + B - An LMVM-type matrix
647: - hist_size - number of past iterates (default 5)
649: Options Database:
650: . -mat_lmvm_hist_size <m>
652: Level: beginner
654: .seealso: MatLMVMGetUpdateCount()
655: @*/
657: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
658: {
659: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
660: PetscErrorCode ierr;
661: PetscBool same;
662: Vec X, F;
666: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
667: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
668: if (hist_size > 0) {
669: lmvm->m = hist_size;
670: if (lmvm->allocated && lmvm->m != lmvm->m_old) {
671: VecDuplicate(lmvm->Xprev, &X);
672: VecDuplicate(lmvm->Fprev, &F);
673: MatLMVMReset(B, PETSC_TRUE);
674: MatLMVMAllocate(B, X, F);
675: VecDestroy(&X);
676: VecDestroy(&F);
677: }
678: } else {
679: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a positive integer.");
680: }
682: return(0);
683: }
685: /*------------------------------------------------------------*/
687: /*@
688: MatLMVMGetUpdateCount - Returns the number of accepted updates.
689: This number may be greater than the total number of update vectors
690: stored in the matrix. The counters are reset when MatLMVMReset()
691: is called.
693: Input Parameters:
694: . B - An LMVM-type matrix
696: Output Parameter:
697: . nupdates - number of accepted updates
699: Level: intermediate
701: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
702: @*/
703: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
704: {
705: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
706: PetscErrorCode ierr;
707: PetscBool same;
711: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
712: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
713: *nupdates = lmvm->nupdates;
714: return(0);
715: }
717: /*------------------------------------------------------------*/
719: /*@
720: MatLMVMGetRejectCount - Returns the number of rejected updates.
721: The counters are reset when MatLMVMReset() is called.
723: Input Parameters:
724: . B - An LMVM-type matrix
726: Output Parameter:
727: . nrejects - number of rejected updates
729: Level: intermediate
731: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
732: @*/
733: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
734: {
735: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
736: PetscErrorCode ierr;
737: PetscBool same;
741: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
742: if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
743: *nrejects = lmvm->nrejects;
744: return(0);
745: }