Actual source code: lmvmutils.c

petsc-3.13.3 2020-07-01
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:   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;
522:   
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.
607:    
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;
663:   
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:   }
681:   
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;
708:   
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;
738:   
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: }