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: }