Actual source code: lmvmutils.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  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.
  8:    
  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.
162:    
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:   if (!J0->assembled) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "J0 is not assembled.");
190:   if (B->symmetric && (!J0->symmetric)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "J0 and J0pre must be symmetric when B is symmetric");
191:   if (lmvm->allocated) {
192:     MatCheckSameSize(B, 1, J0, 2);
193:   }
194:   MatLMVMClearJ0(B);
195:   MatDestroy(&lmvm->J0);
196:   PetscObjectReference((PetscObject)J0);
197:   lmvm->J0 = J0;
198:   PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
199:   if (!same && lmvm->square) {
200:     KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0);
201:   }
202:   return(0);
203: }

205: /*------------------------------------------------------------*/

207: /*@
208:    MatLMVMSetJ0PC - Allows the user to define a PC object that 
209:    acts as the initial inverse-Jacobian matrix. This PC should 
210:    already contain all the operators necessary for its application. 
211:    The LMVM matrix only calls PCApply() without changing any other 
212:    options.

214:    Input Parameters:
215: +  B - An LMVM-type matrix
216: -  J0pc - PC object where PCApply defines an inverse application for J0

218:    Level: advanced

220: .seealso: MatLMVMGetJ0PC()
221: @*/
222: PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
223: {
224:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
225:   PetscErrorCode    ierr;
226:   PetscBool         same;
227:   MPI_Comm          comm = PetscObjectComm((PetscObject)B);
228:   Mat               J0, J0pre;

233:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
234:   if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
235:   if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
236:   PCGetOperators(J0pc, &J0, &J0pre);
237:   if (B->symmetric && (!J0->symmetric || !J0pre->symmetric)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "J0 and J0pre must be symmetric when B is symmetric");
238:   if (lmvm->allocated) {
239:     MatCheckSameSize(B, 1, J0, 2);
240:     MatCheckSameSize(B, 1, J0pre, 3);
241:   }
242:   MatDestroy(&J0);
243:   MatDestroy(&J0pre);
244:   MatLMVMClearJ0(B);
245:   PetscObjectReference((PetscObject)J0pc);
246:   lmvm->J0pc = J0pc;
247:   lmvm->user_pc = PETSC_TRUE;
248:   return(0);
249: }

251: /*------------------------------------------------------------*/

253: /*@
254:    MatLMVMSetJ0KSP - Allows the user to provide a pre-configured 
255:    KSP solver for the initial inverse-Jacobian approximation. 
256:    This KSP solver should already contain all the operators 
257:    necessary to perform the inversion. The LMVM matrix only 
258:    calls KSPSolve() without changing any other options.

260:    Input Parameters:
261: +  B - An LMVM-type matrix
262: -  J0ksp - KSP solver for the initial inverse-Jacobian application

264:    Level: advanced

266: .seealso: MatLMVMGetJ0KSP()
267: @*/
268: PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
269: {
270:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
271:   PetscErrorCode    ierr;
272:   PetscBool         same;
273:   MPI_Comm          comm = PetscObjectComm((PetscObject)B);
274:   Mat               J0, J0pre;

279:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
280:   if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
281:   if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
282:   KSPGetOperators(J0ksp, &J0, &J0pre);
283:   if (B->symmetric && (!J0->symmetric || !J0pre->symmetric)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "J0 and J0pre must be symmetric when B is symmetric");
284:   if (lmvm->allocated) {
285:     MatCheckSameSize(B, 1, J0, 2);
286:     MatCheckSameSize(B, 1, J0pre, 3);
287:   }
288:   MatLMVMClearJ0(B);
289:   KSPDestroy(&lmvm->J0ksp);
290:   PetscObjectReference((PetscObject)J0ksp);
291:   lmvm->J0ksp = J0ksp;
292:   lmvm->user_ksp = PETSC_TRUE;
293:   return(0);
294: }

296: /*------------------------------------------------------------*/

298: /*@
299:    MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.

301:    Input Parameters:
302: .  B - An LMVM-type matrix

304:    Output Parameter:
305: .  J0 - Mat object for defining the initial Jacobian

307:    Level: advanced

309: .seealso: MatLMVMSetJ0()
310: @*/
311: PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
312: {
313:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
314:   PetscErrorCode    ierr;
315:   PetscBool         same;

319:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
320:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
321:   *J0 = lmvm->J0;
322:   return(0);
323: }

325: /*------------------------------------------------------------*/

327: /*@
328:    MatLMVMGetJ0PC - Returns a pointer to the internal PC object 
329:    associated with the initial Jacobian.

331:    Input Parameters:
332: .  B - An LMVM-type matrix

334:    Output Parameter:
335: .  J0pc - PC object for defining the initial inverse-Jacobian

337:    Level: advanced

339: .seealso: MatLMVMSetJ0PC()
340: @*/
341: PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
342: {
343:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
344:   PetscErrorCode    ierr;
345:   PetscBool         same;

349:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
350:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
351:   if (lmvm->J0pc) {
352:     *J0pc = lmvm->J0pc;
353:   } else {
354:     KSPGetPC(lmvm->J0ksp, J0pc);
355:   }
356:   return(0);
357: }

