Actual source code: hpddm.cxx

petsc-master 2019-11-16
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/matimpl.h>
  3:  #include <petsc/private/petschpddm.h>
  4: #include <petsc/private/pcimpl.h> /* this must be included after petschpddm.h so that _PCIMPL_H is not defined            */
  5:                                   /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */

  7: static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar>* const, IS const, IS, const Mat, Mat, std::vector<Vec>, PetscInt* const, PC_HPDDM_Level** const) = NULL;

  9: static PetscBool citePC = PETSC_FALSE;
 10: static const char hpddmCitationPC[] = "@inproceedings{jolivet2013scalable,\n\tTitle = {{Scalable Domain Decomposition Preconditioners For Heterogeneous Elliptic Problems}},\n\tAuthor = {Jolivet, Pierre and Hecht, Fr\'ed\'eric and Nataf, Fr\'ed\'eric and Prud'homme, Christophe},\n\tOrganization = {ACM},\n\tYear = {2013},\n\tSeries = {SC13},\n\tBooktitle = {Proceedings of the 2013 International Conference for High Performance Computing, Networking, Storage and Analysis}\n}\n";

 12: static const char *PCHPDDMCoarseCorrectionTypes[] = { "deflated", "additive", "balanced" };

 14: static PetscErrorCode PCReset_HPDDM(PC pc)
 15: {
 16:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;
 17:   PetscInt       i;

 21:   if (data->levels) {
 22:     for (i = 0; i < PETSC_HPDDM_MAXLEVELS && data->levels[i]; ++i) {
 23:       KSPDestroy(&data->levels[i]->ksp);
 24:       PCDestroy(&data->levels[i]->pc);
 25:       PetscFree(data->levels[i]);
 26:     }
 27:     PetscFree(data->levels);
 28:   }

 30:   ISDestroy(&data->is);
 31:   MatDestroy(&data->aux);
 32:   data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED;
 33:   data->Neumann    = PETSC_FALSE;
 34:   data->setup      = NULL;
 35:   data->setup_ctx  = NULL;
 36:   return(0);
 37: }

 39: static PetscErrorCode PCDestroy_HPDDM(PC pc)
 40: {
 41:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;

 45:   PCReset_HPDDM(pc);
 46:   PetscFree(data);
 47:   PetscObjectChangeTypeName((PetscObject)pc, 0);
 48:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", NULL);
 49:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", NULL);
 50:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", NULL);
 51:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", NULL);
 52:   return(0);
 53: }

 55: static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void* setup_ctx)
 56: {
 57:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;

 61:   if (is) {
 62:     PetscObjectReference((PetscObject)is);
 63:     if (data->is) { /* new overlap definition resets the PC */
 64:       PCReset_HPDDM(pc);
 65:       pc->setfromoptionscalled = 0;
 66:     }
 67:     ISDestroy(&data->is);
 68:     data->is = is;
 69:   }
 70:   if (A) {
 71:     PetscObjectReference((PetscObject)A);
 72:     MatDestroy(&data->aux);
 73:     data->aux = A;
 74:   }
 75:   if (setup) {
 76:     data->setup = setup;
 77:     data->setup_ctx = setup_ctx;
 78:   }
 79:   return(0);
 80: }

 82: /*MC
 83:      PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by PCHPDDM for the concurrent GenEO eigenproblems at the finest level. As an example, in a finite element context with nonoverlapping subdomains plus (overlapping) ghost elements, this could be the unassembled (Neumann) local overlapping operator. As opposed to the assembled (Dirichlet) local overlapping operator obtained by summing neighborhood contributions at the interface of ghost elements.

 85:    Input Parameters:
 86: +     pc - preconditioner context
 87: .     is - index set of the local auxiliary, e.g., Neumann, matrix
 88: .     A - auxiliary sequential matrix
 89: .     setup - function for generating the auxiliary matrix
 90: -     setup_ctx - context for setup

 92:    Level: intermediate

 94: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS
 95: M*/
 96: PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void* setup_ctx)
 97: {

104:   PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode (*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void*), (pc, is, A, setup, setup_ctx));
105:   return(0);
106: }

108: static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has)
109: {
110:   PC_HPDDM *data = (PC_HPDDM*)pc->data;

113:   data->Neumann = has;
114:   return(0);
115: }

