Actual source code: hpddm.cxx

petsc-3.13.0 2020-03-29
Report Typos and Errors
  1:  #include <petsc/private/petschpddm.h>

  3: /* static array length */
  4: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

  6: static const char *HPDDMType[]              = { "gmres", "bgmres", "cg", "bcg", "gcrodr", "bgcrodr", "bfbcg" };
  7: static const char *HPDDMOrthogonalization[] = { "cgs", "mgs" };
  8: static const char *HPDDMQR[]                = { "cholqr", "cgs", "mgs" };
  9: static const char *HPDDMVariant[]           = { "left", "right", "flexible" };
 10: static const char *HPDDMRecycleTarget[]     = { "SM", "LM", "SR", "LR", "SI", "LI" };
 11: static const char *HPDDMRecycleStrategy[]   = { "A", "B" };

 13: static PetscBool citeKSP = PETSC_FALSE;
 14: 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";

 16: static PetscErrorCode KSPSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, KSP ksp)
 17: {
 18:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
 19:   PetscInt       i, j;

 23:   data->scntl[0] = ksp->max_it;
 24:   data->rcntl[0] = ksp->rtol;
 25:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 26:   i = HPDDM_KRYLOV_METHOD_GMRES;
 27:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDM", HPDDMType, ALEN(HPDDMType), HPDDMType[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 28:   data->cntl[5] = i;
 29:   if (data->cntl[5] == HPDDM_KRYLOV_METHOD_RICHARDSON) {
 30:     data->rcntl[0] = 1.0;
 31:     PetscOptionsReal("-ksp_richardson_scale", "Damping factor used in Richardson iterations", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
 32:   } else {
 33:     i = HPDDM_VARIANT_LEFT;
 34:     PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, ALEN(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 35:     if (i != HPDDM_VARIANT_LEFT && (data->cntl[5] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Right and flexible preconditioned (BF)BCG not implemented");
 36:     data->cntl[1] = i;
 37:     if (i > 0) {
 38:       KSPSetPCSide(ksp, PC_RIGHT);
 39:     }
 40:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG) {
 41:       data->rcntl[1] = -1.0;
 42:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[1], data->rcntl + 1, NULL);
 43:       i = 1;
 44:       PetscOptionsRangeInt("-ksp_hpddm_enlarge_krylov_subspace", "Split the initial right-hand side into multiple vectors", "KSPHPDDM", i, &i, NULL, 1, std::numeric_limits<unsigned short>::max() - 1);
 45:       data->scntl[1 + (data->cntl[5] != HPDDM_KRYLOV_METHOD_BFBCG)] = i;
 46:     } else data->scntl[2] = 0;
 47:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 48:       i = HPDDM_ORTHOGONALIZATION_CGS;
 49:       PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, ALEN(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
 50:       j = HPDDM_QR_CHOLQR;
 51:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
 52:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 53:       i = PetscMin(40, ksp->max_it - 1);
 54:       PetscOptionsRangeInt("-ksp_gmres_restart", "Maximum number of Arnoldi vectors generated per cycle", "KSPHPDDM", i, &i, NULL, PetscMin(1, ksp->max_it), PetscMin(ksp->max_it, std::numeric_limits<unsigned short>::max() - 1));
 55:       data->scntl[1] = i;
 56:     }
 57:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG) {
 58:       j = HPDDM_QR_CHOLQR;
 59:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
 60:       data->cntl[1] = j;
 61:     }
 62:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 63:       i = PetscMin(20, data->scntl[1] - 1);
 64:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[1] - 1);
 65:       data->icntl[0] = i;
 66:       i = HPDDM_RECYCLE_TARGET_SM;
 67:       PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, ALEN(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
 68:       data->cntl[3] = i;
 69:       i = HPDDM_RECYCLE_STRATEGY_A;
 70:       PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, ALEN(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
 71:       data->cntl[4] = i;
 72:     }
 73:   }
 74:   PetscOptionsTail();
 75:   return(0);
 76: }

 78: static PetscErrorCode KSPView_HPDDM(KSP ksp, PetscViewer viewer)
 79: {
 80:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
 81:   HPDDM::PETScOperator *op = data->op;
 82:   const PetscScalar    *array = op ? op->storage() : NULL;
 83:   PetscBool            ascii;
 84:   PetscErrorCode       ierr;

 87:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
 88:   if (op && ascii) {
 89:     PetscViewerASCIIPrintf(viewer, "HPDDM type: %s\n", HPDDMType[static_cast<PetscInt>(data->cntl[5])]);
 90:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 91:       PetscViewerASCIIPrintf(viewer, "deflation subspace attached? %s\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
 92:       PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]);
 93:     }
 94:   }
 95:   return(0);
 96: }

 98: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
 99: {
100:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
101:   Mat            A;
102:   PetscInt       n;

106:   KSPGetOperators(ksp, &A, NULL);
107:   MatGetLocalSize(A, &n, NULL);
108:   data->op = new HPDDM::PETScOperator(ksp, n, 1);
109:   return(0);
110: }

112: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
113: {
114:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;
116:   if (data->op) {
117:     delete data->op;
118:     data->op = NULL;
119:   }
120:   return(0);
121: }

123: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
124: {

128:   KSPReset_HPDDM(ksp);
129:   KSPDestroyDefault(ksp);
130:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", NULL);
131:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", NULL);
132:   return(0);
133: }

