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 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: 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) {
40: MatLMVMUpdate(lmvm->J0, X, F);
41: }
42: }
43: (*lmvm->ops->update)(B, X, F);
44: return 0;
45: }
47: /*------------------------------------------------------------*/
49: /*@
50: MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
51: an identity matrix (scale = 1.0).
53: Input Parameters:
54: . B - An LMVM-type matrix
56: Level: advanced
58: .seealso: MatLMVMSetJ0()
59: @*/
60: PetscErrorCode MatLMVMClearJ0(Mat B)
61: {
62: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
63: PetscBool same;
64: MPI_Comm comm = PetscObjectComm((PetscObject)B);
67: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
69: lmvm->user_pc = PETSC_FALSE;
70: lmvm->user_ksp = PETSC_FALSE;
71: lmvm->user_scale = PETSC_FALSE;
72: lmvm->J0scalar = 1.0;
73: VecDestroy(&lmvm->J0diag);
74: MatDestroy(&lmvm->J0);
75: PCDestroy(&lmvm->J0pc);
76: return 0;
77: }
79: /*------------------------------------------------------------*/
81: /*@
82: MatLMVMSetJ0Scale - Allows the user to define a scalar value
83: mu such that J0 = mu*I.
85: Input Parameters:
86: + B - An LMVM-type matrix
87: - scale - Scalar value mu that defines the initial Jacobian
89: Level: advanced
91: .seealso: MatLMVMSetDiagScale(), MatLMVMSetJ0()
92: @*/
93: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
94: {
95: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
96: PetscBool same;
97: MPI_Comm comm = PetscObjectComm((PetscObject)B);
100: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
103: MatLMVMClearJ0(B);
104: lmvm->J0scalar = scale;
105: lmvm->user_scale = PETSC_TRUE;
106: return 0;
107: }
109: /*------------------------------------------------------------*/
111: /*@
112: MatLMVMSetJ0Diag - Allows the user to define a vector
113: V such that J0 = diag(V).
115: Input Parameters:
116: + B - An LMVM-type matrix
117: - V - Vector that defines the diagonal of the initial Jacobian
119: Level: advanced
121: .seealso: MatLMVMSetScale(), MatLMVMSetJ0()
122: @*/
123: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
124: {
125: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
126: PetscBool same;
127: MPI_Comm comm = PetscObjectComm((PetscObject)B);
131: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
135: VecCheckSameSize(V, 2, lmvm->Fprev, 3);
137: MatLMVMClearJ0(B);
138: if (!lmvm->J0diag) {
139: VecDuplicate(V, &lmvm->J0diag);
140: }
141: VecCopy(V, lmvm->J0diag);
142: lmvm->user_scale = PETSC_TRUE;
143: return 0;
144: }
146: /*------------------------------------------------------------*/
148: /*@
149: MatLMVMSetJ0 - Allows the user to define the initial
150: Jacobian matrix from which the LMVM approximation is
151: built up. Inverse of this initial Jacobian is applied
152: using an internal KSP solver, which defaults to GMRES.
153: This internal KSP solver has the "mat_lmvm_" option
154: prefix.
156: Note that another LMVM matrix can be used in place of
157: J0, in which case updating the outer LMVM matrix will
158: also trigger the update for the inner LMVM matrix. This
159: is useful in cases where a full-memory diagonal approximation
160: such as MATLMVMDIAGBRDN is used in place of J0.
162: Input Parameters:
163: + B - An LMVM-type matrix
164: - J0 - The initial Jacobian matrix
166: Level: advanced
168: .seealso: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP()
169: @*/
170: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
171: {
172: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
173: PetscBool same;
174: MPI_Comm comm = PetscObjectComm((PetscObject)B);
178: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
180: MatLMVMClearJ0(B);
181: MatDestroy(&lmvm->J0);
182: PetscObjectReference((PetscObject)J0);
183: lmvm->J0 = J0;
184: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
185: if (!same && lmvm->square) {
186: KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0);
187: }
188: return 0;
189: }
191: /*------------------------------------------------------------*/
193: /*@
194: MatLMVMSetJ0PC - Allows the user to define a PC object that
195: acts as the initial inverse-Jacobian matrix. This PC should
196: already contain all the operators necessary for its application.
197: The LMVM matrix only calls PCApply() without changing any other
198: options.
200: Input Parameters:
201: + B - An LMVM-type matrix
202: - J0pc - PC object where PCApply defines an inverse application for J0
204: Level: advanced
206: .seealso: MatLMVMGetJ0PC()
207: @*/
208: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
209: {
210: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
211: PetscBool same;
212: MPI_Comm comm = PetscObjectComm((PetscObject)B);
216: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
219: MatLMVMClearJ0(B);
220: PetscObjectReference((PetscObject)J0pc);
221: lmvm->J0pc = J0pc;
222: lmvm->user_pc = PETSC_TRUE;
223: return 0;
224: }
226: /*------------------------------------------------------------*/
228: /*@
229: MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
230: KSP solver for the initial inverse-Jacobian approximation.
231: This KSP solver should already contain all the operators
232: necessary to perform the inversion. The LMVM matrix only
233: calls KSPSolve() without changing any other options.
235: Input Parameters:
236: + B - An LMVM-type matrix
237: - J0ksp - KSP solver for the initial inverse-Jacobian application
239: Level: advanced
241: .seealso: MatLMVMGetJ0KSP()
242: @*/
243: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
244: {
245: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
246: PetscBool same;
247: MPI_Comm comm = PetscObjectComm((PetscObject)B);
251: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
254: MatLMVMClearJ0(B);
255: KSPDestroy(&lmvm->J0ksp);
256: PetscObjectReference((PetscObject)J0ksp);
257: lmvm->J0ksp = J0ksp;
258: lmvm->user_ksp = PETSC_TRUE;
259: return 0;
260: }
262: /*------------------------------------------------------------*/
264: /*@
265: MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.
267: Input Parameters:
268: . B - An LMVM-type matrix
270: Output Parameter:
271: . J0 - Mat object for defining the initial Jacobian
273: Level: advanced
275: .seealso: MatLMVMSetJ0()
276: @*/
277: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
278: {
279: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
280: PetscBool same;
283: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
285: *J0 = lmvm->J0;
286: return 0;
287: }
289: /*------------------------------------------------------------*/
291: /*@
292: MatLMVMGetJ0PC - Returns a pointer to the internal PC object
293: associated with the initial Jacobian.
295: Input Parameters:
296: . B - An LMVM-type matrix
298: Output Parameter:
299: . J0pc - PC object for defining the initial inverse-Jacobian
301: Level: advanced
303: .seealso: MatLMVMSetJ0PC()
304: @*/
305: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
306: {
307: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
308: PetscBool same;
311: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
313: if (lmvm->J0pc) {
314: *J0pc = lmvm->J0pc;
315: } else {
316: KSPGetPC(lmvm->J0ksp, J0pc);
317: }
318: return 0;
319: }
321: /*------------------------------------------------------------*/
323: /*@
324: MatLMVMGetJ0KSP - Returns a pointer to the internal KSP solver
325: associated with the initial Jacobian.
327: Input Parameters:
328: . B - An LMVM-type matrix
330: Output Parameter:
331: . J0ksp - KSP solver for defining the initial inverse-Jacobian
333: Level: advanced
335: .seealso: MatLMVMSetJ0KSP()
336: @*/
337: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
338: {
339: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
340: PetscBool same;
343: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
345: *J0ksp = lmvm->J0ksp;
346: return 0;
347: }
349: /*------------------------------------------------------------*/
351: /*@
352: MatLMVMApplyJ0Fwd - Applies an approximation of the forward
353: matrix-vector product with the initial Jacobian.
355: Input Parameters:
356: + B - An LMVM-type matrix
357: - X - vector to multiply with J0
359: Output Parameter:
360: . Y - resulting vector for the operation
362: Level: advanced
364: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
365: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Inv()
366: @*/
367: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
368: {
369: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
370: PetscBool same, hasMult;
371: MPI_Comm comm = PetscObjectComm((PetscObject)B);
372: Mat Amat, Pmat;
377: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
380: VecCheckMatCompatible(B, X, 2, Y, 3);
381: if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
382: /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
383: if (lmvm->user_pc) {
384: PCGetOperators(lmvm->J0pc, &Amat, &Pmat);
385: } else if (lmvm->user_ksp) {
386: KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat);
387: } else {
388: Amat = lmvm->J0;
389: }
390: MatHasOperation(Amat, MATOP_MULT, &hasMult);
391: if (hasMult) {
392: /* product is available, use it */
393: MatMult(Amat, X, Y);
394: } else {
395: /* there's no product, so treat J0 as identity */
396: VecCopy(X, Y);
397: }
398: } else if (lmvm->user_scale) {
399: if (lmvm->J0diag) {
400: /* User has defined a diagonal vector for J0 */
401: VecPointwiseMult(X, lmvm->J0diag, Y);
402: } else {
403: /* User has defined a scalar value for J0 */
404: VecCopy(X, Y);
405: VecScale(Y, lmvm->J0scalar);
406: }
407: } else {
408: /* There is no J0 representation so just apply an identity matrix */
409: VecCopy(X, Y);
410: }
411: return 0;
412: }
414: /*------------------------------------------------------------*/
416: /*@
417: MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
418: inverse to the given vector. The specific form of the application
419: depends on whether the user provided a scaling factor, a J0 matrix,
420: a J0 PC, or a J0 KSP object. If no form of the initial Jacobian is
421: provided, the function simply does an identity matrix application
422: (vector copy).
424: Input Parameters:
425: + B - An LMVM-type matrix
426: - X - vector to "multiply" with J0^{-1}
428: Output Parameter:
429: . Y - resulting vector for the operation
431: Level: advanced
433: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
434: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Fwd()
435: @*/
436: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
437: {
438: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
439: PetscBool same, hasSolve;
440: MPI_Comm comm = PetscObjectComm((PetscObject)B);
445: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
448: VecCheckMatCompatible(B, X, 2, Y, 3);
450: /* Invert the initial Jacobian onto q (or apply scaling) */
451: if (lmvm->user_pc) {
452: /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
453: PCApply(lmvm->J0pc, X, Y);
454: } else if (lmvm->user_ksp) {
455: /* User has defined a J0 or a custom KSP so just perform a solution */
456: KSPSolve(lmvm->J0ksp, X, Y);
457: } else if (lmvm->J0) {
458: MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve);
459: if (hasSolve) {
460: MatSolve(lmvm->J0, X, Y);
461: } else {
462: KSPSolve(lmvm->J0ksp, X, Y);
463: }
464: } else if (lmvm->user_scale) {
465: if (lmvm->J0diag) {
466: VecPointwiseDivide(X, Y, lmvm->J0diag);
467: } else {
468: VecCopy(X, Y);
469: VecScale(Y, 1.0/lmvm->J0scalar);
470: }
471: } else {
472: /* There is no J0 representation so just apply an identity matrix */
473: VecCopy(X, Y);
474: }
475: return 0;
476: }
478: /*------------------------------------------------------------*/
480: /*@
481: MatLMVMIsAllocated - Returns a boolean flag that shows whether
482: the necessary data structures for the underlying matrix is allocated.
484: Input Parameters:
485: . B - An LMVM-type matrix
487: Output Parameter:
488: . flg - PETSC_TRUE if allocated, PETSC_FALSE otherwise
490: Level: intermediate
492: .seealso: MatLMVMAllocate(), MatLMVMReset()
493: @*/
495: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
496: {
497: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
498: PetscBool same;
501: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
503: *flg = PETSC_FALSE;
504: if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
505: return 0;
506: }
508: /*------------------------------------------------------------*/
510: /*@
511: MatLMVMAllocate - Produces all necessary common memory for
512: LMVM approximations based on the solution and function vectors
513: provided. If MatSetSizes() and MatSetUp() have not been called
514: before MatLMVMAllocate(), the allocation will read sizes from
515: the provided vectors and update the matrix.
517: Input Parameters:
518: + B - An LMVM-type matrix
519: . X - Solution vector
520: - F - Function vector
522: Level: intermediate
524: .seealso: MatLMVMReset(), MatLMVMUpdate()
525: @*/
526: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
527: {
528: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
529: PetscBool same;
534: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
536: (*lmvm->ops->allocate)(B, X, F);
537: if (lmvm->J0) {
538: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
539: if (same) {
540: MatLMVMAllocate(lmvm->J0, X, F);
541: }
542: }
543: return 0;
544: }
546: /*------------------------------------------------------------*/
548: /*@
549: MatLMVMResetShift - Zero the shift factor.
551: Input Parameters:
552: . B - An LMVM-type matrix
554: Level: intermediate
556: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
557: @*/
558: PetscErrorCode MatLMVMResetShift(Mat B)
559: {
560: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
561: PetscBool same;
564: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
566: lmvm->shift = 0.0;
567: return 0;
568: }
570: /*------------------------------------------------------------*/
572: /*@
573: MatLMVMReset - Flushes all of the accumulated updates out of
574: the LMVM approximation. In practice, this will not actually
575: destroy the data associated with the updates. It simply resets
576: counters, which leads to existing data being overwritten, and
577: MatSolve() being applied as if there are no updates. A boolean
578: flag is available to force destruction of the update vectors.
580: If the user has provided another LMVM matrix as J0, the J0
581: matrix is also reset in this function.
583: Input Parameters:
584: + B - An LMVM-type matrix
585: - destructive - flag for enabling destruction of data structures
587: Level: intermediate
589: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
590: @*/
591: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
592: {
593: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
594: PetscBool same;
597: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
599: (*lmvm->ops->reset)(B, destructive);
600: if (lmvm->J0) {
601: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
602: if (same) {
603: MatLMVMReset(lmvm->J0, destructive);
604: }
605: }
606: return 0;
607: }
609: /*------------------------------------------------------------*/
611: /*@
612: MatLMVMSetHistorySize - Set the number of past iterates to be
613: stored for the construction of the limited-memory QN update.
615: Input Parameters:
616: + B - An LMVM-type matrix
617: - hist_size - number of past iterates (default 5)
619: Options Database:
620: . -mat_lmvm_hist_size <m> - set number of past iterates
622: Level: beginner
624: .seealso: MatLMVMGetUpdateCount()
625: @*/
627: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
628: {
629: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
630: PetscBool same;
631: Vec X, F;
634: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
636: if (hist_size > 0) {
637: lmvm->m = hist_size;
638: if (lmvm->allocated && lmvm->m != lmvm->m_old) {
639: VecDuplicate(lmvm->Xprev, &X);
640: VecDuplicate(lmvm->Fprev, &F);
641: MatLMVMReset(B, PETSC_TRUE);
642: MatLMVMAllocate(B, X, F);
643: VecDestroy(&X);
644: VecDestroy(&F);
645: }
647: return 0;
648: }
650: /*------------------------------------------------------------*/
652: /*@
653: MatLMVMGetUpdateCount - Returns the number of accepted updates.
654: This number may be greater than the total number of update vectors
655: stored in the matrix. The counters are reset when MatLMVMReset()
656: is called.
658: Input Parameters:
659: . B - An LMVM-type matrix
661: Output Parameter:
662: . nupdates - number of accepted updates
664: Level: intermediate
666: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
667: @*/
668: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
669: {
670: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
671: PetscBool same;
674: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
676: *nupdates = lmvm->nupdates;
677: return 0;
678: }
680: /*------------------------------------------------------------*/
682: /*@
683: MatLMVMGetRejectCount - Returns the number of rejected updates.
684: The counters are reset when MatLMVMReset() is called.
686: Input Parameters:
687: . B - An LMVM-type matrix
689: Output Parameter:
690: . nrejects - number of rejected updates
692: Level: intermediate
694: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
695: @*/
696: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
697: {
698: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
699: PetscBool same;
702: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
704: *nrejects = lmvm->nrejects;
705: return 0;
706: }