117: /*MC
118:      PCHPDDMHasNeumannMat - Informs PCHPDDM that the Mat passed to PCHPDDMSetAuxiliaryMat is the local Neumann matrix. This may be used to bypass a call to MatCreateSubMatrices and to MatConvert for MATMPISBAIJ matrices. If a DMCreateNeumannOverlap implementation is available in the DM attached to the Pmat, or the Amat, or the PC, the flag is internally set to PETSC_TRUE. Its default value is otherwise PETSC_FALSE.

120:    Input Parameters:
121: +     pc - preconditioner context
122: -     has - Boolean value

124:    Level: intermediate

126: .seealso:  PCHPDDM, PCHPDDMSetAuxiliaryMat()
127: M*/
128: PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has)
129: {

134:   PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has));
135:   return(0);
136: }

138: static PetscErrorCode PCSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, PC pc)
139: {
140:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;
141:   PC_HPDDM_Level **levels = data->levels;
142:   MPI_Comm       comm;
143:   char           prefix[256];
144:   int            i = 1;
145:   PetscMPIInt    size, previous;
146:   PetscInt       n;
147:   PetscBool      flg = PETSC_TRUE;

151:   if (!data->levels) {
152:     PetscCalloc1(PETSC_HPDDM_MAXLEVELS, &levels);
153:     data->levels = levels;
154:   }
155:   PetscOptionsHead(PetscOptionsObject, "PCHPDDM options");
156:   PetscObjectGetComm((PetscObject)pc, &comm);
157:   MPI_Comm_size(comm, &size);
158:   previous = size;
159:   while (i < PETSC_HPDDM_MAXLEVELS) {
160:     PetscInt p = 1;

162:     if (!data->levels[i - 1]) {
163:       PetscNewLog(pc, data->levels + i - 1);
164:     }
165:     data->levels[i - 1]->parent = data;
166:     /* if the previous level has a single process, it is not possible to coarsen further */
167:     if (previous == 1 || !flg) break;
168:     data->levels[i - 1]->nu = 0;
169:     data->levels[i - 1]->threshold = -1.0;
170:     PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
171:     PetscOptionsInt(prefix, "Local number of deflation vectors computed by SLEPc", "none", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, NULL);
172:     PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
173:     PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "none", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, NULL);
174:     /* if there is no prescribed coarsening, just break out of the loop */
175:     if (data->levels[i - 1]->threshold <= 0.0 && data->levels[i - 1]->nu <= 0) break;
176:     else {
177:       ++i;
178:       PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
179:       PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
180:       if (!flg) {
181:         PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
182:         PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
183:       }
184:       if (flg) {
185:         /* if there are coarsening options for the next level, then register it  */
186:         /* otherwise, don't to avoid having both options levels_N_p and coarse_p */
187:         PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i);
188:         PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "none", p, &p, &flg, 1, PetscMax(1, previous / 2));
189:         previous = p;
190:       }
191:     }
192:   }
193:   data->N = i;
194:   n = 1;
195:   if (i > 1) {
196:     PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p");
197:     PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "none", n, &n, NULL, 1, PetscMax(1, previous / 2));
198:     PetscOptionsEList("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, 3, PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_DEFLATED], &n, &flg);
199:     if (flg) data->correction = PCHPDDMCoarseCorrectionType(n);
200:     PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann");
201:     PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", data->Neumann, &data->Neumann, NULL);
202:   }
203:   PetscOptionsTail();
204:   while (i < PETSC_HPDDM_MAXLEVELS && data->levels[i]) {
205:     PetscFree(data->levels[i++]);
206:   }
207:   return(0);
208: }