135: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
136: {
137:   KSP_HPDDM         *data = (KSP_HPDDM*)ksp->data;
138:   PetscScalar       *x;
139:   const PetscScalar *b;
140:   MPI_Comm          comm;
141:   PetscErrorCode    ierr;

144:   PetscCitationsRegister(hpddmCitationKSP, &citeKSP);
145:   VecGetArray(ksp->vec_sol, &x);
146:   VecGetArrayRead(ksp->vec_rhs, &b);
147:   PetscObjectGetComm((PetscObject)ksp, &comm);
148:   ksp->its = HPDDM::IterativeMethod::solve(*data->op, b, x, 1, comm);
149:   VecRestoreArrayRead(ksp->vec_rhs, &b);
150:   VecRestoreArray(ksp->vec_sol, &x);
151:   if (ksp->its < ksp->max_it) ksp->reason = KSP_CONVERGED_RTOL;
152:   else ksp->reason = KSP_DIVERGED_ITS;
153:   return(0);
154: }

156: /*@
157:      KSPHPDDMSetDeflationSpace - Sets the deflation space used by Krylov methods with recycling. This space is viewed as a set of vectors stored in a MATDENSE (column major).

159:    Input Parameters:
160: +     ksp - iterative context
161: -     U - deflation space to be used during KSPSolve()

163:    Level: intermediate

165: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetDeflationSpace()
166: @*/
167: PetscErrorCode KSPHPDDMSetDeflationSpace(KSP ksp, Mat U)
168: {

175:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationSpace_C", (KSP, Mat), (ksp, U));
176:   return(0);
177: }

179: /*@
180:      KSPHPDDMGetDeflationSpace - Gets the deflation space computed by Krylov methods with recycling or NULL if KSPSolve() has not been called yet. This space is viewed as a set of vectors stored in a MATDENSE (column major). It is the responsibility of the user to free the returned Mat.

182:    Input Parameter:
183: .     ksp - iterative context

185:    Output Parameter:
186: .     U - deflation space generated during KSPSolve()

188:    Level: intermediate

190: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetDeflationSpace()
191: @*/
192: PetscErrorCode KSPHPDDMGetDeflationSpace(KSP ksp, Mat *U)
193: {

198:   PetscUseMethod(ksp, "KSPHPDDMGetDeflationSpace_C", (KSP, Mat*), (ksp, U));
199:   return(0);
200: }

202: static PetscErrorCode KSPHPDDMSetDeflationSpace_HPDDM(KSP ksp, Mat U)
203: {
204:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
205:   HPDDM::PETScOperator *op = data->op;
206:   Mat                  A, local;
207:   const PetscScalar    *array;
208:   PetscScalar          *copy;
209:   PetscInt             m1, M1, m2, M2, n2, N2, ldu;
210:   PetscBool            match;
211:   PetscErrorCode       ierr;

214:   KSPGetOperators(ksp, &A, NULL);
215:   MatGetLocalSize(A, &m1, NULL);
216:   MatGetLocalSize(U, &m2, &n2);
217:   MatGetSize(A, &M1, NULL);
218:   MatGetSize(U, &M2, &N2);
219:   if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a deflation space with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
220:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
221:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
222:   MatDenseGetLocalMatrix(U, &local);
223:   MatDenseGetArrayRead(local, &array);
224:   copy = op->allocate(m2, 1, N2);
225:   if (!copy) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
226:   MatDenseGetLDA(local, &ldu);
227:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
228:   MatDenseRestoreArrayRead(local, &array);
229:   return(0);
230: }

232: static PetscErrorCode KSPHPDDMGetDeflationSpace_HPDDM(KSP ksp, Mat *U)
233: {
234:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
235:   HPDDM::PETScOperator *op = data->op;
236:   Mat                  A;
237:   const PetscScalar    *array = op->storage();
238:   PetscScalar          *copy;
239:   PetscInt             m1, M1, N2 = op->k();
240:   PetscErrorCode       ierr;

243:   if (!array) *U = NULL;
244:   else {
245:     KSPGetOperators(ksp, &A, NULL);
246:     MatGetLocalSize(A, &m1, NULL);
247:     MatGetSize(A, &M1, NULL);
248:     MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U);
249:     MatDenseGetArray(*U, &copy);
250:     PetscArraycpy(copy, array, m1 * N2);
251:     MatDenseRestoreArray(*U, &copy);
252:   }
253:   return(0);
254: }

256: /*MC
257:      KSPHPDDM - Interface with the HPDDM library.

259:    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]. A chronological bibliography of relevant publications linked with KSP available in HPDDM through KSPHPDDM, and not available directly in PETSc, may be found below.

261:    Options Database Keys:
262: +   -ksp_richardson_scale <scale, default=1.0> - see KSPRICHARDSON
263: .   -ksp_gmres_restart <restart, default=40> - see KSPGMRES
264: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, or bfbcg
265: .   -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)
266: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
267: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
268: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
269: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible
270: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
271: .   -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)
272: -   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)

274:    References:
275: +   1980 - The Block Conjugate Gradient Algorithm and Related Methods. O'Leary. Linear Algebra and its Applications.
276: .   2006 - Recycling Krylov Subspaces for Sequences of Linear Systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
277: .   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.
278: .   2016 - Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers. Jolivet and Tournier. SC16.
279: -   2017 - A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.

281:    Level: intermediate

283: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
284: M*/
285: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
286: {
287:   KSP_HPDDM      *data;

291:   PetscNewLog(ksp, &data);
292:   ksp->data = (void*)data;
293:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
294:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
295:   ksp->ops->setup          = KSPSetUp_HPDDM;
296:   ksp->ops->solve          = KSPSolve_HPDDM;
297:   ksp->ops->reset          = KSPReset_HPDDM;
298:   ksp->ops->destroy        = KSPDestroy_HPDDM;
299:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
300:   ksp->ops->view           = KSPView_HPDDM;
301:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", KSPHPDDMSetDeflationSpace_HPDDM);
302:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", KSPHPDDMGetDeflationSpace_HPDDM);
303:   return(0);
304: }