Actual source code: hpddm.cxx
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/matimpl.h>
3: #include <petsc/private/petschpddm.h>
4: #include <petsc/private/pcimpl.h>
5: /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */
6: #if defined(PETSC_HAVE_FORTRAN)
7: #include <petsc/private/fortranimpl.h>
8: #endif
10: static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar>* const, IS, Mat, Mat, Mat, std::vector<Vec>, PC_HPDDM_Level** const) = NULL;
12: static PetscBool PCHPDDMPackageInitialized = PETSC_FALSE;
14: PetscLogEvent PC_HPDDM_Strc;
15: PetscLogEvent PC_HPDDM_PtAP;
16: PetscLogEvent PC_HPDDM_PtBP;
17: PetscLogEvent PC_HPDDM_Next;
18: PetscLogEvent PC_HPDDM_SetUp[PETSC_HPDDM_MAXLEVELS];
19: PetscLogEvent PC_HPDDM_Solve[PETSC_HPDDM_MAXLEVELS];
21: static const char *PCHPDDMCoarseCorrectionTypes[] = { "deflated", "additive", "balanced" };
23: static PetscErrorCode PCReset_HPDDM(PC pc)
24: {
25: PC_HPDDM *data = (PC_HPDDM*)pc->data;
26: PetscInt i;
30: if (data->levels) {
31: for (i = 0; i < PETSC_HPDDM_MAXLEVELS && data->levels[i]; ++i) {
32: KSPDestroy(&data->levels[i]->ksp);
33: PCDestroy(&data->levels[i]->pc);
34: PetscFree(data->levels[i]);
35: }
36: PetscFree(data->levels);
37: }
39: ISDestroy(&data->is);
40: MatDestroy(&data->aux);
41: MatDestroy(&data->B);
42: VecDestroy(&data->normal);
43: data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED;
44: data->Neumann = PETSC_FALSE;
45: data->setup = NULL;
46: data->setup_ctx = NULL;
47: return(0);
48: }
50: static PetscErrorCode PCDestroy_HPDDM(PC pc)
51: {
52: PC_HPDDM *data = (PC_HPDDM*)pc->data;
56: PCReset_HPDDM(pc);
57: PetscFree(data);
58: PetscObjectChangeTypeName((PetscObject)pc, NULL);
59: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", NULL);
60: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", NULL);
61: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", NULL);
62: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", NULL);
63: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", NULL);
64: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", NULL);
65: return(0);
66: }
68: static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void* setup_ctx)
69: {
70: PC_HPDDM *data = (PC_HPDDM*)pc->data;
74: if (is) {
75: PetscObjectReference((PetscObject)is);
76: if (data->is) { /* new overlap definition resets the PC */
77: PCReset_HPDDM(pc);
78: pc->setfromoptionscalled = 0;
79: }
80: ISDestroy(&data->is);
81: data->is = is;
82: }
83: if (A) {
84: PetscObjectReference((PetscObject)A);
85: MatDestroy(&data->aux);
86: data->aux = A;
87: }
88: if (setup) {
89: data->setup = setup;
90: data->setup_ctx = setup_ctx;
91: }
92: return(0);
93: }
95: /*@
96: 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.
98: Input Parameters:
99: + pc - preconditioner context
100: . is - index set of the local auxiliary, e.g., Neumann, matrix
101: . A - auxiliary sequential matrix
102: . setup - function for generating the auxiliary matrix
103: - setup_ctx - context for setup
105: Level: intermediate
107: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCHPDDMSetRHSMat(), MATIS
108: @*/
109: PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void* setup_ctx)
110: {
117: #if defined(PETSC_HAVE_FORTRAN)
118: if (reinterpret_cast<void*>(setup) == reinterpret_cast<void*>(PETSC_NULL_FUNCTION_Fortran)) setup = NULL;
119: if (setup_ctx == PETSC_NULL_INTEGER_Fortran) setup_ctx = NULL;
120: #endif
121: PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode (*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void*), (pc, is, A, setup, setup_ctx));
122: return(0);
123: }
125: static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has)
126: {
127: PC_HPDDM *data = (PC_HPDDM*)pc->data;
130: data->Neumann = has;
131: return(0);
132: }
134: /*@
135: 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.
137: Input Parameters:
138: + pc - preconditioner context
139: - has - Boolean value
141: Level: intermediate
143: .seealso: PCHPDDM, PCHPDDMSetAuxiliaryMat()
144: @*/
145: PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has)
146: {
151: PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has));
152: return(0);
153: }
155: static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B)
156: {
157: PC_HPDDM *data = (PC_HPDDM*)pc->data;
161: PetscObjectReference((PetscObject)B);
162: MatDestroy(&data->B);
163: data->B = B;
164: return(0);
165: }
167: /*@
168: PCHPDDMSetRHSMat - Sets the right-hand side matrix used by PCHPDDM for the concurrent GenEO eigenproblems at the finest level. Must be used in conjuction with PCHPDDMSetAuxiliaryMat(N), so that Nv = lambda Bv is solved using EPSSetOperators(N, B). It is assumed that N and B are provided using the same numbering. This provides a means to try more advanced methods such as GenEO-II or H-GenEO.
170: Input Parameters:
171: + pc - preconditioner context
172: - B - right-hand side sequential matrix
174: Level: advanced
176: .seealso: PCHPDDMSetAuxiliaryMat(), PCHPDDM
177: @*/
178: PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B)
179: {
184: if (B) {
186: PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B));
187: }
188: return(0);
189: }
191: static PetscErrorCode PCSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, PC pc)
192: {
193: PC_HPDDM *data = (PC_HPDDM*)pc->data;
194: PC_HPDDM_Level **levels = data->levels;
195: char prefix[256];
196: int i = 1;
197: PetscMPIInt size, previous;
198: PetscInt n;
199: PetscBool flg = PETSC_TRUE;
203: if (!data->levels) {
204: PetscCalloc1(PETSC_HPDDM_MAXLEVELS, &levels);
205: data->levels = levels;
206: }
207: PetscOptionsHead(PetscOptionsObject, "PCHPDDM options");
208: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
209: previous = size;
210: while (i < PETSC_HPDDM_MAXLEVELS) {
211: PetscInt p = 1;
213: if (!data->levels[i - 1]) {
214: PetscNewLog(pc, data->levels + i - 1);
215: }
216: data->levels[i - 1]->parent = data;
217: /* if the previous level has a single process, it is not possible to coarsen further */
218: if (previous == 1 || !flg) break;
219: data->levels[i - 1]->nu = 0;
220: data->levels[i - 1]->threshold = -1.0;
221: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
222: PetscOptionsInt(prefix, "Local number of deflation vectors computed by SLEPc", "none", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, NULL);
223: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
224: PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "none", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, NULL);
225: if (i == 1) {
226: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp");
227: PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "none", PETSC_FALSE, &data->share, NULL);
228: }
229: /* if there is no prescribed coarsening, just break out of the loop */
230: if (data->levels[i - 1]->threshold <= 0.0 && data->levels[i - 1]->nu <= 0) break;
231: else {
232: ++i;
233: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
234: PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
235: if (!flg) {
236: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
237: PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
238: }
239: if (flg) {
240: /* if there are coarsening options for the next level, then register it */
241: /* otherwise, don't to avoid having both options levels_N_p and coarse_p */
242: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i);
243: PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "none", p, &p, &flg, 1, PetscMax(1, previous / 2));
244: previous = p;
245: }
246: }
247: }
248: data->N = i;
249: n = 1;
250: if (i > 1) {
251: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p");
252: PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "none", n, &n, NULL, 1, PetscMax(1, previous / 2));
253: PetscOptionsEList("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, 3, PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_DEFLATED], &n, &flg);
254: if (flg) data->correction = PCHPDDMCoarseCorrectionType(n);
255: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann");
256: PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", data->Neumann, &data->Neumann, NULL);
257: data->log_separate = PETSC_FALSE;
258: if (PetscDefined(USE_LOG)) {
259: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate");
260: PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", NULL, data->log_separate, &data->log_separate, NULL);
261: }
262: }
263: PetscOptionsTail();
264: while (i < PETSC_HPDDM_MAXLEVELS && data->levels[i]) {
265: PetscFree(data->levels[i++]);
266: }
267: return(0);
268: }
270: static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y)
271: {
272: PC_HPDDM *data = (PC_HPDDM*)pc->data;
276: PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
277: if (data->levels[0]->ksp) {
278: if (data->log_separate) { /* coarser-level events are directly triggered in HPDDM */
279: PetscLogEventBegin(PC_HPDDM_Solve[0], data->levels[0]->ksp, 0, 0, 0);
280: }
281: KSPSolve(data->levels[0]->ksp, x, y);
282: if (data->log_separate) {
283: PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, 0, 0, 0);
284: }
285: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM");
286: return(0);
287: }
289: static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y)
290: {
291: PC_HPDDM *data = (PC_HPDDM*)pc->data;
295: PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
296: if (data->levels[0]->ksp) {
297: KSPMatSolve(data->levels[0]->ksp, X, Y);
298: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM");
299: return(0);
300: }
302: /*@C
303: PCHPDDMGetComplexities - Computes the grid and operator complexities.
305: Input Parameter:
306: . pc - preconditioner context
308: Output Parameters:
309: + gc - grid complexity = sum_i(m_i) / m_1
310: - oc - operator complexity = sum_i(nnz_i) / nnz_1
312: Notes:
313: 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.
315: Level: advanced
317: .seealso: PCMGGetGridComplexity(), PCHPDDM
318: @*/
319: static PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc)
320: {
321: PC_HPDDM *data = (PC_HPDDM*)pc->data;
322: MatInfo info;
323: PetscInt n, m;
324: PetscLogDouble accumulate[2] { }, nnz1 = 1.0, m1 = 1.0;
328: for (n = 0, *gc = 0, *oc = 0; n < data->N; ++n) {
329: if (data->levels[n]->ksp) {
330: Mat P;
331: KSPGetOperators(data->levels[n]->ksp, NULL, &P);
332: MatGetSize(P, &m, NULL);
333: accumulate[0] += m;
334: if (P->ops->getinfo) {
335: MatGetInfo(P, MAT_GLOBAL_SUM, &info);
336: accumulate[1] += info.nz_used;
337: }
338: if (n == 0) {
339: m1 = m;
340: if (P->ops->getinfo) nnz1 = info.nz_used;
341: }
342: }
343: }
344: *gc = static_cast<PetscReal>(accumulate[0]/m1);
345: *oc = static_cast<PetscReal>(accumulate[1]/nnz1);
346: return(0);
347: }
349: static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer)
350: {
351: PC_HPDDM *data = (PC_HPDDM*)pc->data;
352: PetscViewer subviewer;
353: PetscSubcomm subcomm;
354: PetscReal oc, gc;
355: PetscInt i, tabs;
356: PetscMPIInt size, color, rank;
357: PetscBool ascii;
361: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
362: if (ascii) {
363: PetscViewerASCIIPrintf(viewer, "level%s: %D\n", data->N > 1 ? "s" : "", data->N);
364: PCHPDDMGetComplexities(pc, &gc, &oc);
365: if (data->N > 1) {
366: PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[data->Neumann]);
367: PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction]);
368: PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "");
369: PetscViewerASCIIGetTab(viewer, &tabs);
370: PetscViewerASCIISetTab(viewer, 0);
371: for (i = 1; i < data->N; ++i) {
372: PetscViewerASCIIPrintf(viewer, " %D", data->levels[i - 1]->nu);
373: if (data->levels[i - 1]->threshold > -0.1) {
374: PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold);
375: }
376: }
377: PetscViewerASCIIPrintf(viewer, "\n");
378: PetscViewerASCIISetTab(viewer, tabs);
379: }
380: PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc);
381: if (data->levels[0]->ksp) {
382: KSPView(data->levels[0]->ksp, viewer);
383: if (data->levels[0]->pc) {
384: PCView(data->levels[0]->pc, viewer);
385: }
386: for (i = 1; i < data->N; ++i) {
387: if (data->levels[i]->ksp) color = 1;
388: else color = 0;
389: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
390: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
391: PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm);
392: PetscSubcommSetNumber(subcomm, PetscMin(size, 2));
393: PetscSubcommSetTypeGeneral(subcomm, color, rank);
394: PetscViewerASCIIPushTab(viewer);
395: PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
396: if (color == 1) {
397: KSPView(data->levels[i]->ksp, subviewer);
398: if (data->levels[i]->pc) {
399: PCView(data->levels[i]->pc, subviewer);
400: }
401: PetscViewerFlush(subviewer);
402: }
403: PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
404: PetscViewerASCIIPopTab(viewer);
405: PetscSubcommDestroy(&subcomm);
406: PetscViewerFlush(viewer);
407: }
408: }
409: }
410: return(0);
411: }
413: static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec)
414: {
415: PC_HPDDM *data = (PC_HPDDM*)pc->data;
416: PetscBool flg;
417: Mat A;
421: if (ksp) {
422: PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg);
423: if (flg && !data->normal) {
424: KSPGetOperators(ksp, &A, NULL);
425: MatCreateVecs(A, NULL, &data->normal); /* temporary Vec used in PCHPDDMShellApply() for coarse grid corrections */
426: }
427: }
428: return(0);
429: }
431: static PetscErrorCode PCHPDDMShellSetUp(PC pc)
432: {
433: PC_HPDDM_Level *ctx;
434: Mat A, P;
435: Vec x;
436: const char *pcpre;
440: PCShellGetContext(pc, (void**)&ctx);
441: KSPGetOptionsPrefix(ctx->ksp, &pcpre);
442: KSPGetOperators(ctx->ksp, &A, &P);
443: /* smoother */
444: PCSetOptionsPrefix(ctx->pc, pcpre);
445: PCSetOperators(ctx->pc, A, P);
446: if (!ctx->v[0]) {
447: VecDuplicateVecs(ctx->D, 1, &ctx->v[0]);
448: if (!std::is_same<PetscScalar, PetscReal>::value) {
449: VecDestroy(&ctx->D);
450: }
451: MatCreateVecs(A, &x, NULL);
452: VecDuplicateVecs(x, 2, &ctx->v[1]);
453: VecDestroy(&x);
454: }
455: std::fill_n(ctx->V, 3, nullptr);
456: return(0);
457: }
459: template<class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type* = nullptr>
460: PETSC_STATIC_INLINE PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y)
461: {
462: PC_HPDDM_Level *ctx;
463: PetscScalar *out;
467: PCShellGetContext(pc, (void**)&ctx);
468: /* going from PETSc to HPDDM numbering */
469: VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
470: VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
471: VecGetArrayWrite(ctx->v[0][0], &out);
472: ctx->P->deflation<false>(NULL, out, 1); /* y = Q x */
473: VecRestoreArrayWrite(ctx->v[0][0], &out);
474: /* going from HPDDM to PETSc numbering */
475: VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
476: VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
477: return(0);
478: }
480: template<class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type* = nullptr>
481: PETSC_STATIC_INLINE PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y)
482: {
483: PC_HPDDM_Level *ctx;
484: Vec vX, vY, vC;
485: PetscScalar *out;
486: PetscInt i, N;
490: PCShellGetContext(pc, (void**)&ctx);
491: MatGetSize(X, NULL, &N);
492: /* going from PETSc to HPDDM numbering */
493: for (i = 0; i < N; ++i) {
494: MatDenseGetColumnVecRead(X, i, &vX);
495: MatDenseGetColumnVecWrite(ctx->V[0], i, &vC);
496: VecScatterBegin(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD);
497: VecScatterEnd(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD);
498: MatDenseRestoreColumnVecWrite(ctx->V[0], i, &vC);
499: MatDenseRestoreColumnVecRead(X, i, &vX);
500: }
501: MatDenseGetArrayWrite(ctx->V[0], &out);
502: ctx->P->deflation<false>(NULL, out, N); /* Y = Q X */
503: MatDenseRestoreArrayWrite(ctx->V[0], &out);
504: /* going from HPDDM to PETSc numbering */
505: for (i = 0; i < N; ++i) {
506: MatDenseGetColumnVecRead(ctx->V[0], i, &vC);
507: MatDenseGetColumnVecWrite(Y, i, &vY);
508: VecScatterBegin(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE);
509: VecScatterEnd(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE);
510: MatDenseRestoreColumnVecWrite(Y, i, &vY);
511: MatDenseRestoreColumnVecRead(ctx->V[0], i, &vC);
512: }
513: return(0);
514: }
516: /*@C
517: 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.
519: .vb
520: (1) y = Pmat^-1 x + Q x,
521: (2) y = Pmat^-1 (I - Amat Q) x + Q x (default),
522: (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x.
523: .ve
525: Input Parameters:
526: + pc - preconditioner context
527: - x - input vector
529: Output Parameter:
530: . y - output vector
532: Application Interface Routine: PCApply()
534: Notes:
535: 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.
536: 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.
537: (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.
539: Level: advanced
541: .seealso: PCHPDDM, PCHPDDMCoarseCorrectionType
542: @*/
543: static PetscErrorCode PCHPDDMShellApply(PC pc, Vec x, Vec y)
544: {
545: PC_HPDDM_Level *ctx;
546: Mat A;
550: PCShellGetContext(pc, (void**)&ctx);
551: if (ctx->P) {
552: KSPGetOperators(ctx->ksp, &A, NULL);
553: PCHPDDMDeflate_Private(pc, x, y); /* y = Q x */
554: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
555: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) {
556: MatMult(A, y, ctx->v[1][0]); /* y = A Q x */
557: } else { /* KSPLSQR and finest level */
558: MatMult(A, y, ctx->parent->normal); /* y = A Q x */
559: MatMultTranspose(A, ctx->parent->normal, ctx->v[1][0]); /* y = A^T A Q x */
560: }
561: VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x); /* y = (I - A Q) x */
562: PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0]); /* y = M^-1 (I - A Q) x */
563: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
564: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) {
565: MatMultTranspose(A, ctx->v[1][0], ctx->v[1][1]); /* z = A^T y */
566: } else {
567: MatMult(A, ctx->v[1][0], ctx->parent->normal);
568: MatMultTranspose(A, ctx->parent->normal, ctx->v[1][1]); /* z = A^T A y */
569: }
570: PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1]);
571: VecAXPBYPCZ(y, -1.0, 1.0, 1.0, ctx->v[1][1], ctx->v[1][0]); /* y = (I - Q A^T) y + Q x */
572: } else {
573: VecAXPY(y, 1.0, ctx->v[1][0]); /* y = Q M^-1 (I - A Q) x + Q x */
574: }
575: } else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE) {
576: PCApply(ctx->pc, x, ctx->v[1][0]);
577: VecAXPY(y, 1.0, ctx->v[1][0]); /* y = M^-1 x + Q x */
578: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
579: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
580: return(0);
581: }
583: /*@C
584: PCHPDDMShellMatApply - Variant of PCHPDDMShellApply() for blocks of vectors
586: Input Parameters:
587: + pc - preconditioner context
588: - X - block of input vectors
590: Output Parameter:
591: . Y - block of output vectors
593: Application Interface Routine: PCApply()
595: Level: advanced
597: .seealso: PCHPDDM, PCHPDDMShellMatApply(), PCHPDDMCoarseCorrectionType
598: @*/
599: static PetscErrorCode PCHPDDMShellMatApply(PC pc, Mat X, Mat Y)
600: {
601: PC_HPDDM_Level *ctx;
602: Mat A;
603: PetscScalar *array;
604: PetscInt m, M, N, prev = 0;
608: PCShellGetContext(pc, (void**)&ctx);
609: if (ctx->P) {
610: MatGetSize(X, NULL, &N);
611: if (ctx->V[0]) {
612: MatGetSize(ctx->V[0], NULL, &prev);
613: }
614: KSPGetOperators(ctx->ksp, &A, NULL);
615: if (N != prev) {
616: MatDestroy(ctx->V);
617: MatDestroy(ctx->V + 1);
618: MatDestroy(ctx->V + 2);
619: VecGetLocalSize(ctx->v[0][0], &m);
620: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, ctx->V);
621: MatGetLocalSize(X, &m, NULL);
622: MatGetSize(X, &M, NULL);
623: if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) {
624: MatDenseGetArrayWrite(ctx->V[0], &array);
625: } else array = NULL;
626: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, array, ctx->V + 1);
627: if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) {
628: MatDenseRestoreArrayWrite(ctx->V[0], &array);
629: }
630: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, NULL, ctx->V + 2);
631: MatAssemblyBegin(ctx->V[1], MAT_FINAL_ASSEMBLY);
632: MatAssemblyEnd(ctx->V[1], MAT_FINAL_ASSEMBLY);
633: MatAssemblyBegin(ctx->V[2], MAT_FINAL_ASSEMBLY);
634: MatAssemblyEnd(ctx->V[2], MAT_FINAL_ASSEMBLY);
635: MatProductCreateWithMat(A, Y, NULL, ctx->V[1]);
636: MatProductSetType(ctx->V[1], MATPRODUCT_AB);
637: MatProductSetFromOptions(ctx->V[1]);
638: MatProductSymbolic(ctx->V[1]);
639: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
640: MatProductCreateWithMat(A, ctx->V[1], NULL, ctx->V[2]);
641: MatProductSetType(ctx->V[2], MATPRODUCT_AtB);
642: MatProductSetFromOptions(ctx->V[2]);
643: MatProductSymbolic(ctx->V[2]);
644: }
645: ctx->P->start(N);
646: } else {
647: MatProductReplaceMats(NULL, Y, NULL, ctx->V[1]);
648: }
649: PCHPDDMDeflate_Private(pc, X, Y);
650: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
651: MatProductNumeric(ctx->V[1]);
652: MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN);
653: MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN);
654: PCMatApply(ctx->pc, ctx->V[2], ctx->V[1]);
655: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
656: MatProductNumeric(ctx->V[2]);
657: PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2]);
658: MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN);
659: }
660: MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN);
661: } else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE) {
662: PCMatApply(ctx->pc, X, ctx->V[1]);
663: MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN);
664: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
665: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
666: return(0);
667: }
669: static PetscErrorCode PCHPDDMShellDestroy(PC pc)
670: {
671: PC_HPDDM_Level *ctx;
675: PCShellGetContext(pc, (void**)&ctx);
676: HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE);
677: VecDestroyVecs(1, &ctx->v[0]);
678: VecDestroyVecs(2, &ctx->v[1]);
679: MatDestroy(ctx->V);
680: MatDestroy(ctx->V + 1);
681: MatDestroy(ctx->V + 2);
682: VecDestroy(&ctx->D);
683: VecScatterDestroy(&ctx->scatter);
684: PCDestroy(&ctx->pc);
685: return(0);
686: }
688: static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short& mu)
689: {
690: Mat B, X;
691: PetscInt n, N, j = 0;
695: KSPGetOperators(ctx->ksp, &B, NULL);
696: MatGetLocalSize(B, &n, NULL);
697: MatGetSize(B, &N, NULL);
698: if (ctx->parent->log_separate) {
699: j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx));
700: PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, 0, 0, 0);
701: }
702: if (mu == 1) {
703: if (!ctx->ksp->vec_rhs) {
704: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, NULL, &ctx->ksp->vec_rhs);
705: VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol);
706: }
707: VecPlaceArray(ctx->ksp->vec_rhs, rhs);
708: KSPSolve(ctx->ksp, NULL, NULL);
709: VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs);
710: VecResetArray(ctx->ksp->vec_rhs);
711: } else {
712: MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B);
713: MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, NULL, &X);
714: KSPMatSolve(ctx->ksp, B, X);
715: MatCopy(X, B, SAME_NONZERO_PATTERN);
716: MatDestroy(&X);
717: MatDestroy(&B);
718: }
719: if (ctx->parent->log_separate) {
720: PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, 0, 0, 0);
721: }
722: return(0);
723: }
725: static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc)
726: {
727: PC_HPDDM *data = (PC_HPDDM*)pc->data;
731: if (data->setup) {
732: Mat P;
733: Vec x, xt = NULL;
734: PetscReal t = 0.0, s = 0.0;
736: PCGetOperators(pc, NULL, &P);
737: PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject*)&x);
738: PetscStackPush("PCHPDDM Neumann callback");
739: (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx);
740: PetscStackPop;
741: }
742: return(0);
743: }
745: static PetscErrorCode PCSetUp_HPDDM(PC pc)
746: {
747: PC_HPDDM *data = (PC_HPDDM*)pc->data;
748: PC inner;
749: KSP *ksp;
750: Mat *sub, A, P, N, C = NULL, uaux = NULL, weighted, subA[2];
751: Vec xin, v;
752: std::vector<Vec> initial;
753: IS is[1], loc, uis = data->is;
754: ISLocalToGlobalMapping l2g;
755: char prefix[256];
756: const char *pcpre;
757: const PetscScalar *const *ev;
758: PetscInt n, requested = data->N, reused = 0;
759: PetscBool subdomains = PETSC_FALSE, flag = PETSC_FALSE, ismatis, swap = PETSC_FALSE;
760: DM dm;
761: PetscErrorCode ierr;
764: if (!data->levels || !data->levels[0]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated");
765: PCGetOptionsPrefix(pc, &pcpre);
766: PCGetOperators(pc, &A, &P);
767: if (!data->levels[0]->ksp) {
768: KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp);
769: PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse");
770: KSPSetOptionsPrefix(data->levels[0]->ksp, prefix);
771: KSPSetType(data->levels[0]->ksp, KSPPREONLY);
772: } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) {
773: /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */
774: /* then just propagate the appropriate flag to the coarser levels */
775: for (n = 0; n < PETSC_HPDDM_MAXLEVELS && data->levels[n]; ++n) {
776: /* the following KSP and PC may be NULL for some processes, hence the check */
777: if (data->levels[n]->ksp) {
778: KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE);
779: }
780: if (data->levels[n]->pc) {
781: PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE);
782: }
783: }
784: /* early bail out because there is nothing to do */
785: return(0);
786: } else {
787: /* reset coarser levels */
788: for (n = 1; n < PETSC_HPDDM_MAXLEVELS && data->levels[n]; ++n) {
789: if (data->levels[n]->ksp && data->levels[n]->ksp->pc && data->levels[n]->ksp->pc->setupcalled == 1 && data->levels[n]->ksp->pc->reusepreconditioner && n < data->N) {
790: reused = data->N - n;
791: break;
792: }
793: KSPDestroy(&data->levels[n]->ksp);
794: PCDestroy(&data->levels[n]->pc);
795: }
796: /* check if some coarser levels are being reused */
797: MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc));
798: const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0;
800: if (addr != &HPDDM::i__0 && reused != data->N - 1) {
801: /* reuse previously computed eigenvectors */
802: ev = data->levels[0]->P->getVectors();
803: if (ev) {
804: initial.reserve(*addr);
805: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin);
806: for (n = 0; n < *addr; ++n) {
807: VecDuplicate(xin, &v);
808: VecPlaceArray(xin, ev[n]);
809: VecCopy(xin, v);
810: initial.emplace_back(v);
811: VecResetArray(xin);
812: }
813: VecDestroy(&xin);
814: }
815: }
816: }
817: data->N -= reused;
818: KSPSetOperators(data->levels[0]->ksp, A, P);
820: PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis);
821: if (!data->is && !ismatis) {
822: PetscErrorCode (*create)(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void**) = NULL;
823: PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*) = NULL;
824: void *uctx = NULL;
826: /* first see if we can get the data from the DM */
827: MatGetDM(P, &dm);
828: if (!dm) {
829: MatGetDM(A, &dm);
830: }
831: if (!dm) {
832: PCGetDM(pc, &dm);
833: }
834: if (dm) { /* this is the hook for DMPLEX and DMDA for which the auxiliary Mat is the local Neumann matrix */
835: PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create);
836: if (create) {
837: (*create)(dm, &uis, &uaux, &usetup, &uctx);
838: data->Neumann = PETSC_TRUE;
839: }
840: }
841: if (!create) {
842: if (!uis) {
843: PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject*)&uis);
844: PetscObjectReference((PetscObject)uis);
845: }
846: if (!uaux) {
847: PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject*)&uaux);
848: PetscObjectReference((PetscObject)uaux);
849: }
850: /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */
851: if (!uis) {
852: PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject*)&uis);
853: PetscObjectReference((PetscObject)uis);
854: }
855: if (!uaux) {
856: PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject*)&uaux);
857: PetscObjectReference((PetscObject)uaux);
858: }
859: }
860: PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx);
861: MatDestroy(&uaux);
862: ISDestroy(&uis);
863: }
865: if (!ismatis) {
866: PCHPDDMSetUpNeumannOverlap_Private(pc);
867: }
869: if (data->is || (ismatis && data->N > 1)) {
870: if (ismatis) {
871: std::initializer_list<std::string> list = { MATSEQBAIJ, MATSEQSBAIJ };
872: MatISGetLocalMat(P, &N);
873: std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name);
874: MatISRestoreLocalMat(P, &N);
875: switch (std::distance(list.begin(), it)) {
876: case 0:
877: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
878: break;
879: case 1:
880: /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */
881: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
882: MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE);
883: break;
884: default:
885: MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C);
886: }
887: MatGetLocalToGlobalMapping(P, &l2g, NULL);
888: PetscObjectReference((PetscObject)P);
889: KSPSetOperators(data->levels[0]->ksp, A, C);
890: std::swap(C, P);
891: ISLocalToGlobalMappingGetSize(l2g, &n);
892: ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc);
893: ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0]);
894: ISDestroy(&loc);
895: /* the auxiliary Mat is _not_ the local Neumann matrix */
896: /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */
897: data->Neumann = PETSC_FALSE;
898: } else {
899: is[0] = data->is;
900: PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_define_subdomains", &subdomains, NULL);
901: PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_has_neumann", &data->Neumann, NULL);
902: ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc);
903: }
904: if (data->N > 1 && (data->aux || ismatis)) {
905: if (!loadedSym) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level");
906: MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
907: if (ismatis) {
908: /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */
909: MatIncreaseOverlap(P, 1, is, 1);
910: ISDestroy(&data->is);
911: data->is = is[0];
912: } else {
913: if (PetscDefined(USE_DEBUG)) {
914: PetscBool equal;
915: IS intersect;
917: ISIntersect(data->is, loc, &intersect);
918: ISEqualUnsorted(loc, intersect, &equal);
919: ISDestroy(&intersect);
920: if (!equal) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS of the auxiliary Mat does not include all local rows of A");
921: }
922: if (!data->Neumann) {
923: PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flag);
924: if (flag) {
925: /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */
926: /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */
927: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux);
928: flag = PETSC_FALSE;
929: }
930: }
931: }
932: if (!uaux) {
933: if (data->Neumann) sub = &data->aux;
934: else {
935: MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub);
936: }
937: } else {
938: MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub);
939: MatDestroy(&uaux);
940: MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub);
941: }
942: /* Vec holding the partition of unity */
943: if (!data->levels[0]->D) {
944: ISGetLocalSize(data->is, &n);
945: VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D);
946: }
947: /* it is possible to share the PC only given specific conditions, otherwise there is not warranty that the matrices have the same nonzero pattern */
948: if (!ismatis && sub == &data->aux && !data->B && subdomains && data->share) {
949: IS perm;
950: PetscInt size;
951: const PetscInt *ptr;
952: PetscInt *concatenate;
953: ISGetLocalSize(*is, &size);
954: ISGetIndices(*is, &ptr);
955: std::map<PetscInt, PetscInt> order;
956: /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */
957: for (n = 0; n < size; ++n)
958: order.insert(std::make_pair(ptr[n], n));
959: ISRestoreIndices(*is, &ptr);
960: PetscMalloc1(size, &concatenate);
961: for (const std::pair<const PetscInt, PetscInt>& i : order)
962: *concatenate++ = i.second;
963: concatenate -= size;
964: ISCreateGeneral(PETSC_COMM_SELF, size, concatenate, PETSC_OWN_POINTER, &perm);
965: ISSetPermutation(perm);
966: /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */
967: MatPermute(data->aux, perm, perm, &C);
968: ISDestroy(&perm);
969: PetscMalloc1(size, &concatenate);
970: for (const std::pair<const PetscInt, PetscInt>& i : order)
971: *concatenate++ = i.first;
972: concatenate -= size;
973: /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */
974: ISCreateGeneral(PetscObjectComm((PetscObject)data->is), size, concatenate, PETSC_OWN_POINTER, &uis);
975: if (!data->levels[0]->pc) {
976: PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "");
977: PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc);
978: PCSetOptionsPrefix(data->levels[0]->pc, prefix);
979: PCSetOperators(data->levels[0]->pc, A, P);
980: }
981: PCSetType(data->levels[0]->pc, PCASM);
982: if (!data->levels[0]->pc->setupcalled) {
983: PCASMSetLocalSubdomains(data->levels[0]->pc, 1, is, &loc);
984: }
985: PCSetFromOptions(data->levels[0]->pc);
986: PCSetUp(data->levels[0]->pc);
987: size = -1;
988: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt*, PetscInt*, KSP**), (data->levels[0]->pc, &size, NULL, &ksp));
989: if (size != 1) {
990: PCDestroy(&data->levels[0]->pc);
991: MatDestroy(&C);
992: ISDestroy(&uis);
993: data->share = PETSC_FALSE;
994: if (size == -1) {
995: PetscInfo(pc, "Cannot share PC between ST and subdomain solver since PCASMGetSubKSP() not found in fine-level PC\n");
996: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %D != 1", size);
997: } else {
998: Mat D;
999: const char *matpre;
1000: KSPGetOperators(ksp[0], subA, subA + 1);
1001: MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D);
1002: MatGetOptionsPrefix(subA[1], &matpre);
1003: MatSetOptionsPrefix(D, matpre);
1004: MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN);
1005: MatPropagateSymmetryOptions(C, D);
1006: MatDestroy(&C);
1007: C = D;
1008: /* swap pointers so that variables stay consistent throughout PCSetUp() */
1009: std::swap(C, data->aux);
1010: std::swap(uis, data->is);
1011: swap = PETSC_TRUE;
1012: }
1013: } else if (data->share) {
1014: PetscInfo(pc, "Cannot share PC between ST and subdomain solver\n");
1015: data->share = PETSC_FALSE;
1016: }
1017: if (!data->levels[0]->scatter) {
1018: MatCreateVecs(P, &xin, NULL);
1019: if (ismatis) {
1020: MatDestroy(&P);
1021: }
1022: VecScatterCreate(xin, data->is, data->levels[0]->D, NULL, &data->levels[0]->scatter);
1023: VecDestroy(&xin);
1024: }
1025: if (data->levels[0]->P) {
1026: /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */
1027: HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE);
1028: }
1029: if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>();
1030: if (data->log_separate) {
1031: PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, 0, 0, 0);
1032: } else {
1033: PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, 0, 0, 0);
1034: }
1035: /* HPDDM internal data structure */
1036: data->levels[0]->P->structure(loc, data->is, sub[0], ismatis ? C : data->aux, data->levels);
1037: if (!data->log_separate) {
1038: PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, 0, 0, 0);
1039: }
1040: /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */
1041: if (!data->B) {
1042: MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted);
1043: MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D);
1044: /* neither MatDuplicate() nor MatDiagonaleScale() handles the symmetry options, so propagate the options explicitly */
1045: /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */
1046: MatPropagateSymmetryOptions(sub[0], weighted);
1047: } else weighted = data->B;
1048: /* SLEPc is used inside the loaded symbol */
1049: (*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : data->aux, weighted, data->B, initial, data->levels);
1050: if (data->share) {
1051: Mat st[2];
1052: KSPGetOperators(ksp[0], st, st + 1);
1053: MatCopy(subA[0], st[0], SAME_NONZERO_PATTERN);
1054: if (subA[1] != subA[0] || st[1] != st[0]) {
1055: MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN);
1056: }
1057: }
1058: if (data->log_separate) {
1059: PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, 0, 0, 0);
1060: }
1061: if (ismatis) {
1062: MatISGetLocalMat(C, &N);
1063: } else N = data->aux;
1064: P = sub[0];
1065: /* going through the grid hierarchy */
1066: for (n = 1; n < data->N; ++n) {
1067: if (data->log_separate) {
1068: PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, 0, 0, 0);
1069: }
1070: /* method composed in the loaded symbol since there, SLEPc is used as well */
1071: PetscUseMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat*, Mat*, PetscInt, PetscInt* const, PC_HPDDM_Level** const), (&P, &N, n, &data->N, data->levels));
1072: if (data->log_separate) {
1073: PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, 0, 0, 0);
1074: }
1075: }
1076: /* reset to NULL to avoid any faulty use */
1077: PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", NULL);
1078: if (ismatis) {
1079: /* matching PetscObjectReference() above */
1080: PetscObjectDereference((PetscObject)C);
1081: }
1082: for (n = 0; n < data->N - 1; ++n)
1083: if (data->levels[n]->P) {
1084: /* HPDDM internal work buffers */
1085: data->levels[n]->P->setBuffer();
1086: data->levels[n]->P->super::start();
1087: }
1088: if (!data->Neumann) {
1089: MatDestroySubMatrices(1, &sub);
1090: }
1091: if (ismatis) data->is = NULL;
1092: for (n = 0; n < data->N - 1 + (reused > 0); ++n) {
1093: if (data->levels[n]->P) {
1094: PC spc;
1096: /* force the PC to be PCSHELL to do the coarse grid corrections */
1097: KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE);
1098: KSPGetPC(data->levels[n]->ksp, &spc);
1099: PCSetType(spc, PCSHELL);
1100: PCShellSetContext(spc, data->levels[n]);
1101: PCShellSetSetUp(spc, PCHPDDMShellSetUp);
1102: PCShellSetApply(spc, PCHPDDMShellApply);
1103: PCShellSetMatApply(spc, PCHPDDMShellMatApply);
1104: PCShellSetDestroy(spc, PCHPDDMShellDestroy);
1105: if (!data->levels[n]->pc) {
1106: PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc);
1107: }
1108: if (n < reused) {
1109: PCSetReusePreconditioner(spc, PETSC_TRUE);
1110: PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE);
1111: }
1112: PCSetUp(spc);
1113: }
1114: }
1115: } else flag = reused ? PETSC_FALSE : PETSC_TRUE;
1116: if (!ismatis && subdomains) {
1117: if (flag) {
1118: KSPGetPC(data->levels[0]->ksp, &inner);
1119: } else inner = data->levels[0]->pc;
1120: if (inner) {
1121: PCSetType(inner, PCASM);
1122: if (!inner->setupcalled) {
1123: PCASMSetLocalSubdomains(inner, 1, is, &loc);
1124: }
1125: }
1126: }
1127: ISDestroy(&loc);
1128: } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */
1129: if (requested != data->N + reused) {
1130: PetscInfo5(pc, "%D levels requested, only %D built + %D reused. Options for level(s) > %D, including -%spc_hpddm_coarse_ will not be taken into account\n", requested, data->N, reused, data->N, pcpre ? pcpre : "");
1131: PetscInfo2(pc, "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);
1132: /* cannot use PCHPDDMShellDestroy() because PCSHELL not set for unassembled levels */
1133: for (n = data->N - 1; n < requested - 1; ++n) {
1134: if (data->levels[n]->P) {
1135: HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE);
1136: VecDestroyVecs(1, &data->levels[n]->v[0]);
1137: VecDestroyVecs(2, &data->levels[n]->v[1]);
1138: MatDestroy(data->levels[n]->V);
1139: MatDestroy(data->levels[n]->V + 1);
1140: MatDestroy(data->levels[n]->V + 2);
1141: VecDestroy(&data->levels[n]->D);
1142: VecScatterDestroy(&data->levels[n]->scatter);
1143: }
1144: }
1145: if (reused) {
1146: for (n = reused; n < PETSC_HPDDM_MAXLEVELS && data->levels[n]; ++n) {
1147: KSPDestroy(&data->levels[n]->ksp);
1148: PCDestroy(&data->levels[n]->pc);
1149: }
1150: }
1151: if (PetscDefined(USE_DEBUG)) SETERRQ7(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "%D levels requested, only %D built + %D reused. 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, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N);
1152: }
1153: /* these solvers are created after PCSetFromOptions() is called */
1154: if (pc->setfromoptionscalled) {
1155: for (n = 0; n < data->N; ++n) {
1156: if (data->levels[n]->ksp) {
1157: KSPSetFromOptions(data->levels[n]->ksp);
1158: }
1159: if (data->levels[n]->pc) {
1160: PCSetFromOptions(data->levels[n]->pc);
1161: }
1162: }
1163: pc->setfromoptionscalled = 0;
1164: }
1165: data->N += reused;
1166: if (data->share && swap) {
1167: /* swap back pointers so that variables follow the user-provided numbering */
1168: std::swap(C, data->aux);
1169: std::swap(uis, data->is);
1170: MatDestroy(&C);
1171: ISDestroy(&uis);
1172: }
1173: return(0);
1174: }
1176: /*@
1177: PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type.
1179: Input Parameters:
1180: + pc - preconditioner context
1181: - type - PC_HPDDM_COARSE_CORRECTION_DEFLATED, PC_HPDDM_COARSE_CORRECTION_ADDITIVE, or PC_HPDDM_COARSE_CORRECTION_BALANCED
1183: Options Database Key:
1184: . -pc_hpddm_coarse_correction <deflated, additive, balanced> - type of coarse correction to apply
1186: Level: intermediate
1188: .seealso: PCHPDDMGetCoarseCorrectionType(), PCHPDDM, PCHPDDMCoarseCorrectionType
1189: @*/
1190: PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type)
1191: {
1196: PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type));
1197: return(0);
1198: }
1200: /*@
1201: PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type.
1203: Input Parameter:
1204: . pc - preconditioner context
1206: Output Parameter:
1207: . type - PC_HPDDM_COARSE_CORRECTION_DEFLATED, PC_HPDDM_COARSE_CORRECTION_ADDITIVE, or PC_HPDDM_COARSE_CORRECTION_BALANCED
1209: Level: intermediate
1211: .seealso: PCHPDDMSetCoarseCorrectionType(), PCHPDDM, PCHPDDMCoarseCorrectionType
1212: @*/
1213: PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type)
1214: {
1219: PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType*), (pc, type));
1220: return(0);
1221: }
1223: static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type)
1224: {
1225: PC_HPDDM *data = (PC_HPDDM*)pc->data;
1228: if (type < 0 || type > 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PCHPDDMCoarseCorrectionType %d", type);
1229: data->correction = type;
1230: return(0);
1231: }
1233: static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type)
1234: {
1235: PC_HPDDM *data = (PC_HPDDM*)pc->data;
1238: *type = data->correction;
1239: return(0);
1240: }
1242: /*@
1243: PCHPDDMGetSTShareSubKSP - Gets whether the KSP in SLEPc ST and the fine-level subdomain solver is shared.
1245: Input Parameter:
1246: . pc - preconditioner context
1248: Output Parameter:
1249: . share - whether the KSP is shared or not
1251: Notes:
1252: This is not the same as PCGetReusePreconditioner(). The return value is unlikely to be true, but when it is, a symbolic factorization can be skipped
1253: when using a subdomain PCType such as PCLU or PCCHOLESKY.
1255: Level: advanced
1257: @*/
1258: PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share)
1259: {
1264: PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool*), (pc, share));
1265: return(0);
1266: }
1268: static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share)
1269: {
1270: PC_HPDDM *data = (PC_HPDDM*)pc->data;
1273: *share = data->share;
1274: return(0);
1275: }
1277: PetscErrorCode HPDDMLoadDL_Private(PetscBool *found)
1278: {
1279: char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN];
1283: PetscStrcpy(dir, "${PETSC_LIB_DIR}");
1284: PetscOptionsGetString(NULL, NULL, "-hpddm_dir", dir, sizeof(dir), NULL);
1285: PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir);
1286: PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found);
1287: if (*found) {
1288: PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib);
1289: #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure since */
1290: } else { /* slepcconf.h is not yet built (and thus can't be included) */
1291: PetscStrcpy(dir, HPDDM_STR(SLEPC_LIB_DIR));
1292: PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir);
1293: PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found);
1294: if (*found) {
1295: PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib);
1296: #endif
1297: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib);
1298: #if defined(SLEPC_LIB_DIR)
1299: }
1300: #endif
1301: return(0);
1302: }
1304: /*MC
1305: PCHPDDM - Interface with the HPDDM library.
1307: 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. A chronological bibliography of relevant publications linked with PC available in HPDDM through PCHPDDM may be found below. The interface is explained in details in [2021].
1309: 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).
1311: Options Database Keys:
1312: + -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)
1313: . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the PC that the local Neumann matrix is supplied in PCHPDDMSetAuxiliaryMat()
1314: - -pc_hpddm_coarse_correction <type, default=deflated> - determines the PCHPDDMCoarseCorrectionType when calling PCApply
1316: Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set with
1317: .vb
1318: -pc_hpddm_levels_%d_pc_
1319: -pc_hpddm_levels_%d_ksp_
1320: -pc_hpddm_levels_%d_eps_
1321: -pc_hpddm_levels_%d_p
1322: -pc_hpddm_levels_%d_mat_type_
1323: -pc_hpddm_coarse_
1324: -pc_hpddm_coarse_p
1325: -pc_hpddm_coarse_mat_type_
1326: .ve
1327: 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).
1329: 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.
1331: 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.
1333: References:
1334: + 2011 - A robust two-level domain decomposition preconditioner for systems of PDEs. Spillane, Dolean, Hauret, Nataf, Pechstein, and Scheichl. Comptes Rendus Mathematique.
1335: . 2013 - Scalable Domain Decomposition Preconditioners For Heterogeneous Elliptic Problems. Jolivet, Hecht, Nataf, and Prud'homme. SC13.
1336: . 2015 - An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation. Dolean, Jolivet, and Nataf. SIAM.
1337: . 2019 - A Multilevel Schwarz Preconditioner Based on a Hierarchy of Robust Coarse Spaces. Al Daas, Grigori, Jolivet, and Tournier.
1338: - 2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini. Computer & Mathematics with Applications.
1340: Level: intermediate
1342: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCHPDDMSetAuxiliaryMat(), MATIS, PCBDDC, PCDEFLATION, PCTELESCOPE
1343: M*/
1344: PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc)
1345: {
1346: PC_HPDDM *data;
1347: PetscBool found;
1351: if (!loadedSym) {
1352: HPDDMLoadDL_Private(&found);
1353: if (found) {
1354: PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, NULL, "PCHPDDM_Internal", (void**)&loadedSym);
1355: }
1356: }
1357: if (!loadedSym) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc");
1358: PetscNewLog(pc, &data);
1359: pc->data = data;
1360: pc->ops->reset = PCReset_HPDDM;
1361: pc->ops->destroy = PCDestroy_HPDDM;
1362: pc->ops->setfromoptions = PCSetFromOptions_HPDDM;
1363: pc->ops->setup = PCSetUp_HPDDM;
1364: pc->ops->apply = PCApply_HPDDM;
1365: pc->ops->matapply = PCMatApply_HPDDM;
1366: pc->ops->view = PCView_HPDDM;
1367: pc->ops->applytranspose = 0;
1368: pc->ops->applysymmetricleft = 0;
1369: pc->ops->applysymmetricright = 0;
1370: pc->ops->presolve = PCPreSolve_HPDDM;
1371: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM);
1372: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM);
1373: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM);
1374: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM);
1375: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM);
1376: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM);
1377: return(0);
1378: }
1380: /*@C
1381: PCHPDDMInitializePackage - This function initializes everything in the PCHPDDM package. It is called from PCInitializePackage().
1383: Level: intermediate
1385: .seealso: PetscInitialize()
1386: @*/
1387: PetscErrorCode PCHPDDMInitializePackage(void)
1388: {
1389: char ename[32];
1390: PetscInt i;
1394: if (PCHPDDMPackageInitialized) return(0);
1395: PCHPDDMPackageInitialized = PETSC_TRUE;
1396: PetscRegisterFinalize(PCHPDDMFinalizePackage);
1397: /* general events registered once during package initialization */
1398: /* some of these events are not triggered in libpetsc, */
1399: /* but rather directly in libhpddm_petsc, */
1400: /* which is in charge of performing the following operations */
1402: /* domain decomposition structure from Pmat sparsity pattern */
1403: PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc);
1404: /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */
1405: PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP);
1406: /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */
1407: PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP);
1408: /* next level construction using PtAP and PtBP (not triggered in libpetsc) */
1409: PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next);
1410: static_assert(PETSC_HPDDM_MAXLEVELS <= 9, "PETSC_HPDDM_MAXLEVELS value is too high");
1411: for (i = 1; i < PETSC_HPDDM_MAXLEVELS; ++i) {
1412: PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1d", i);
1413: /* events during a PCSetUp() at level #i _except_ the assembly */
1414: /* of the Galerkin operator of the coarser level #(i + 1) */
1415: PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1]);
1416: PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1d", i);
1417: /* events during a PCApply() at level #i _except_ */
1418: /* the KSPSolve() of the coarser level #(i + 1) */
1419: PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1]);
1420: }
1421: return(0);
1422: }
1424: /*@C
1425: PCHPDDMFinalizePackage - This function frees everything from the PCHPDDM package. It is called from PetscFinalize().
1427: Level: intermediate
1429: .seealso: PetscFinalize()
1430: @*/
1431: PetscErrorCode PCHPDDMFinalizePackage(void)
1432: {
1434: PCHPDDMPackageInitialized = PETSC_FALSE;
1435: return(0);
1436: }