210: static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y)
211: {
212:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;

216:   PetscCitationsRegister(hpddmCitationPC, &citePC);
217:   if (data->levels[0]->ksp) {
218:     KSPSolve(data->levels[0]->ksp, x, y);
219:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM");
220:   return(0);
221: }

223: /*
224:      PCHPDDMGetComplexities - Computes the grid and operator complexities.

226:    Input Parameter:
227: .     pc - preconditioner context

229:    Output Parameters:
230: +     gc - grid complexity = sum_i(m_i) / m_1
231: .     oc - operator complexity = sum_i(nnz_i) / nnz_1

233:    Notes:
234:      PCGAMG does not follow the usual convention and names the grid complexity what is usually referred to as the operator complexity. PCHPDDM follows what is found in the literature, and in particular, what you get with PCHYPRE and -pc_hypre_boomeramg_print_statistics.

236:    Level: advanced

238: .seealso:  PCMGGetGridComplexity(), PCHPDDM
239: */
240: static PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc)
241: {
242:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;
243:   MatInfo        info;
244:   PetscInt       n, m;
245:   PetscLogDouble accumulate[2] { }, nnz1 = 1.0, m1 = 1.0;

249:   for (n = 0, *gc = 0, *oc = 0; n < data->N; ++n) {
250:     if (data->levels[n]->ksp) {
251:       Mat P;
252:       KSPGetOperators(data->levels[n]->ksp, NULL, &P);
253:       MatGetSize(P, &m, NULL);
254:       accumulate[0] += m;
255:       MatGetInfo(P, MAT_GLOBAL_SUM, &info);
256:       accumulate[1] += info.nz_used;
257:       if (n == 0) {
258:         m1 = m;
259:         nnz1 = info.nz_used;
260:       }
261:     }
262:   }
263:   *gc = static_cast<PetscReal>(accumulate[0]/m1);
264:   *oc = static_cast<PetscReal>(accumulate[1]/nnz1);
265:   return(0);
266: }

268: static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer)
269: {
270:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;
271:   PetscViewer    subviewer;
272:   PetscSubcomm   subcomm;
273:   PetscReal      oc, gc;
274:   PetscInt       i, tabs;
275:   PetscMPIInt    size, color, rank;
276:   PetscBool      ascii;

280:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
281:   if (ascii) {
282:     PetscViewerASCIIPrintf(viewer, "level%s: %D\n", data->N > 1 ? "s" : "", data->N);
283:     PCHPDDMGetComplexities(pc, &gc, &oc);
284:     if (data->N > 1) {
285:       PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", data->Neumann ? "PETSC_TRUE" : "PETSC_FALSE");
286:       PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction]);
287:       PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "");
288:       PetscViewerASCIIGetTab(viewer, &tabs);
289:       PetscViewerASCIISetTab(viewer, 0);
290:       for (i = 1; i < data->N; ++i) {
291:         PetscViewerASCIIPrintf(viewer, " %D", data->levels[i - 1]->nu);
292:         if (data->levels[i - 1]->threshold > -0.1) {
293:           PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold);
294:         }
295:       }
296:       PetscViewerASCIIPrintf(viewer, "\n");
297:       PetscViewerASCIISetTab(viewer, tabs);
298:     }
299:     PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc);
300:     if (data->levels[0]->ksp) {
301:       KSPView(data->levels[0]->ksp, viewer);
302:       if (data->levels[0]->pc) {
303:         PCView(data->levels[0]->pc, viewer);
304:       }
305:       for (i = 1; i < data->N; ++i) {
306:         if (data->levels[i]->ksp) color = 1;
307:         else color = 0;
308:         MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
309:         MPI_Comm_size(PetscObjectComm((PetscObject)pc), &rank);
310:         PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm);
311:         PetscSubcommSetNumber(subcomm, PetscMin(size, 2));
312:         PetscSubcommSetTypeGeneral(subcomm, color, rank);
313:         PetscViewerASCIIPushTab(viewer);
314:         PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
315:         if (color == 1) {
316:           KSPView(data->levels[i]->ksp, subviewer);
317:           if (data->levels[i]->pc) {
318:             PCView(data->levels[i]->pc, subviewer);
319:           }
320:           PetscViewerFlush(subviewer);
321:         }
322:         PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
323:         PetscViewerASCIIPopTab(viewer);
324:         PetscSubcommDestroy(&subcomm);
325:         PetscViewerFlush(viewer);
326:       }
327:     }
328:   }
329:   return(0);
330: }

332: static PetscErrorCode PCHPDDMShellSetUp(PC pc)
333: {
334:   PC_HPDDM_Level *ctx;
335:   Mat            A, P;
336:   Vec            x;
337:   const char     *pcpre;

341:   PCShellGetContext(pc, (void**)&ctx);
342:   KSPGetOptionsPrefix(ctx->ksp, &pcpre);
343:   KSPGetOperators(ctx->ksp, &A, &P);
344:   /* smoother */
345:   PCSetOptionsPrefix(ctx->pc, pcpre);
346:   PCSetOperators(ctx->pc, A, P);
347:   if (!ctx->v[0]) {
348:     VecDuplicateVecs(ctx->D, 1, &ctx->v[0]);
349:     if (!std::is_same<PetscScalar, PetscReal>::value) {
350:       VecDestroy(&ctx->D);
351:     }
352:     MatCreateVecs(A, &x, NULL);
353:     VecDuplicateVecs(x, 2, &ctx->v[1]);
354:     VecDestroy(&x);
355:   }
356:   return(0);
357: }

