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: [](chapter_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;
29: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
31: if (!lmvm->allocated) {
32: MatLMVMAllocate(B, X, F);
33: } else {
34: VecCheckMatCompatible(B, X, 2, F, 3);
35: }
36: if (lmvm->J0) {
37: /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
38: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
39: if (same) MatLMVMUpdate(lmvm->J0, X, F);
40: }
41: (*lmvm->ops->update)(B, X, F);
42: return 0;
43: }
45: /*@
46: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
47: an identity matrix (scale = 1.0).
49: Input Parameters:
50: . B - An LMVM-type matrix
52: Level: advanced
54: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
55: @*/
56: PetscErrorCode MatLMVMClearJ0(Mat B)
57: {
58: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
59: PetscBool same;
60: MPI_Comm comm = PetscObjectComm((PetscObject)B);
63: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
65: lmvm->user_pc = PETSC_FALSE;
66: lmvm->user_ksp = PETSC_FALSE;
67: lmvm->user_scale = PETSC_FALSE;
68: lmvm->J0scalar = 1.0;
69: VecDestroy(&lmvm->J0diag);
70: MatDestroy(&lmvm->J0);
71: PCDestroy(&lmvm->J0pc);
72: return 0;
73: }
75: /*@
76: MatLMVMSetJ0Scale - Allows the user to define a scalar value
77: mu such that J0 = mu*I.
79: Input Parameters:
80: + B - An LMVM-type matrix
81: - scale - Scalar value mu that defines the initial Jacobian
83: Level: advanced
85: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
86: @*/
87: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
88: {
89: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
90: PetscBool same;
91: MPI_Comm comm = PetscObjectComm((PetscObject)B);
94: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
97: MatLMVMClearJ0(B);
98: lmvm->J0scalar = scale;
99: lmvm->user_scale = PETSC_TRUE;
100: return 0;
101: }
103: /*@
104: MatLMVMSetJ0Diag - Allows the user to define a vector
105: V such that J0 = diag(V).
107: Input Parameters:
108: + B - An LMVM-type matrix
109: - V - Vector that defines the diagonal of the initial Jacobian
111: Level: advanced
113: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
114: @*/
115: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
116: {
117: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
118: PetscBool same;
119: MPI_Comm comm = PetscObjectComm((PetscObject)B);
123: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
127: VecCheckSameSize(V, 2, lmvm->Fprev, 3);
129: MatLMVMClearJ0(B);
130: if (!lmvm->J0diag) VecDuplicate(V, &lmvm->J0diag);
131: VecCopy(V, lmvm->J0diag);
132: lmvm->user_scale = PETSC_TRUE;
133: return 0;
134: }
136: /*@
137: MatLMVMSetJ0 - Allows the user to define the initial
138: Jacobian matrix from which the LMVM-type approximation is
139: built up. Inverse of this initial Jacobian is applied
140: using an internal `KSP` solver, which defaults to `KSPGMRES`.
141: This internal `KSP` solver has the "mat_lmvm_" option
142: prefix.
144: Note that another LMVM-type matrix can be used in place of
145: J0, in which case updating the outer LMVM-type matrix will
146: also trigger the update for the inner LMVM-type matrix. This
147: is useful in cases where a full-memory diagonal approximation
148: such as `MATLMVMDIAGBRDN` is used in place of J0.
150: Input Parameters:
151: + B - An LMVM-type matrix
152: - J0 - The initial Jacobian matrix
154: Level: advanced
156: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
157: @*/
158: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
159: {
160: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
161: PetscBool same;
162: MPI_Comm comm = PetscObjectComm((PetscObject)B);
166: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
168: MatLMVMClearJ0(B);
169: MatDestroy(&lmvm->J0);
170: PetscObjectReference((PetscObject)J0);
171: lmvm->J0 = J0;
172: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
173: if (!same && lmvm->square) KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0);
174: return 0;
175: }
177: /*@
178: MatLMVMSetJ0PC - Allows the user to define a `PC` object that
179: acts as the initial inverse-Jacobian matrix. This `PC` should
180: already contain all the operators necessary for its application.
181: The LMVM-type matrix only calls `PCApply()` without changing any other
182: options.
184: Input Parameters:
185: + B - An LMVM-type matrix
186: - J0pc - `PC` object where `PCApply()` defines an inverse application for J0
188: Level: advanced
190: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
191: @*/
192: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
193: {
194: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
195: PetscBool same;
196: MPI_Comm comm = PetscObjectComm((PetscObject)B);
200: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
203: MatLMVMClearJ0(B);
204: PetscObjectReference((PetscObject)J0pc);
205: lmvm->J0pc = J0pc;
206: lmvm->user_pc = PETSC_TRUE;
207: return 0;
208: }
210: /*@
211: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
212: KSP solver for the initial inverse-Jacobian approximation.
213: This `KSP` solver should already contain all the operators
214: necessary to perform the inversion. The LMVM-type matrix only
215: calls `KSPSolve()` without changing any other options.
217: Input Parameters:
218: + B - An LMVM-type matrix
219: - J0ksp - `KSP` solver for the initial inverse-Jacobian application
221: Level: advanced
223: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
224: @*/
225: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
226: {
227: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
228: PetscBool same;
229: MPI_Comm comm = PetscObjectComm((PetscObject)B);
233: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
236: MatLMVMClearJ0(B);
237: KSPDestroy(&lmvm->J0ksp);
238: PetscObjectReference((PetscObject)J0ksp);
239: lmvm->J0ksp = J0ksp;
240: lmvm->user_ksp = PETSC_TRUE;
241: return 0;
242: }
244: /*@
245: MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.
247: Input Parameters:
248: . B - An LMVM-type matrix
250: Output Parameter:
251: . J0 - `Mat` object for defining the initial Jacobian
253: Level: advanced
255: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
256: @*/
257: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
258: {
259: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
260: PetscBool same;
263: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
265: *J0 = lmvm->J0;
266: return 0;
267: }
269: /*@
270: MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
271: associated with the initial Jacobian.
273: Input Parameters:
274: . B - An LMVM-type matrix
276: Output Parameter:
277: . J0pc - `PC` object for defining the initial inverse-Jacobian
279: Level: advanced
281: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
282: @*/
283: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
284: {
285: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
286: PetscBool same;
289: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
291: if (lmvm->J0pc) {
292: *J0pc = lmvm->J0pc;
293: } else {
294: KSPGetPC(lmvm->J0ksp, J0pc);
295: }
296: return 0;
297: }
299: /*@
300: MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
301: associated with the initial Jacobian.
303: Input Parameters:
304: . B - An LMVM-type matrix
306: Output Parameter:
307: . J0ksp - `KSP` solver for defining the initial inverse-Jacobian
309: Level: advanced
311: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
312: @*/
313: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
314: {
315: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
316: PetscBool same;
319: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
321: *J0ksp = lmvm->J0ksp;
322: return 0;
323: }
325: /*@
326: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
327: matrix-vector product with the initial Jacobian.
329: Input Parameters:
330: + B - An LMVM-type matrix
331: - X - vector to multiply with J0
333: Output Parameter:
334: . Y - resulting vector for the operation
336: Level: advanced
338: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
339: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
340: @*/
341: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
342: {
343: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
344: PetscBool same, hasMult;
345: MPI_Comm comm = PetscObjectComm((PetscObject)B);
346: Mat Amat, Pmat;
351: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
354: VecCheckMatCompatible(B, X, 2, Y, 3);
355: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
356: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
357: if (lmvm->user_pc) {
358: PCGetOperators(lmvm->J0pc, &Amat, &Pmat);
359: } else if (lmvm->user_ksp) {
360: KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat);
361: } else {
362: Amat = lmvm->J0;
363: }
364: MatHasOperation(Amat, MATOP_MULT, &hasMult);
365: if (hasMult) {
366: /* product is available, use it */
367: MatMult(Amat, X, Y);
368: } else {
369: /* there's no product, so treat J0 as identity */
370: VecCopy(X, Y);
371: }
372: } else if (lmvm->user_scale) {
373: if (lmvm->J0diag) {
374: /* User has defined a diagonal vector for J0 */
375: VecPointwiseMult(X, lmvm->J0diag, Y);
376: } else {
377: /* User has defined a scalar value for J0 */
378: VecCopy(X, Y);
379: VecScale(Y, lmvm->J0scalar);
380: }
381: } else {
382: /* There is no J0 representation so just apply an identity matrix */
383: VecCopy(X, Y);
384: }
385: return 0;
386: }
388: /*@
389: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
390: inverse to the given vector. The specific form of the application
391: depends on whether the user provided a scaling factor, a J0 matrix,
392: a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
393: provided, the function simply does an identity matrix application
394: (vector copy).
396: Input Parameters:
397: + B - An LMVM-type matrix
398: - X - vector to "multiply" with J0^{-1}
400: Output Parameter:
401: . Y - resulting vector for the operation
403: Level: advanced
405: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
406: `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
407: @*/
408: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
409: {
410: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
411: PetscBool same, hasSolve;
412: MPI_Comm comm = PetscObjectComm((PetscObject)B);
417: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
420: VecCheckMatCompatible(B, X, 2, Y, 3);
422: /* Invert the initial Jacobian onto q (or apply scaling) */
423: if (lmvm->user_pc) {
424: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
425: PCApply(lmvm->J0pc, X, Y);
426: } else if (lmvm->user_ksp) {
427: /* User has defined a J0 or a custom KSP so just perform a solution */
428: KSPSolve(lmvm->J0ksp, X, Y);
429: } else if (lmvm->J0) {
430: MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve);
431: if (hasSolve) {
432: MatSolve(lmvm->J0, X, Y);
433: } else {
434: KSPSolve(lmvm->J0ksp, X, Y);
435: }
436: } else if (lmvm->user_scale) {
437: if (lmvm->J0diag) {
438: VecPointwiseDivide(X, Y, lmvm->J0diag);
439: } else {
440: VecCopy(X, Y);
441: VecScale(Y, 1.0 / lmvm->J0scalar);
442: }
443: } else {
444: /* There is no J0 representation so just apply an identity matrix */
445: VecCopy(X, Y);
446: }
447: return 0;
448: }
450: /*@
451: MatLMVMIsAllocated - Returns a boolean flag that shows whether
452: the necessary data structures for the underlying matrix is allocated.
454: Input Parameters:
455: . B - An LMVM-type matrix
457: Output Parameter:
458: . flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise
460: Level: intermediate
462: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
463: @*/
465: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
466: {
467: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
468: PetscBool same;
471: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
473: *flg = PETSC_FALSE;
474: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
475: return 0;
476: }
478: /*@
479: MatLMVMAllocate - Produces all necessary common memory for
480: LMVM approximations based on the solution and function vectors
481: provided. If `MatSetSizes()` and `MatSetUp()` have not been called
482: before `MatLMVMAllocate()`, the allocation will read sizes from
483: the provided vectors and update the matrix.
485: Input Parameters:
486: + B - An LMVM-type matrix
487: . X - Solution vector
488: - F - Function vector
490: Level: intermediate
492: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
493: @*/
494: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
495: {
496: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
497: PetscBool same;
502: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
504: (*lmvm->ops->allocate)(B, X, F);
505: if (lmvm->J0) {
506: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
507: if (same) MatLMVMAllocate(lmvm->J0, X, F);
508: }
509: return 0;
510: }
512: /*@
513: MatLMVMResetShift - Zero the shift factor.
515: Input Parameters:
516: . B - An LMVM-type matrix
518: Level: intermediate
520: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
521: @*/
522: PetscErrorCode MatLMVMResetShift(Mat B)
523: {
524: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
525: PetscBool same;
528: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
530: lmvm->shift = 0.0;
531: return 0;
532: }
534: /*@
535: MatLMVMReset - Flushes all of the accumulated updates out of
536: the LMVM approximation. In practice, this will not actually
537: destroy the data associated with the updates. It simply resets
538: counters, which leads to existing data being overwritten, and
539: `MatSolve()` being applied as if there are no updates. A boolean
540: flag is available to force destruction of the update vectors.
542: If the user has provided another LMVM matrix as J0, the J0
543: matrix is also reset in this function.
545: Input Parameters:
546: + B - An LMVM-type matrix
547: - destructive - flag for enabling destruction of data structures
549: Level: intermediate
551: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
552: @*/
553: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
554: {
555: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
556: PetscBool same;
559: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
561: (*lmvm->ops->reset)(B, destructive);
562: if (lmvm->J0) {
563: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
564: if (same) MatLMVMReset(lmvm->J0, destructive);
565: }
566: return 0;
567: }
569: /*@
570: MatLMVMSetHistorySize - Set the number of past iterates to be
571: stored for the construction of the limited-memory QN update.
573: Input Parameters:
574: + B - An LMVM-type matrix
575: - hist_size - number of past iterates (default 5)
577: Options Database Key:
578: . -mat_lmvm_hist_size <m> - set number of past iterates
580: Level: beginner
582: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
583: @*/
584: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
585: {
586: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
587: PetscBool same;
588: Vec X, F;
591: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
593: if (hist_size > 0) {
594: lmvm->m = hist_size;
595: if (lmvm->allocated && lmvm->m != lmvm->m_old) {
596: VecDuplicate(lmvm->Xprev, &X);
597: VecDuplicate(lmvm->Fprev, &F);
598: MatLMVMReset(B, PETSC_TRUE);
599: MatLMVMAllocate(B, X, F);
600: VecDestroy(&X);
601: VecDestroy(&F);
602: }
604: return 0;
605: }
607: /*@
608: MatLMVMGetUpdateCount - Returns the number of accepted updates.
609: This number may be greater than the total number of update vectors
610: stored in the matrix. The counters are reset when `MatLMVMReset()`
611: is called.
613: Input Parameters:
614: . B - An LMVM-type matrix
616: Output Parameter:
617: . nupdates - number of accepted updates
619: Level: intermediate
621: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
622: @*/
623: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
624: {
625: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
626: PetscBool same;
629: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
631: *nupdates = lmvm->nupdates;
632: return 0;
633: }
635: /*@
636: MatLMVMGetRejectCount - Returns the number of rejected updates.
637: The counters are reset when `MatLMVMReset()` is called.
639: Input Parameters:
640: . B - An LMVM-type matrix
642: Output Parameter:
643: . nrejects - number of rejected updates
645: Level: intermediate
647: .seealso: [](chapter_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
648: @*/
649: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
650: {
651: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
652: PetscBool same;
655: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
657: *nrejects = lmvm->nrejects;
658: return 0;
659: }