Actual source code: hpddm.cxx

petsc-master 2019-11-16
Report Typos and Errors
  1:  #include <petsc/private/petschpddm.h>

  3: static PetscBool citeKSP = PETSC_FALSE;
  4: static const char hpddmCitationKSP[] = "@inproceedings{jolivet2016block,\n\tTitle = {{Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers}},\n\tAuthor = {Jolivet, Pierre and Tournier, Pierre-Henri},\n\tOrganization = {IEEE},\n\tYear = {2016},\n\tSeries = {SC16},\n\tBooktitle = {Proceedings of the 2016 International Conference for High Performance Computing, Networking, Storage and Analysis}\n}\n";

  6: static PetscErrorCode KSPSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, KSP ksp)
  7: {
  8:   PetscReal      r;
  9:   PetscInt       i;

 13:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 14:   PetscOptionsReal("-ksp_richardson_scale", "Damping factor used in Richardson iterations", "none", 1.0, &r, NULL);
 15:   PetscOptionsInt("-ksp_gmres_restart", "Maximum number of Arnoldi vectors generated per cycle", "none", 40, &i, NULL);
 16:   PetscOptionsEList("-ksp_hpddm_krylov_method", "Type of Krylov method", "none", HPDDMKrylovMethod, 7, HPDDMKrylovMethod[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 17:   PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "none", -1.0, &r, NULL);
 18:   PetscOptionsInt("-ksp_hpddm_enlarge_krylov_subspace", "Split the initial right-hand side into multiple vectors", "none", 1, &i, NULL);
 19:   PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "none", HPDDMOrthogonalization, 2, HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
 20:   PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "none", HPDDMQR, 3, HPDDMQR[HPDDM_QR_CHOLQR], &i, NULL);
 21:   i = HPDDM_VARIANT_LEFT;
 22:   PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "none", HPDDMVariant, 3, HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 23:   if (i > 0) {
 24:     KSPSetPCSide(ksp, PC_RIGHT);
 25:   }
 26:   PetscOptionsInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "none", 0, &i, NULL);
 27:   PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "none", HPDDMRecycleTarget, 6, HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
 28:   PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "none", HPDDMRecycleStrategy, 2, HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
 29:   PetscOptionsTail();
 30:   return(0);
 31: }

 33: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
 34: {
 35:   Mat            A;
 36:   PetscInt       n;

 40:   KSPGetOperators(ksp, &A, NULL);
 41:   MatGetLocalSize(A, &n, NULL);
 42:   HPDDM::PETScOperator *op = new HPDDM::PETScOperator(ksp, n, 1);
 43:   ksp->data = reinterpret_cast<void*>(op);
 44:   return(0);
 45: }

 47: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
 48: {
 50:   if (ksp->data) {
 51:     delete reinterpret_cast<HPDDM::PETScOperator*>(ksp->data);
 52:     ksp->data = NULL;
 53:   }
 54:   return(0);
 55: }

 57: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
 58: {

 62:   KSPReset_HPDDM(ksp);
 63:   KSPDestroyDefault(ksp);
 64:   return(0);
 65: }

 67: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
 68: {
 69:   PetscScalar       *x;
 70:   const PetscScalar *b;
 71:   MPI_Comm          comm;
 72:   PetscErrorCode    ierr;

 75:   PetscCitationsRegister(hpddmCitationKSP, &citeKSP);
 76:   VecGetArray(ksp->vec_sol, &x);
 77:   VecGetArrayRead(ksp->vec_rhs, &b);
 78:   PetscObjectGetComm((PetscObject)ksp, &comm);
 79:   const HPDDM::PETScOperator& op = *reinterpret_cast<HPDDM::PETScOperator*>(ksp->data);
 80:   ksp->its = HPDDM::IterativeMethod::solve(op, b, x, 1, comm);
 81:   VecRestoreArrayRead(ksp->vec_rhs, &b);
 82:   VecRestoreArray(ksp->vec_sol, &x);
 83:   if (ksp->its < ksp->max_it) ksp->reason = KSP_CONVERGED_RTOL;
 84:   else ksp->reason = KSP_DIVERGED_ITS;
 85:   return(0);
 86: }

 88: /*MC
 89:      KSPHPDDM - Interface with the HPDDM library.

 91:    This KSP may be used to further select methods that are currently not implemented natively in PETSc, e.g., GCRODR [2006], a recycled Krylov method which is similar to KSPLGMRES, see [2016] for a comparison. ex75.c shows how to reproduce the results from the aforementioned paper [2006]. Here is a chronological bibliography of relevant publications linked with KSP available in HPDDM through KSPHPDDM, and not available directly in PETSc, like KSPCG and KSPGMRES.

 93: .vb
 94:    [1980] The Block Conjugate Gradient Algorithm and Related Methods. O'Leary. Linear Algebra and its Applications.
 95:    [2006] Recycling Krylov Subspaces for Sequences of Linear Systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
 96:    [2013] A Modified Block Flexible GMRES Method with Deflation at Each Iteration for the Solution of Non-Hermitian Linear Systems with Multiple Right-Hand Sides. Calandra, Gratton, Lago, Vasseur, and Carvalho. SIAM Journal on Scientific Computing.
 97:    [2016] Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers. Jolivet and Tournier. SC16.
 98:    [2017] A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.
 99: .ve

101:    Options Database Keys:
102: +   -ksp_richardson_scale <scale, default=1.0> - see KSPRICHARDSON
103: .   -ksp_gmres_restart <restart, default=40> - see KSPGMRES
104: .   -ksp_hpddm_krylov_method <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, or bfbcg
105: .   -ksp_hpddm_deflation_tol <eps, default=\-1.0> - tolerance when deflating right-hand sides inside block methods (no deflation by default, only relevant with block methods)
106: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
107: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
108: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
109: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible
110: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
111: .   -ksp_hpddm_recycle_target <type, default=SM> - criterion to select harmonic Ritz vectors using either SM, LM, SR, LR, SI, or LI (only relevant with GCRODR or BGCRODR)
112: -   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)

114:    Level: intermediate

116: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
117: M*/
118: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
119: {

123:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
124:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
125:   ksp->ops->setup          = KSPSetUp_HPDDM;
126:   ksp->ops->solve          = KSPSolve_HPDDM;
127:   ksp->ops->reset          = KSPReset_HPDDM;
128:   ksp->ops->destroy        = KSPDestroy_HPDDM;
129:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
130:   return(0);
131: }