359: PETSC_STATIC_INLINE PetscErrorCode PCHPDDMDeflate_Private(PC pc, Vec x, Vec y)
360: {
361:   PC_HPDDM_Level *ctx;
362:   PetscScalar    *out;

366:   PCShellGetContext(pc, (void**)&ctx);
367:   /* going from PETSc to HPDDM numbering */
368:   VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
369:   VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
370:   VecGetArray(ctx->v[0][0], &out);
371:   ctx->P->deflation<false>(NULL, out, 1); /* y = Q x */
372:   VecRestoreArray(ctx->v[0][0], &out);
373:   /* going from HPDDM to PETSc numbering */
374:   VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
375:   VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
376:   return(0);
377: }

379: /*MC
380:      PCHPDDMShellApply - Applies a (2) deflated, (1) additive, or (3) balanced coarse correction. In what follows, E = Z Pmat Z^T and Q = Z^T E^-1 Z.

382: .vb
383:    (1) y =                Pmat^-1              x + Q x,
384:    (2) y =                Pmat^-1 (I - Amat Q) x + Q x (default),
385:    (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x.
386: .ve

388:    Input Parameters:
389: +     pc - preconditioner context
390: .     x - input vector
391: -     y - output vector

393:    Application Interface Routine: PCApply()

395:    Notes:
396:      The options of Pmat^1 = pc(Pmat) are prefixed by -pc_hpddm_levels_1_pc_. Z is a tall-and-skiny matrix assembled by HPDDM. The number of processes on which (Z Pmat Z^T) is aggregated is set via -pc_hpddm_coarse_p.
397:      The options of (Z Pmat Z^T)^-1 = ksp(Z Pmat Z^T) are prefixed by -pc_hpddm_coarse_ (KSPPREONLY and PCCHOLESKY by default), unless a multilevel correction is turned on, in which case, this function is called recursively at each level except the coarsest one.
398:      (1) and (2) visit the "next" level (in terms of coarsening) once per application, while (3) visits it twice, so it is asymptotically twice costlier. (2) is not symmetric even if both Amat and Pmat are symmetric.

400: .seealso:  PCHPDDM, PCHPDDMCoarseCorrectionType
401: M*/
402: static PetscErrorCode PCHPDDMShellApply(PC pc, Vec x, Vec y)
403: {
404:   PC_HPDDM_Level *ctx;
405:   Mat            A;

409:   PCShellGetContext(pc, (void**)&ctx);
410:   if (ctx->P) {
411:     KSPGetOperators(ctx->ksp, &A, NULL);
412:     PCHPDDMDeflate_Private(pc, x, y);                    /* y = Q x                          */
413:     if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
414:       MatMult(A, y, ctx->v[1][0]);                       /* y = A Q x                        */
415:       VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x);     /* y = (I - A Q) x                  */
416:       PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0]);      /* y = M^-1 (I - A Q) x             */
417:       if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
418:         MatMultTranspose(A, ctx->v[1][0], ctx->v[1][1]); /* z = A^T M^-1 (I - A Q) x         */
419:         PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1]);
420:         VecAXPY(ctx->v[1][0], -1.0, ctx->v[1][1]);       /* y = (I - Q A^T) M^-1 (I - A Q) x */
421:       }
422:       VecAXPY(y, 1.0, ctx->v[1][0]);                     /* y = y + Q x                      */
423:     } else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE) {
424:       PCApply(ctx->pc, x, ctx->v[1][0]);
425:       VecAXPY(y, 1.0, ctx->v[1][0]);                     /* y = M^-1 x + Q x                 */
426:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
427:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
428:   return(0);
429: }

431: static PetscErrorCode PCHPDDMShellDestroy(PC pc)
432: {
433:   PC_HPDDM_Level *ctx;

437:   PCShellGetContext(pc, (void**)&ctx);
438:   HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE);
439:   VecDestroyVecs(1, &ctx->v[0]);
440:   VecDestroyVecs(2, &ctx->v[1]);
441:   VecDestroy(&ctx->D);
442:   VecScatterDestroy(&ctx->scatter);
443:   PCDestroy(&ctx->pc);
444:   return(0);
445: }