359: /*------------------------------------------------------------*/

361: /*@
362:    MatLMVMGetJ0KSP - Returns a pointer to the internal KSP solver 
363:    associated with the initial Jacobian.

365:    Input Parameters:
366: .  B - An LMVM-type matrix

368:    Output Parameter:
369: .  J0ksp - KSP solver for defining the initial inverse-Jacobian

371:    Level: advanced

373: .seealso: MatLMVMSetJ0KSP()
374: @*/
375: PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
376: {
377:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
378:   PetscErrorCode    ierr;
379:   PetscBool         same;

383:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
384:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
385:   *J0ksp = lmvm->J0ksp;
386:   return(0);
387: }

389: /*------------------------------------------------------------*/

391: /*@
392:    MatLMVMApplyJ0Fwd - Applies an approximation of the forward 
393:    matrix-vector product with the initial Jacobian.

395:    Input Parameters:
396: +  B - An LMVM-type matrix
397: -  X - vector to multiply with J0

399:    Output Parameter:
400: .  Y - resulting vector for the operation

402:    Level: advanced

404: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(), 
405:           MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Inv()
406: @*/
407: PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
408: {
409:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
410:   PetscErrorCode    ierr;
411:   PetscBool         same, hasMult;
412:   MPI_Comm          comm = PetscObjectComm((PetscObject)B);
413:   Mat               Amat, Pmat;

419:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
420:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
421:   if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
422:   VecCheckMatCompatible(B, X, 2, Y, 3);
423:   if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
424:     /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
425:     if (lmvm->user_pc){
426:       PCGetOperators(lmvm->J0pc, &Amat, &Pmat);
427:     } else if (lmvm->user_ksp) {
428:       KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat);
429:     } else {
430:       Amat = lmvm->J0;
431:     }
432:     MatHasOperation(Amat, MATOP_MULT, &hasMult);
433:     if (hasMult) {
434:       /* product is available, use it */
435:       MatMult(Amat, X, Y);
436:     } else {
437:       /* there's no product, so treat J0 as identity */
438:       VecCopy(X, Y);
439:     }
440:   } else if (lmvm->user_scale) {
441:     if (lmvm->J0diag) {
442:       /* User has defined a diagonal vector for J0 */
443:       VecPointwiseMult(X, lmvm->J0diag, Y);
444:     } else {
445:       /* User has defined a scalar value for J0 */
446:       VecCopy(X, Y);
447:       VecScale(Y, lmvm->J0scalar);
448:     }
449:   } else {
450:     /* There is no J0 representation so just apply an identity matrix */
451:     VecCopy(X, Y);
452:   }
453:   return(0);
454: }

456: /*------------------------------------------------------------*/

458: /*@
459:    MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian 
460:    inverse to the given vector. The specific form of the application 
461:    depends on whether the user provided a scaling factor, a J0 matrix, 
462:    a J0 PC, or a J0 KSP object. If no form of the initial Jacobian is 
463:    provided, the function simply does an identity matrix application 
464:    (vector copy).

466:    Input Parameters:
467: +  B - An LMVM-type matrix
468: -  X - vector to "multiply" with J0^{-1}

470:    Output Parameter:
471: .  Y - resulting vector for the operation

473:    Level: advanced

475: .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(), 
476:           MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Fwd()
477: @*/
478: PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
479: {
480:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
481:   PetscErrorCode    ierr;
482:   PetscBool         same, hasSolve;
483:   MPI_Comm          comm = PetscObjectComm((PetscObject)B);

489:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
490:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
491:   if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
492:   VecCheckMatCompatible(B, X, 2, Y, 3);
493:   /* Invert the initial Jacobian onto q (or apply scaling) */
494:   if (lmvm->user_pc) {
495:     /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
496:     PCApply(lmvm->J0pc, X, Y);
497:   } else if (lmvm->user_ksp) {
498:     /* User has defined a J0 or a custom KSP so just perform a solution */
499:     KSPSolve(lmvm->J0ksp, X, Y);
500:   } else if (lmvm->J0) {
501:     MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve);
502:     if (hasSolve) {
503:       MatSolve(lmvm->J0, X, Y);
504:     } else {
505:       KSPSolve(lmvm->J0ksp, X, Y);
506:     }
507:   } else if (lmvm->user_scale) {
508:     if (lmvm->J0diag) {
509:       VecPointwiseDivide(X, Y, lmvm->J0diag);
510:     } else {
511:       VecCopy(X, Y);
512:       VecScale(Y, 1.0/lmvm->J0scalar);
513:     }
514:   } else {
515:     /* There is no J0 representation so just apply an identity matrix */
516:     VecCopy(X, Y);
517:   }
518:   return(0);
519: }

