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: [](ch_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;

 26:   PetscFunctionBegin;
 30:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
 31:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
 32:   if (!lmvm->allocated) {
 33:     PetscCall(MatLMVMAllocate(B, X, F));
 34:   } else {
 35:     VecCheckMatCompatible(B, X, 2, F, 3);
 36:   }
 37:   if (lmvm->J0) {
 38:     /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
 39:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
 40:     if (same) PetscCall(MatLMVMUpdate(lmvm->J0, X, F));
 41:   }
 42:   PetscCall((*lmvm->ops->update)(B, X, F));
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: /*@
 47:    MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
 48:    an identity matrix (scale = 1.0).

 50:    Input Parameter:
 51: .  B - An LMVM-type matrix

 53:    Level: advanced

 55: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
 56: @*/
 57: PetscErrorCode MatLMVMClearJ0(Mat B)
 58: {
 59:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
 60:   PetscBool same;

 62:   PetscFunctionBegin;
 64:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
 65:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
 66:   lmvm->user_pc    = PETSC_FALSE;
 67:   lmvm->user_ksp   = PETSC_FALSE;
 68:   lmvm->user_scale = PETSC_FALSE;
 69:   lmvm->J0scalar   = 1.0;
 70:   PetscCall(VecDestroy(&lmvm->J0diag));
 71:   PetscCall(MatDestroy(&lmvm->J0));
 72:   PetscCall(PCDestroy(&lmvm->J0pc));
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }

 76: /*@
 77:    MatLMVMSetJ0Scale - Allows the user to define a scalar value
 78:    mu such that J0 = mu*I.

 80:    Input Parameters:
 81: +  B - An LMVM-type matrix
 82: -  scale - Scalar value mu that defines the initial Jacobian

 84:    Level: advanced

 86: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetDiagScale()`, `MatLMVMSetJ0()`
 87: @*/
 88: PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
 89: {
 90:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
 91:   PetscBool same;

 93:   PetscFunctionBegin;
 95:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
 96:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
 97:   PetscCheck(lmvm->square, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Scaling is available only for square LMVM matrices");
 98:   PetscCall(MatLMVMClearJ0(B));
 99:   lmvm->J0scalar   = scale;
100:   lmvm->user_scale = PETSC_TRUE;
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: /*@
105:    MatLMVMSetJ0Diag - Allows the user to define a vector
106:    V such that J0 = diag(V).

108:    Input Parameters:
109: +  B - An LMVM-type matrix
110: -  V - Vector that defines the diagonal of the initial Jacobian

112:    Level: advanced

114: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetScale()`, `MatLMVMSetJ0()`
115: @*/
116: PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
117: {
118:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
119:   PetscBool same;
120:   MPI_Comm  comm = PetscObjectComm((PetscObject)B);

122:   PetscFunctionBegin;
125:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
126:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
127:   PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "Matrix must be allocated before setting diagonal scaling");
128:   PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Diagonal scaling is available only for square LMVM matrices");
129:   VecCheckSameSize(V, 2, lmvm->Fprev, 3);

131:   PetscCall(MatLMVMClearJ0(B));
132:   if (!lmvm->J0diag) PetscCall(VecDuplicate(V, &lmvm->J0diag));
133:   PetscCall(VecCopy(V, lmvm->J0diag));
134:   lmvm->user_scale = PETSC_TRUE;
135:   PetscFunctionReturn(PETSC_SUCCESS);
136: }

138: /*@
139:    MatLMVMSetJ0 - Allows the user to define the initial
140:    Jacobian matrix from which the LMVM-type approximation is
141:    built up. Inverse of this initial Jacobian is applied
142:    using an internal `KSP` solver, which defaults to `KSPGMRES`.
143:    This internal `KSP` solver has the "mat_lmvm_" option
144:    prefix.

146:    Note that another LMVM-type matrix can be used in place of
147:    J0, in which case updating the outer LMVM-type matrix will
148:    also trigger the update for the inner LMVM-type matrix. This
149:    is useful in cases where a full-memory diagonal approximation
150:    such as `MATLMVMDIAGBRDN` is used in place of J0.

152:    Input Parameters:
153: +  B - An LMVM-type matrix
154: -  J0 - The initial Jacobian matrix

156:    Level: advanced

158: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`
159: @*/
160: PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
161: {
162:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
163:   PetscBool same;

165:   PetscFunctionBegin;
168:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
169:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
170:   PetscCall(MatLMVMClearJ0(B));
171:   PetscCall(MatDestroy(&lmvm->J0));
172:   PetscCall(PetscObjectReference((PetscObject)J0));
173:   lmvm->J0 = J0;
174:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
175:   if (!same && lmvm->square) PetscCall(KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0));
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: /*@
180:    MatLMVMSetJ0PC - Allows the user to define a `PC` object that
181:    acts as the initial inverse-Jacobian matrix. This `PC` should
182:    already contain all the operators necessary for its application.
183:    The LMVM-type matrix only calls `PCApply()` without changing any other
184:    options.

186:    Input Parameters:
187: +  B - An LMVM-type matrix
188: -  J0pc - `PC` object where `PCApply()` defines an inverse application for J0

190:    Level: advanced

192: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0PC()`
193: @*/
194: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
195: {
196:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
197:   PetscBool same;
198:   MPI_Comm  comm = PetscObjectComm((PetscObject)B);

200:   PetscFunctionBegin;
203:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
204:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
205:   PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
206:   PetscCall(MatLMVMClearJ0(B));
207:   PetscCall(PetscObjectReference((PetscObject)J0pc));
208:   lmvm->J0pc    = J0pc;
209:   lmvm->user_pc = PETSC_TRUE;
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: /*@
214:    MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
215:    KSP solver for the initial inverse-Jacobian approximation.
216:    This `KSP` solver should already contain all the operators
217:    necessary to perform the inversion. The LMVM-type matrix only
218:    calls `KSPSolve()` without changing any other options.

220:    Input Parameters:
221: +  B - An LMVM-type matrix
222: -  J0ksp - `KSP` solver for the initial inverse-Jacobian application

224:    Level: advanced

226: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetJ0KSP()`
227: @*/
228: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
229: {
230:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
231:   PetscBool same;
232:   MPI_Comm  comm = PetscObjectComm((PetscObject)B);

234:   PetscFunctionBegin;
237:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
238:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
239:   PetscCheck(lmvm->square, comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
240:   PetscCall(MatLMVMClearJ0(B));
241:   PetscCall(KSPDestroy(&lmvm->J0ksp));
242:   PetscCall(PetscObjectReference((PetscObject)J0ksp));
243:   lmvm->J0ksp    = J0ksp;
244:   lmvm->user_ksp = PETSC_TRUE;
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: /*@
249:    MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.

251:    Input Parameters:
252: .  B - An LMVM-type matrix

254:    Output Parameter:
255: .  J0 - `Mat` object for defining the initial Jacobian

257:    Level: advanced

259: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`
260: @*/
261: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
262: {
263:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
264:   PetscBool same;

266:   PetscFunctionBegin;
268:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
269:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
270:   *J0 = lmvm->J0;
271:   PetscFunctionReturn(PETSC_SUCCESS);
272: }

274: /*@
275:    MatLMVMGetJ0PC - Returns a pointer to the internal `PC` object
276:    associated with the initial Jacobian.

278:    Input Parameter:
279: .  B - An LMVM-type matrix

281:    Output Parameter:
282: .  J0pc - `PC` object for defining the initial inverse-Jacobian

284:    Level: advanced

286: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0PC()`
287: @*/
288: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
289: {
290:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
291:   PetscBool same;

293:   PetscFunctionBegin;
295:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
296:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
297:   if (lmvm->J0pc) {
298:     *J0pc = lmvm->J0pc;
299:   } else {
300:     PetscCall(KSPGetPC(lmvm->J0ksp, J0pc));
301:   }
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: /*@
306:    MatLMVMGetJ0KSP - Returns a pointer to the internal `KSP` solver
307:    associated with the initial Jacobian.

309:    Input Parameter:
310: .  B - An LMVM-type matrix

312:    Output Parameter:
313: .  J0ksp - `KSP` solver for defining the initial inverse-Jacobian

315:    Level: advanced

317: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0KSP()`
318: @*/
319: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
320: {
321:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
322:   PetscBool same;

324:   PetscFunctionBegin;
326:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
327:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
328:   *J0ksp = lmvm->J0ksp;
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: /*@
333:    MatLMVMApplyJ0Fwd - Applies an approximation of the forward
334:    matrix-vector product with the initial Jacobian.

336:    Input Parameters:
337: +  B - An LMVM-type matrix
338: -  X - vector to multiply with J0

340:    Output Parameter:
341: .  Y - resulting vector for the operation

343:    Level: advanced

345: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
346:           `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Inv()`
347: @*/
348: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
349: {
350:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
351:   PetscBool same, hasMult;
352:   MPI_Comm  comm = PetscObjectComm((PetscObject)B);
353:   Mat       Amat, Pmat;

355:   PetscFunctionBegin;
359:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
360:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
361:   PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
362:   VecCheckMatCompatible(B, X, 2, Y, 3);
363:   if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
364:     /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
365:     if (lmvm->user_pc) {
366:       PetscCall(PCGetOperators(lmvm->J0pc, &Amat, &Pmat));
367:     } else if (lmvm->user_ksp) {
368:       PetscCall(KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat));
369:     } else {
370:       Amat = lmvm->J0;
371:     }
372:     PetscCall(MatHasOperation(Amat, MATOP_MULT, &hasMult));
373:     if (hasMult) {
374:       /* product is available, use it */
375:       PetscCall(MatMult(Amat, X, Y));
376:     } else {
377:       /* there's no product, so treat J0 as identity */
378:       PetscCall(VecCopy(X, Y));
379:     }
380:   } else if (lmvm->user_scale) {
381:     if (lmvm->J0diag) {
382:       /* User has defined a diagonal vector for J0 */
383:       PetscCall(VecPointwiseMult(X, lmvm->J0diag, Y));
384:     } else {
385:       /* User has defined a scalar value for J0 */
386:       PetscCall(VecAXPBY(Y, lmvm->J0scalar, 0.0, X));
387:     }
388:   } else {
389:     /* There is no J0 representation so just apply an identity matrix */
390:     PetscCall(VecCopy(X, Y));
391:   }
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: /*@
396:    MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
397:    inverse to the given vector. The specific form of the application
398:    depends on whether the user provided a scaling factor, a J0 matrix,
399:    a J0 `PC`, or a J0 `KSP` object. If no form of the initial Jacobian is
400:    provided, the function simply does an identity matrix application
401:    (vector copy).

403:    Input Parameters:
404: +  B - An LMVM-type matrix
405: -  X - vector to "multiply" with J0^{-1}

407:    Output Parameter:
408: .  Y - resulting vector for the operation

410:    Level: advanced

412: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMSetJ0()`, `MatLMVMSetJ0Scale()`, `MatLMVMSetJ0ScaleDiag()`,
413:           `MatLMVMSetJ0PC()`, `MatLMVMSetJ0KSP()`, `MatLMVMApplyJ0Fwd()`
414: @*/
415: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
416: {
417:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
418:   PetscBool same, hasSolve;
419:   MPI_Comm  comm = PetscObjectComm((PetscObject)B);

421:   PetscFunctionBegin;
425:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
426:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
427:   PetscCheck(lmvm->allocated, comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
428:   VecCheckMatCompatible(B, X, 2, Y, 3);

430:   /* Invert the initial Jacobian onto q (or apply scaling) */
431:   if (lmvm->user_pc) {
432:     /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
433:     PetscCall(PCApply(lmvm->J0pc, X, Y));
434:   } else if (lmvm->user_ksp) {
435:     /* User has defined a J0 or a custom KSP so just perform a solution */
436:     PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
437:   } else if (lmvm->J0) {
438:     PetscCall(MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve));
439:     if (hasSolve) {
440:       PetscCall(MatSolve(lmvm->J0, X, Y));
441:     } else {
442:       PetscCall(KSPSolve(lmvm->J0ksp, X, Y));
443:     }
444:   } else if (lmvm->user_scale) {
445:     if (lmvm->J0diag) {
446:       PetscCall(VecPointwiseDivide(X, Y, lmvm->J0diag));
447:     } else {
448:       PetscCall(VecAXPBY(Y, 1.0 / lmvm->J0scalar, 0.0, X));
449:     }
450:   } else {
451:     /* There is no J0 representation so just apply an identity matrix */
452:     PetscCall(VecCopy(X, Y));
453:   }
454:   PetscFunctionReturn(PETSC_SUCCESS);
455: }

457: /*@
458:    MatLMVMIsAllocated - Returns a boolean flag that shows whether
459:    the necessary data structures for the underlying matrix is allocated.

461:    Input Parameter:
462: .  B - An LMVM-type matrix

464:    Output Parameter:
465: .  flg - `PETSC_TRUE` if allocated, `PETSC_FALSE` otherwise

467:    Level: intermediate

469: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMReset()`
470: @*/

472: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
473: {
474:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
475:   PetscBool same;

477:   PetscFunctionBegin;
479:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
480:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
481:   *flg = PETSC_FALSE;
482:   if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: /*@
487:    MatLMVMAllocate - Produces all necessary common memory for
488:    LMVM approximations based on the solution and function vectors
489:    provided. If `MatSetSizes()` and `MatSetUp()` have not been called
490:    before `MatLMVMAllocate()`, the allocation will read sizes from
491:    the provided vectors and update the matrix.

493:    Input Parameters:
494: +  B - An LMVM-type matrix
495: .  X - Solution vector
496: -  F - Function vector

498:    Level: intermediate

500: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMReset()`, `MatLMVMUpdate()`
501: @*/
502: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
503: {
504:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
505:   PetscBool same;
506:   VecType   vtype;

508:   PetscFunctionBegin;
512:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
513:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
514:   PetscCall(VecGetType(X, &vtype));
515:   PetscCall(MatSetVecType(B, vtype));
516:   PetscCall((*lmvm->ops->allocate)(B, X, F));
517:   if (lmvm->J0) {
518:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
519:     if (same) PetscCall(MatLMVMAllocate(lmvm->J0, X, F));
520:   }
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: /*@
525:    MatLMVMResetShift - Zero the shift factor.

527:    Input Parameter:
528: .  B - An LMVM-type matrix

530:    Level: intermediate

532: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
533: @*/
534: PetscErrorCode MatLMVMResetShift(Mat B)
535: {
536:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
537:   PetscBool same;

539:   PetscFunctionBegin;
541:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
542:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
543:   lmvm->shift = 0.0;
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: /*@
548:    MatLMVMReset - Flushes all of the accumulated updates out of
549:    the LMVM approximation. In practice, this will not actually
550:    destroy the data associated with the updates. It simply resets
551:    counters, which leads to existing data being overwritten, and
552:    `MatSolve()` being applied as if there are no updates. A boolean
553:    flag is available to force destruction of the update vectors.

555:    If the user has provided another LMVM matrix as J0, the J0
556:    matrix is also reset in this function.

558:    Input Parameters:
559: +  B - An LMVM-type matrix
560: -  destructive - flag for enabling destruction of data structures

562:    Level: intermediate

564: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMAllocate()`, `MatLMVMUpdate()`
565: @*/
566: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
567: {
568:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
569:   PetscBool same;

571:   PetscFunctionBegin;
573:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
574:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
575:   PetscCall((*lmvm->ops->reset)(B, destructive));
576:   if (lmvm->J0) {
577:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same));
578:     if (same) PetscCall(MatLMVMReset(lmvm->J0, destructive));
579:   }
580:   PetscFunctionReturn(PETSC_SUCCESS);
581: }

583: /*@
584:    MatLMVMSetHistorySize - Set the number of past iterates to be
585:    stored for the construction of the limited-memory QN update.

587:    Input Parameters:
588: +  B - An LMVM-type matrix
589: -  hist_size - number of past iterates (default 5)

591:    Options Database Key:
592: .  -mat_lmvm_hist_size <m> - set number of past iterates

594:    Level: beginner

596: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetUpdateCount()`
597: @*/
598: PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
599: {
600:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
601:   PetscBool same;
602:   Vec       X, F;

604:   PetscFunctionBegin;
606:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
607:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
608:   if (hist_size > 0) {
609:     lmvm->m = hist_size;
610:     if (lmvm->allocated && lmvm->m != lmvm->m_old) {
611:       PetscCall(VecDuplicate(lmvm->Xprev, &X));
612:       PetscCall(VecDuplicate(lmvm->Fprev, &F));
613:       PetscCall(MatLMVMReset(B, PETSC_TRUE));
614:       PetscCall(MatLMVMAllocate(B, X, F));
615:       PetscCall(VecDestroy(&X));
616:       PetscCall(VecDestroy(&F));
617:     }
618:   } else PetscCheck(hist_size >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a non-negative integer.");
619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }

622: /*@
623:    MatLMVMGetUpdateCount - Returns the number of accepted updates.
624:    This number may be greater than the total number of update vectors
625:    stored in the matrix. The counters are reset when `MatLMVMReset()`
626:    is called.

628:    Input Parameter:
629: .  B - An LMVM-type matrix

631:    Output Parameter:
632: .  nupdates - number of accepted updates

634:    Level: intermediate

636: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
637: @*/
638: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
639: {
640:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
641:   PetscBool same;

643:   PetscFunctionBegin;
645:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
646:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
647:   *nupdates = lmvm->nupdates;
648:   PetscFunctionReturn(PETSC_SUCCESS);
649: }

651: /*@
652:    MatLMVMGetRejectCount - Returns the number of rejected updates.
653:    The counters are reset when `MatLMVMReset()` is called.

655:    Input Parameter:
656: .  B - An LMVM-type matrix

658:    Output Parameter:
659: .  nrejects - number of rejected updates

661:    Level: intermediate

663: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MATLMVM`, `MatLMVMGetRejectCount()`, `MatLMVMReset()`
664: @*/
665: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
666: {
667:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
668:   PetscBool same;

670:   PetscFunctionBegin;
672:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
673:   PetscCheck(same, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
674:   *nrejects = lmvm->nrejects;
675:   PetscFunctionReturn(PETSC_SUCCESS);
676: }