447: static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc)
448: {
449:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;

453:   if (data->setup) {
454:     Mat       P;
455:     Vec       x, xt = NULL;
456:     PetscReal t = 0.0, s = 0.0;

458:     PCGetOperators(pc, NULL, &P);
459:     PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject*)&x);
460:     PetscStackPush("PCHPDDM Neumann callback");
461:     (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx);
462:     PetscStackPop;
463:   }
464:   return(0);
465: }

467: static PetscErrorCode PCSetUp_HPDDM(PC pc)
468: {
469:   PC_HPDDM                 *data = (PC_HPDDM*)pc->data;
470:   PC                       inner;
471:   Mat                      *sub, A, P, N, C, uaux = NULL;
472:   Vec                      xin, v;
473:   std::vector<Vec>         initial;
474:   IS                       is[1], loc, uis = data->is;
475:   ISLocalToGlobalMapping   l2g;
476:   MPI_Comm                 comm;
477:   char                     prefix[256];
478:   const char               *pcpre;
479:   const PetscScalar* const *ev;
480:   PetscInt                 n, requested = data->N;
481:   PetscBool                subdomains = PETSC_FALSE, flag = PETSC_FALSE, ismatis;
482:   DM                       dm;
483:   PetscErrorCode           ierr;

486:   if (!data->levels[0]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated");
487:   PCGetOptionsPrefix(pc, &pcpre);
488:   PCGetOperators(pc, &A, &P);
489:   PetscObjectGetComm((PetscObject)pc, &comm);
490:   if (!data->levels[0]->ksp) {
491:     KSPCreate(comm, &data->levels[0]->ksp);
492:     PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse");
493:     KSPSetOptionsPrefix(data->levels[0]->ksp, prefix);
494:     KSPSetType(data->levels[0]->ksp, KSPPREONLY);
495:   } else {
496:     const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0;

498:     if (addr != &HPDDM::i__0) {
499:       /* reuse previously computed eigenvectors */
500:       ev = data->levels[0]->P->getVectors();
501:       if (ev) {
502:         initial.reserve(*addr);
503:         VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin);
504:         for (n = 0; n < *addr; ++n) {
505:           VecDuplicate(xin, &v);
506:           VecPlaceArray(xin, ev[n]);
507:           VecCopy(xin, v);
508:           initial.emplace_back(v);
509:           VecResetArray(xin);
510:         }
511:         VecDestroy(&xin);
512:       }
513:     }
514:     /* reset coarser levels */
515:     for (n = 1; n < PETSC_HPDDM_MAXLEVELS && data->levels[n]; ++n) {
516:       KSPDestroy(&data->levels[n]->ksp);
517:       PCDestroy(&data->levels[n]->pc);
518:     }
519:   }
520:   KSPSetOperators(data->levels[0]->ksp, A, P);

522:   PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis);
523:   if (!data->is && !ismatis) {
524:     PetscErrorCode (*create)(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void**) = NULL;
525:     PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*) = NULL;
526:     void           *uctx = NULL;

528:     /* first see if we can get the data from the DM */
529:     MatGetDM(P, &dm);
530:     if (!dm) {
531:       MatGetDM(A, &dm);
532:     }
533:     if (!dm) {
534:       PCGetDM(pc, &dm);
535:     }
536:     if (dm) { /* this is the hook for DMPLEX and DMDA for which the auxiliary Mat is the local Neumann matrix */
537:       PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create);
538:       if (create) {
539:         (*create)(dm, &uis, &uaux, &usetup, &uctx);
540:         data->Neumann = PETSC_TRUE;
541:       }
542:     }
543:     if (!create) {
544:       if (!uis) {
545:         PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject*)&uis);
546:         PetscObjectReference((PetscObject)uis);
547:       }
548:       if (!uaux) {
549:         PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject*)&uaux);
550:         PetscObjectReference((PetscObject)uaux);
551:       }
552:     }
553:     PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx);
554:     MatDestroy(&uaux);
555:     ISDestroy(&uis);
556:   }

558:   if (!ismatis) {
559:     PCHPDDMSetUpNeumannOverlap_Private(pc);
560:   }