521: /*------------------------------------------------------------*/

523: /*@
524:    MatLMVMIsAllocated - Returns a boolean flag that shows whether 
525:    the necessary data structures for the underlying matrix is allocated.

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

530:    Output Parameter:
531: .  flg - PETSC_TRUE if allocated, PETSC_FALSE otherwise

533:    Level: intermediate

535: .seealso: MatLMVMAllocate(), MatLMVMReset()
536: @*/

538: PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
539: {
540:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
541:   PetscErrorCode    ierr;
542:   PetscBool         same;
543: 
546:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
547:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
548:   *flg = PETSC_FALSE;
549:   if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
550:   return(0);
551: }

553: /*------------------------------------------------------------*/

555: /*@
556:    MatLMVMAllocate - Produces all necessary common memory for 
557:    LMVM approximations based on the solution and function vectors
558:    provided. If MatSetSizes() and MatSetUp() have not been called 
559:    before MatLMVMAllocate(), the allocation will read sizes from 
560:    the provided vectors and update the matrix.

562:    Input Parameters:
563: +  B - An LMVM-type matrix
564: .  X - Solution vector
565: -  F - Function vector

567:    Level: intermediate

569: .seealso: MatLMVMReset(), MatLMVMUpdate()
570: @*/
571: PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
572: {
573:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
574:   PetscErrorCode    ierr;
575:   PetscBool         same;

581:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
582:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
583:   (*lmvm->ops->allocate)(B, X, F);
584:   if (lmvm->J0) {
585:     PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
586:     if (same) {
587:       MatLMVMAllocate(lmvm->J0, X, F);
588:     }
589:   }
590:   return(0);
591: }

593: /*------------------------------------------------------------*/

595: /*@
596:    MatLMVMResetShift - Zero the shift factor.

598:    Input Parameters:
599: .  B - An LMVM-type matrix

601:    Level: intermediate

603: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
604: @*/
605: PetscErrorCode MatLMVMResetShift(Mat B)
606: {
607:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
608:   PetscErrorCode    ierr;
609:   PetscBool         same;

613:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
614:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
615:   lmvm->shift = 0.0;
616:   return(0);
617: }

619: /*------------------------------------------------------------*/

621: /*@
622:    MatLMVMReset - Flushes all of the accumulated updates out of 
623:    the LMVM approximation. In practice, this will not actually 
624:    destroy the data associated with the updates. It simply resets 
625:    counters, which leads to existing data being overwritten, and 
626:    MatSolve() being applied as if there are no updates. A boolean 
627:    flag is available to force destruction of the update vectors.
628:    
629:    If the user has provided another LMVM matrix as J0, the J0 
630:    matrix is also reset in this function.

632:    Input Parameters:
633: +  B - An LMVM-type matrix
634: -  destructive - flag for enabling destruction of data structures

636:    Level: intermediate

638: .seealso: MatLMVMAllocate(), MatLMVMUpdate()
639: @*/
640: PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
641: {
642:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
643:   PetscErrorCode    ierr;
644:   PetscBool         same;

648:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
649:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
650:   (*lmvm->ops->reset)(B, destructive);
651:   if (lmvm->J0) {
652:     PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);
653:     if (same) {
654:       MatLMVMReset(lmvm->J0, destructive);
655:     }
656:   }
657:   return(0);
658: }

660: /*------------------------------------------------------------*/

662: /*@
663:    MatLMVMGetUpdateCount - Returns the number of accepted updates.
664:    This number may be greater than the total number of update vectors 
665:    stored in the matrix. The counters are reset when MatLMVMReset() 
666:    is called.

668:    Input Parameters:
669: .  B - An LMVM-type matrix

671:    Output Parameter:
672: .  nupdates - number of accepted updates

674:    Level: intermediate

676: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
677: @*/
678: PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
679: {
680:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
681:   PetscErrorCode    ierr;
682:   PetscBool         same;
683: 
686:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
687:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
688:   *nupdates = lmvm->nupdates;
689:   return(0);
690: }

692: /*------------------------------------------------------------*/

694: /*@
695:    MatLMVMGetRejectCount - Returns the number of rejected updates. 
696:    The counters are reset when MatLMVMReset() is called.

698:    Input Parameters:
699: .  B - An LMVM-type matrix

701:    Output Parameter:
702: .  nrejects - number of rejected updates

704:    Level: intermediate

706: .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
707: @*/
708: PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
709: {
710:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
711:   PetscErrorCode    ierr;
712:   PetscBool         same;
713: 
716:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
717:   if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
718:   *nrejects = lmvm->nrejects;
719:   return(0);
720: }