562:   if (data->is || (ismatis && data->N > 1)) {
563:     if (ismatis) {
564:       std::initializer_list<std::string> list = { MATSEQBAIJ, MATSEQSBAIJ };
565:       MatISGetLocalMat(P, &N);
566:       std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name);
567:       MatISRestoreLocalMat(P, &N);
568:       switch (std::distance(list.begin(), it)) {
569:       case 0:
570:         MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
571:         break;
572:       case 1:
573:         /* MatCreateSubMatrices does not work with [Seq|MPI]SBAIJ and unsorted ISes, so convert to MPIBAIJ */
574:         MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
575:         MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE);
576:         break;
577:       default:
578:         MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C);
579:       }
580:       MatGetLocalToGlobalMapping(P, &l2g, NULL);
581:       PetscObjectReference((PetscObject)P);
582:       KSPSetOperators(data->levels[0]->ksp, A, C);
583:       std::swap(C, P);
584:       ISLocalToGlobalMappingGetSize(l2g, &n);
585:       ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc);
586:       ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0]);
587:       ISDestroy(&loc);
588:       /* the auxiliary Mat is _not_ the local Neumann matrix */
589:       /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap */
590:       data->Neumann = PETSC_FALSE;
591:     } else {
592:       is[0] = data->is;
593:       PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_define_subdomains", &subdomains, NULL);
594:       PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_has_neumann", &data->Neumann, NULL);
595:       ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc);
596:     }
597:     if (data->N > 1 && (data->aux || ismatis)) {
598:       if (!loadedSym) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level");
599:       MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
600:       if (ismatis) {
601:         /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */
602:         MatIncreaseOverlap(P, 1, is, 1);
603:         ISDestroy(&data->is);
604:         data->is = is[0];
605:       } else {
606: #if defined(PETSC_USE_DEBUG)
607:         PetscBool equal;
608:         IS        intersect;

610:         ISIntersect(data->is, loc, &intersect);
611:         ISEqualUnsorted(loc, intersect, &equal);
612:         ISDestroy(&intersect);
613:         if (!equal) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS of the auxiliary Mat does not include all local rows of A");
614: #endif
615:         if (!data->Neumann) {
616:           PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flag);
617:           if (flag) {
618:             /* maybe better to ISSort(is[0]), MatCreateSubMatrices, and then MatPermute */
619:             /* but there is no MatPermute_SeqSBAIJ, so as before, just use MATMPIBAIJ   */
620:             MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux);
621:             flag = PETSC_FALSE;
622:           }
623:         }
624:       }
625:       if (!uaux) {
626:         if (data->Neumann) sub = &data->aux;
627:         else {
628:           MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub);
629:         }
630:       } else {
631:         MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub);
632:         MatDestroy(&uaux);
633:         MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub);
634:       }
635:       /* Vec holding the partition of unity */
636:       if (!data->levels[0]->D) {
637:         ISGetLocalSize(data->is, &n);
638:         VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D);
639:       }
640:       if (!data->levels[0]->scatter) {
641:         MatCreateVecs(P, &xin, NULL);
642:         if (ismatis) {
643:           MatDestroy(&P);
644:         }
645:         VecScatterCreate(xin, data->is, data->levels[0]->D, NULL, &data->levels[0]->scatter);
646:         VecDestroy(&xin);
647:       }
648:       if (data->levels[0]->P) {
649:         /* if the pattern is the same and PCSetUp has previously succeeded, reuse HPDDM buffers and connectivity */
650:         HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE);
651:       }
652:       if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>();
653:       (*loadedSym)(data->levels[0]->P, loc, data->is, sub[0], ismatis ? C : data->aux, initial, &data->N, data->levels);
654:       if (!data->Neumann)
655:         MatDestroySubMatrices(1, &sub);
656:       if (ismatis) data->is = NULL;
657:       for (n = 0; n < data->N - 1; ++n) {
658:         if (data->levels[n]->P) {
659:           PC spc;

661:           /* force the PC to be PCSHELL to do the coarse grid corrections */
662:           KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE);
663:           KSPGetPC(data->levels[n]->ksp, &spc);
664:           PCSetType(spc, PCSHELL);
665:           PCShellSetContext(spc, data->levels[n]);
666:           PCShellSetSetUp(spc, PCHPDDMShellSetUp);
667:           PCShellSetApply(spc, PCHPDDMShellApply);
668:           PCShellSetDestroy(spc, PCHPDDMShellDestroy);
669:           if (!data->levels[n]->pc) {
670:             PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc);
671:           }
672:           PCSetUp(spc);
673:         }
674:       }
675:     } else flag = PETSC_TRUE;
676:     if (!ismatis && subdomains) {
677:       if (flag) {
678:         KSPGetPC(data->levels[0]->ksp, &inner);
679:       } else inner = data->levels[0]->pc;
680:       PCSetType(inner, PCASM);
681:       PCASMSetLocalSubdomains(inner, 1, is, &loc);
682:     }
683:     ISDestroy(&loc);
684:   } else data->N = 1; /* enforce this value to 1 if there is no way to build another level */
685:   if (requested != data->N) {
686:     PetscInfo4(pc, "HPDDM: %D levels requested, only %D built. Options for level(s) > %D, including -%spc_hpddm_coarse_ will not be taken into account.\n", requested, data->N, data->N, pcpre ? pcpre : "");
687:     PetscInfo2(pc, "HPDDM: It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%D_eps_threshold so that at least one local deflation vector will be selected.\n", pcpre ? pcpre : "", data->N);
688:     /* cannot use PCHPDDMShellDestroy because PCSHELL not set for unassembled levels */
689:     for (n = data->N - 1; n < requested - 1; ++n) {
690:       if (data->levels[n]->P) {
691:         HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE);
692:         VecDestroyVecs(1, &data->levels[n]->v[0]);
693:         VecDestroyVecs(2, &data->levels[n]->v[1]);
694:         VecDestroy(&data->levels[n]->D);
695:         VecScatterDestroy(&data->levels[n]->scatter);
696:       }
697:     }
698: #if defined(PETSC_USE_DEBUG)
699:     SETERRQ6(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "%D levels requested, only %D built. Options for level(s) > %D, including -%spc_hpddm_coarse_ will not be taken into account. It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%D_eps_threshold so that at least one local deflation vector will be selected. If you don't want this to error out, compile --with-debugging=0", requested, data->N, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N);
700: #endif
701:   }

703:   /* these solvers are created after PCSetFromOptions is called */
704:   if (pc->setfromoptionscalled) {
705:     for (n = 0; n < data->N; ++n) {
706:       if (data->levels[n]->ksp) {
707:         KSPSetFromOptions(data->levels[n]->ksp);
708:       }
709:       if (data->levels[n]->pc) {
710:         PCSetFromOptions(data->levels[n]->pc);
711:       }
712:     }
713:     pc->setfromoptionscalled = 0;
714:   }
715:   return(0);
716: }

718: /*@C
719:      PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type.

721:    Input Parameters:
722: +     pc - preconditioner context
723: -     type - PC_HPDDM_COARSE_CORRECTION_DEFLATED, PC_HPDDM_COARSE_CORRECTION_ADDITIVE, or PC_HPDDM_COARSE_CORRECTION_BALANCED

725:    Options Database Key:
726: .   -pc_hpddm_coarse_correction <deflated, additive, balanced> - type of coarse correction to apply

728:    Level: intermediate

730: .seealso:  PCHPDDMGetCoarseCorrectionType(), PCHPDDM, PCHPDDMCoarseCorrectionType
731: @*/
732: PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type)
733: {

738:   PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type));
739:   return(0);
740: }

742: /*@C
743:      PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type.

745:    Input Parameter:
746: .     pc - preconditioner context

748:    Output Parameter:
749: -     type - PC_HPDDM_COARSE_CORRECTION_DEFLATED, PC_HPDDM_COARSE_CORRECTION_ADDITIVE, or PC_HPDDM_COARSE_CORRECTION_BALANCED

751:    Level: intermediate

753: .seealso:  PCHPDDMSetCoarseCorrectionType(), PCHPDDM, PCHPDDMCoarseCorrectionType
754: @*/
755: PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type)
756: {

761:   PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType*), (pc, type));
762:   return(0);
763: }

765: static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type)
766: {
767:   PC_HPDDM       *data = (PC_HPDDM*)pc->data;

770:   if (type < 0 || type > 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PCHPDDMCoarseCorrectionType %d", type);
771:   data->correction = type;
772:   return(0);
773: }

775: static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type)
776: {
777:   PC_HPDDM *data = (PC_HPDDM*)pc->data;

780:   *type = data->correction;
781:   return(0);
782: }

784: /*MC
785:      PCHPDDM - Interface with the HPDDM library.

787:    This PC may be used to build multilevel spectral domain decomposition methods based on the GenEO framework [2011, 2019]. It may be viewed as an alternative to spectral AMGE or PCBDDC with adaptive selection of constraints. Here is a chronological bibliography of relevant publications linked with PC available in HPDDM through PCHPDDM.

789: .vb
790:    [2011] A robust two-level domain decomposition preconditioner for systems of PDEs. Spillane, Dolean, Hauret, Nataf, Pechstein, and Scheichl. Comptes Rendus Mathematique.
791:    [2013] Scalable Domain Decomposition Preconditioners For Heterogeneous Elliptic Problems. Jolivet, Hecht, Nataf, and Prud'homme. SC13.
792:    [2015] An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation. Dolean, Jolivet, and Nataf. SIAM.
793:    [2019] A Multilevel Schwarz Preconditioner Based on a Hierarchy of Robust Coarse Spaces. Al Daas, Grigori, Jolivet, and Tournier.
794: .ve

796:    The matrix to be preconditioned (Pmat) may be unassembled (MATIS) or assembled (MATMPIAIJ, MATMPIBAIJ, or MATMPISBAIJ). For multilevel preconditioning, when using an assembled Pmat, one must provide an auxiliary local Mat (unassembled local operator for GenEO) using PCHPDDMSetAuxiliaryMat. Calling this routine is not needed when using a MATIS Pmat (assembly done internally using MatConvert).

798:    Options Database Keys:
799: +   -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls PCASMSetLocalSubdomains with the IS supplied in PCHPDDMSetAuxiliaryMat (only relevant with an assembled Pmat)
800: .   -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the PC that the local Neumann matrix is supplied in PCHPDDMSetAuxiliaryMat 
801: -   -pc_hpddm_coarse_correction <type, default=deflated> - determines the PCHPDDMCoarseCorrectionType when calling PCApply

803:    Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set with
804: .vb
805:       -pc_hpddm_levels_%d_pc_
806:       -pc_hpddm_levels_%d_ksp_
807:       -pc_hpddm_levels_%d_eps_
808:       -pc_hpddm_levels_%d_p
809:       -pc_hpddm_levels_%d_mat_type_
810:       -pc_hpddm_coarse_
811:       -pc_hpddm_coarse_p
812:       -pc_hpddm_coarse_mat_type_
813: .ve
814:    e.g., -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_levels_2_p 4 -pc_hpddm_levels_2_sub_pc_type lu -pc_hpddm_levels_2_eps_nev 10 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a MATMPIBAIJ (default is MATMPISBAIJ).

816:    In order to activate a "level N+1" coarse correction, it is mandatory to call -pc_hpddm_levels_N_eps_nev <nu> or -pc_hpddm_levels_N_eps_threshold <val>. The default -pc_hpddm_coarse_p value is 1, meaning that the coarse operator is aggregated on a single process.

818:    This preconditioner requires that you build PETSc with SLEPc (--download-slepc=1). By default, the underlying concurrent eigenproblems are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. [2011, 2013]. As stated above, SLEPc options are available through -pc_hpddm_levels_%d_, e.g., -pc_hpddm_levels_1_eps_type arpack -pc_hpddm_levels_1_eps_threshold 0.1 -pc_hpddm_levels_1_st_type sinvert.

820:    Level: intermediate

822: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCHPDDMSetAuxiliaryMat(), MATIS, PCBDDC, PCDEFLATION, PCTELESCOPE
823: M*/
824: PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc)
825: {
826:   PetscBool      found;
827:   char           lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN];
828:   PC_HPDDM       *data;

832:   if (!loadedSym) {
833:     PetscStrcpy(dir, "${PETSC_LIB_DIR}");
834:     PetscOptionsGetString(NULL, NULL, "-hpddm_dir", dir, sizeof(dir), NULL);
835:     PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir);
836:     PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, &found);
837:     if (found) {
838:       PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib);
839:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib);
840:     PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, NULL, "PCHPDDM_Internal", (void**)&loadedSym);
841:   }
842:   if (!loadedSym) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in %s", lib);
843:   PetscNewLog(pc, &data);
844:   pc->data                     = data;
845:   pc->ops->reset               = PCReset_HPDDM;
846:   pc->ops->destroy             = PCDestroy_HPDDM;
847:   pc->ops->setfromoptions      = PCSetFromOptions_HPDDM;
848:   pc->ops->setup               = PCSetUp_HPDDM;
849:   pc->ops->apply               = PCApply_HPDDM;
850:   pc->ops->view                = PCView_HPDDM;
851:   pc->ops->applytranspose      = 0;
852:   pc->ops->applysymmetricleft  = 0;
853:   pc->ops->applysymmetricright = 0;
854:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM);
855:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM);
856:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM);
857:   PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM);
858:   return(0);
859: }