Actual source code: pchpddm.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_PCHPDDM_MAXLEVELS];
19: PetscLogEvent PC_HPDDM_Solve[PETSC_PCHPDDM_MAXLEVELS];
21: const char *const PCHPDDMCoarseCorrectionTypes[] = {"DEFLATED", "ADDITIVE", "BALANCED", "PCHPDDMCoarseCorrectionType", "PC_HPDDM_COARSE_CORRECTION_", NULL};
23: static PetscErrorCode PCReset_HPDDM(PC pc)
24: {
25: PC_HPDDM *data = (PC_HPDDM *)pc->data;
26: PetscInt i;
28: if (data->levels) {
29: for (i = 0; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) {
30: KSPDestroy(&data->levels[i]->ksp);
31: PCDestroy(&data->levels[i]->pc);
32: PetscFree(data->levels[i]);
33: }
34: PetscFree(data->levels);
35: }
37: ISDestroy(&data->is);
38: MatDestroy(&data->aux);
39: MatDestroy(&data->B);
40: VecDestroy(&data->normal);
41: data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED;
42: data->Neumann = PETSC_FALSE;
43: data->deflation = PETSC_FALSE;
44: data->setup = NULL;
45: data->setup_ctx = NULL;
46: return 0;
47: }
49: static PetscErrorCode PCDestroy_HPDDM(PC pc)
50: {
51: PC_HPDDM *data = (PC_HPDDM *)pc->data;
53: PCReset_HPDDM(pc);
54: PetscFree(data);
55: PetscObjectChangeTypeName((PetscObject)pc, NULL);
56: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", NULL);
57: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", NULL);
58: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", NULL);
59: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", NULL);
60: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", NULL);
61: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", NULL);
62: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", NULL);
63: return 0;
64: }
66: static inline PetscErrorCode PCHPDDMSetAuxiliaryMat_Private(PC pc, IS is, Mat A, PetscBool deflation)
67: {
68: PC_HPDDM *data = (PC_HPDDM *)pc->data;
70: if (is) {
71: PetscObjectReference((PetscObject)is);
72: if (data->is) { /* new overlap definition resets the PC */
73: PCReset_HPDDM(pc);
74: pc->setfromoptionscalled = 0;
75: pc->setupcalled = 0;
76: }
77: ISDestroy(&data->is);
78: data->is = is;
79: }
80: if (A) {
81: PetscObjectReference((PetscObject)A);
82: MatDestroy(&data->aux);
83: data->aux = A;
84: }
85: data->deflation = deflation;
86: return 0;
87: }
89: static inline PetscErrorCode PCHPDDMSetAuxiliaryMatNormal_Private(PC pc, Mat A, Mat N, Mat *B, const char *pcpre)
90: {
91: PC_HPDDM *data = (PC_HPDDM *)pc->data;
92: Mat *splitting, *sub, aux;
93: IS is, cols[2], rows;
94: PetscReal norm;
95: PetscBool flg;
96: char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */
98: MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, B);
99: ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, cols);
100: ISCreateStride(PETSC_COMM_SELF, A->rmap->N, 0, 1, &rows);
101: MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
102: MatIncreaseOverlap(*B, 1, cols, 1);
103: MatCreateSubMatrices(A, 1, &rows, cols, MAT_INITIAL_MATRIX, &splitting);
104: ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, &is);
105: ISEmbed(*cols, is, PETSC_TRUE, cols + 1);
106: ISDestroy(&is);
107: MatCreateSubMatrices(*splitting, 1, &rows, cols + 1, MAT_INITIAL_MATRIX, &sub);
108: ISDestroy(cols + 1);
109: MatFindZeroRows(*sub, &is);
110: MatDestroySubMatrices(1, &sub);
111: ISDestroy(&rows);
112: MatSetOption(*splitting, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
113: MatZeroRowsIS(*splitting, is, 0.0, NULL, NULL);
114: ISDestroy(&is);
115: PetscOptionsGetString(NULL, pcpre, "-pc_hpddm_levels_1_sub_pc_type", type, sizeof(type), NULL);
116: PetscStrcmp(type, PCQR, &flg);
117: if (!flg) {
118: Mat conjugate = *splitting;
119: if (PetscDefined(USE_COMPLEX)) {
120: MatDuplicate(*splitting, MAT_COPY_VALUES, &conjugate);
121: MatConjugate(conjugate);
122: }
123: MatTransposeMatMult(conjugate, *splitting, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &aux);
124: if (PetscDefined(USE_COMPLEX)) MatDestroy(&conjugate);
125: MatNorm(aux, NORM_FROBENIUS, &norm);
126: MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
127: MatShift(aux, PETSC_SMALL * norm);
128: } else {
129: PetscBool flg;
130: PetscObjectTypeCompare((PetscObject)N, MATNORMAL, &flg);
131: if (flg) MatCreateNormal(*splitting, &aux);
132: else MatCreateNormalHermitian(*splitting, &aux);
133: }
134: MatDestroySubMatrices(1, &splitting);
135: PCHPDDMSetAuxiliaryMat(pc, *cols, aux, NULL, NULL);
136: data->Neumann = PETSC_TRUE;
137: ISDestroy(cols);
138: MatDestroy(&aux);
139: return 0;
140: }
142: static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx)
143: {
144: PC_HPDDM *data = (PC_HPDDM *)pc->data;
146: PCHPDDMSetAuxiliaryMat_Private(pc, is, A, PETSC_FALSE);
147: if (setup) {
148: data->setup = setup;
149: data->setup_ctx = setup_ctx;
150: }
151: return 0;
152: }
154: /*@
155: PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by `PCHPDDM` for the concurrent GenEO problems 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.
157: Input Parameters:
158: + pc - preconditioner context
159: . is - index set of the local auxiliary, e.g., Neumann, matrix
160: . A - auxiliary sequential matrix
161: . setup - function for generating the auxiliary matrix
162: - setup_ctx - context for setup
164: Level: intermediate
166: .seealso: `PCHPDDM`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetRHSMat()`, `MATIS`
167: @*/
168: PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx)
169: {
173: #if defined(PETSC_HAVE_FORTRAN)
174: if (reinterpret_cast<void *>(setup) == reinterpret_cast<void *>(PETSC_NULL_FUNCTION_Fortran)) setup = NULL;
175: if (setup_ctx == PETSC_NULL_INTEGER_Fortran) setup_ctx = NULL;
176: #endif
177: PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode(*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *), (pc, is, A, setup, setup_ctx));
178: return 0;
179: }
181: static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has)
182: {
183: PC_HPDDM *data = (PC_HPDDM *)pc->data;
185: data->Neumann = has;
186: return 0;
187: }
189: /*@
190: PCHPDDMHasNeumannMat - Informs `PCHPDDM` that the `Mat` passed to `PCHPDDMSetAuxiliaryMat()` is the local Neumann matrix.
192: Input Parameters:
193: + pc - preconditioner context
194: - has - Boolean value
196: Level: intermediate
198: Notes:
199: This may be used to bypass a call to `MatCreateSubMatrices()` and to `MatConvert()` for `MATSBAIJ` matrices.
201: 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`.
203: .seealso: `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()`
204: @*/
205: PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has)
206: {
208: PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has));
209: return 0;
210: }
212: static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B)
213: {
214: PC_HPDDM *data = (PC_HPDDM *)pc->data;
216: PetscObjectReference((PetscObject)B);
217: MatDestroy(&data->B);
218: data->B = B;
219: return 0;
220: }
222: /*@
223: PCHPDDMSetRHSMat - Sets the right-hand side matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. Must be used in conjunction 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.
225: Input Parameters:
226: + pc - preconditioner context
227: - B - right-hand side sequential matrix
229: Level: advanced
231: .seealso: `PCHPDDMSetAuxiliaryMat()`, `PCHPDDM`
232: @*/
233: PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B)
234: {
236: if (B) {
238: PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B));
239: }
240: return 0;
241: }
243: static PetscErrorCode PCSetFromOptions_HPDDM(PC pc, PetscOptionItems *PetscOptionsObject)
244: {
245: PC_HPDDM *data = (PC_HPDDM *)pc->data;
246: PC_HPDDM_Level **levels = data->levels;
247: char prefix[256];
248: int i = 1;
249: PetscMPIInt size, previous;
250: PetscInt n;
251: PCHPDDMCoarseCorrectionType type;
252: PetscBool flg = PETSC_TRUE;
254: if (!data->levels) {
255: PetscCalloc1(PETSC_PCHPDDM_MAXLEVELS, &levels);
256: data->levels = levels;
257: }
258: PetscOptionsHeadBegin(PetscOptionsObject, "PCHPDDM options");
259: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
260: previous = size;
261: while (i < PETSC_PCHPDDM_MAXLEVELS) {
262: PetscInt p = 1;
264: if (!data->levels[i - 1]) PetscNew(data->levels + i - 1);
265: data->levels[i - 1]->parent = data;
266: /* if the previous level has a single process, it is not possible to coarsen further */
267: if (previous == 1 || !flg) break;
268: data->levels[i - 1]->nu = 0;
269: data->levels[i - 1]->threshold = -1.0;
270: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
271: PetscOptionsInt(prefix, "Local number of deflation vectors computed by SLEPc", "EPSSetDimensions", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, NULL);
272: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
273: PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, NULL);
274: if (i == 1) {
275: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp");
276: PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "PCHPDDMGetSTShareSubKSP", PETSC_FALSE, &data->share, NULL);
277: }
278: /* if there is no prescribed coarsening, just break out of the loop */
279: if (data->levels[i - 1]->threshold <= 0.0 && data->levels[i - 1]->nu <= 0 && !(data->deflation && i == 1)) break;
280: else {
281: ++i;
282: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i);
283: PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
284: if (!flg) {
285: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i);
286: PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg);
287: }
288: if (flg) {
289: /* if there are coarsening options for the next level, then register it */
290: /* otherwise, don't to avoid having both options levels_N_p and coarse_p */
291: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i);
292: PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "PCHPDDM", p, &p, &flg, 1, PetscMax(1, previous / 2));
293: previous = p;
294: }
295: }
296: }
297: data->N = i;
298: n = 1;
299: if (i > 1) {
300: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p");
301: PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "PCHPDDM", n, &n, NULL, 1, PetscMax(1, previous / 2));
302: #if defined(PETSC_HAVE_MUMPS)
303: PetscSNPrintf(prefix, sizeof(prefix), "pc_hpddm_coarse_");
304: PetscOptionsHasName(NULL, prefix, "-mat_mumps_use_omp_threads", &flg);
305: if (flg) {
306: char type[64]; /* same size as in src/ksp/pc/impls/factor/factimpl.c */
307: if (n == 1) PetscStrcpy(type, MATSOLVERPETSC); /* default solver for a sequential Mat */
308: PetscOptionsGetString(NULL, prefix, "-pc_factor_mat_solver_type", type, sizeof(type), &flg);
309: if (flg) PetscStrcmp(type, MATSOLVERMUMPS, &flg);
311: size = n;
312: n = -1;
313: PetscOptionsGetInt(NULL, prefix, "-mat_mumps_use_omp_threads", &n, NULL);
316: }
317: #endif
318: PetscOptionsEnum("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, (PetscEnum)data->correction, (PetscEnum *)&type, &flg);
319: if (flg) PCHPDDMSetCoarseCorrectionType(pc, type);
320: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann");
321: PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", data->Neumann, &data->Neumann, NULL);
322: data->log_separate = PETSC_FALSE;
323: if (PetscDefined(USE_LOG)) {
324: PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate");
325: PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", NULL, data->log_separate, &data->log_separate, NULL);
326: }
327: }
328: PetscOptionsHeadEnd();
329: while (i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]) PetscFree(data->levels[i++]);
330: return 0;
331: }
333: static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y)
334: {
335: PC_HPDDM *data = (PC_HPDDM *)pc->data;
337: PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
339: if (data->log_separate) PetscLogEventBegin(PC_HPDDM_Solve[0], data->levels[0]->ksp, 0, 0, 0); /* coarser-level events are directly triggered in HPDDM */
340: KSPSolve(data->levels[0]->ksp, x, y);
341: if (data->log_separate) PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, 0, 0, 0);
342: return 0;
343: }
345: static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y)
346: {
347: PC_HPDDM *data = (PC_HPDDM *)pc->data;
349: PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
351: KSPMatSolve(data->levels[0]->ksp, X, Y);
352: return 0;
353: }
355: /*@C
356: PCHPDDMGetComplexities - Computes the grid and operator complexities.
358: Input Parameter:
359: . pc - preconditioner context
361: Output Parameters:
362: + gc - grid complexity = sum_i(m_i) / m_1
363: - oc - operator complexity = sum_i(nnz_i) / nnz_1
365: Level: advanced
367: .seealso: `PCMGGetGridComplexity()`, `PCHPDDM`, `PCHYPRE`, `PCGAMG`
368: @*/
369: static PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc)
370: {
371: PC_HPDDM *data = (PC_HPDDM *)pc->data;
372: MatInfo info;
373: PetscInt n, m;
374: PetscLogDouble accumulate[2]{}, nnz1 = 1.0, m1 = 1.0;
376: for (n = 0, *gc = 0, *oc = 0; n < data->N; ++n) {
377: if (data->levels[n]->ksp) {
378: Mat P, A;
379: KSPGetOperators(data->levels[n]->ksp, NULL, &P);
380: MatGetSize(P, &m, NULL);
381: accumulate[0] += m;
382: if (n == 0) {
383: PetscBool flg;
384: PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "");
385: if (flg) {
386: MatConvert(P, MATAIJ, MAT_INITIAL_MATRIX, &A);
387: P = A;
388: } else PetscObjectReference((PetscObject)P);
389: }
390: if (P->ops->getinfo) {
391: MatGetInfo(P, MAT_GLOBAL_SUM, &info);
392: accumulate[1] += info.nz_used;
393: }
394: if (n == 0) {
395: m1 = m;
396: if (P->ops->getinfo) nnz1 = info.nz_used;
397: MatDestroy(&P);
398: }
399: }
400: }
401: *gc = static_cast<PetscReal>(accumulate[0] / m1);
402: *oc = static_cast<PetscReal>(accumulate[1] / nnz1);
403: return 0;
404: }
406: static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer)
407: {
408: PC_HPDDM *data = (PC_HPDDM *)pc->data;
409: PetscViewer subviewer;
410: PetscSubcomm subcomm;
411: PetscReal oc, gc;
412: PetscInt i, tabs;
413: PetscMPIInt size, color, rank;
414: PetscBool ascii;
416: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
417: if (ascii) {
418: PetscViewerASCIIPrintf(viewer, "level%s: %" PetscInt_FMT "\n", data->N > 1 ? "s" : "", data->N);
419: PCHPDDMGetComplexities(pc, &gc, &oc);
420: if (data->N > 1) {
421: if (!data->deflation) {
422: PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[data->Neumann]);
423: PetscViewerASCIIPrintf(viewer, "shared subdomain KSP between SLEPc and PETSc? %s\n", PetscBools[data->share]);
424: } else PetscViewerASCIIPrintf(viewer, "user-supplied deflation matrix\n");
425: PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction]);
426: PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "");
427: PetscViewerASCIIGetTab(viewer, &tabs);
428: PetscViewerASCIISetTab(viewer, 0);
429: for (i = 1; i < data->N; ++i) {
430: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, data->levels[i - 1]->nu);
431: if (data->levels[i - 1]->threshold > -0.1) PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold);
432: }
433: PetscViewerASCIIPrintf(viewer, "\n");
434: PetscViewerASCIISetTab(viewer, tabs);
435: }
436: PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc);
437: if (data->levels[0]->ksp) {
438: KSPView(data->levels[0]->ksp, viewer);
439: if (data->levels[0]->pc) PCView(data->levels[0]->pc, viewer);
440: for (i = 1; i < data->N; ++i) {
441: if (data->levels[i]->ksp) color = 1;
442: else color = 0;
443: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
444: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
445: PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm);
446: PetscSubcommSetNumber(subcomm, PetscMin(size, 2));
447: PetscSubcommSetTypeGeneral(subcomm, color, rank);
448: PetscViewerASCIIPushTab(viewer);
449: PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
450: if (color == 1) {
451: KSPView(data->levels[i]->ksp, subviewer);
452: if (data->levels[i]->pc) PCView(data->levels[i]->pc, subviewer);
453: PetscViewerFlush(subviewer);
454: }
455: PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer);
456: PetscViewerASCIIPopTab(viewer);
457: PetscSubcommDestroy(&subcomm);
458: PetscViewerFlush(viewer);
459: }
460: }
461: }
462: return 0;
463: }
465: static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec)
466: {
467: PC_HPDDM *data = (PC_HPDDM *)pc->data;
468: PetscBool flg;
469: Mat A;
471: if (ksp) {
472: PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg);
473: if (flg && !data->normal) {
474: KSPGetOperators(ksp, &A, NULL);
475: MatCreateVecs(A, NULL, &data->normal); /* temporary Vec used in PCApply_HPDDMShell() for coarse grid corrections */
476: }
477: }
478: return 0;
479: }
481: static PetscErrorCode PCSetUp_HPDDMShell(PC pc)
482: {
483: PC_HPDDM_Level *ctx;
484: Mat A, P;
485: Vec x;
486: const char *pcpre;
488: PCShellGetContext(pc, &ctx);
489: KSPGetOptionsPrefix(ctx->ksp, &pcpre);
490: KSPGetOperators(ctx->ksp, &A, &P);
491: /* smoother */
492: PCSetOptionsPrefix(ctx->pc, pcpre);
493: PCSetOperators(ctx->pc, A, P);
494: if (!ctx->v[0]) {
495: VecDuplicateVecs(ctx->D, 1, &ctx->v[0]);
496: if (!std::is_same<PetscScalar, PetscReal>::value) VecDestroy(&ctx->D);
497: MatCreateVecs(A, &x, NULL);
498: VecDuplicateVecs(x, 2, &ctx->v[1]);
499: VecDestroy(&x);
500: }
501: std::fill_n(ctx->V, 3, nullptr);
502: return 0;
503: }
505: template <class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
506: static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y)
507: {
508: PC_HPDDM_Level *ctx;
509: PetscScalar *out;
511: PCShellGetContext(pc, &ctx);
512: /* going from PETSc to HPDDM numbering */
513: VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
514: VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD);
515: VecGetArrayWrite(ctx->v[0][0], &out);
516: ctx->P->deflation<false>(NULL, out, 1); /* y = Q x */
517: VecRestoreArrayWrite(ctx->v[0][0], &out);
518: /* going from HPDDM to PETSc numbering */
519: VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
520: VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE);
521: return 0;
522: }
524: template <class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
525: static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y)
526: {
527: PC_HPDDM_Level *ctx;
528: Vec vX, vY, vC;
529: PetscScalar *out;
530: PetscInt i, N;
532: PCShellGetContext(pc, &ctx);
533: MatGetSize(X, NULL, &N);
534: /* going from PETSc to HPDDM numbering */
535: for (i = 0; i < N; ++i) {
536: MatDenseGetColumnVecRead(X, i, &vX);
537: MatDenseGetColumnVecWrite(ctx->V[0], i, &vC);
538: VecScatterBegin(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD);
539: VecScatterEnd(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD);
540: MatDenseRestoreColumnVecWrite(ctx->V[0], i, &vC);
541: MatDenseRestoreColumnVecRead(X, i, &vX);
542: }
543: MatDenseGetArrayWrite(ctx->V[0], &out);
544: ctx->P->deflation<false>(NULL, out, N); /* Y = Q X */
545: MatDenseRestoreArrayWrite(ctx->V[0], &out);
546: /* going from HPDDM to PETSc numbering */
547: for (i = 0; i < N; ++i) {
548: MatDenseGetColumnVecRead(ctx->V[0], i, &vC);
549: MatDenseGetColumnVecWrite(Y, i, &vY);
550: VecScatterBegin(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE);
551: VecScatterEnd(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE);
552: MatDenseRestoreColumnVecWrite(Y, i, &vY);
553: MatDenseRestoreColumnVecRead(ctx->V[0], i, &vC);
554: }
555: return 0;
556: }
558: /*
559: PCApply_HPDDMShell - 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.
561: .vb
562: (1) y = Pmat^-1 x + Q x,
563: (2) y = Pmat^-1 (I - Amat Q) x + Q x (default),
564: (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x.
565: .ve
567: Input Parameters:
568: + pc - preconditioner context
569: - x - input vector
571: Output Parameter:
572: . y - output vector
574: Notes:
575: 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.
576: 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.
577: (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.
579: Level: advanced
581: Developer Note:
582: Since this is not an actual manual page the material below should be moved to an appropriate manual page with the appropriate context, i.e. explaining when it is used and how
583: to trigger it. Likely the manual page is `PCHPDDM`
585: .seealso: `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
586: */
587: static PetscErrorCode PCApply_HPDDMShell(PC pc, Vec x, Vec y)
588: {
589: PC_HPDDM_Level *ctx;
590: Mat A;
592: PCShellGetContext(pc, &ctx);
594: KSPGetOperators(ctx->ksp, &A, NULL);
595: PCHPDDMDeflate_Private(pc, x, y); /* y = Q x */
596: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
597: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) MatMult(A, y, ctx->v[1][0]); /* y = A Q x */
598: else { /* KSPLSQR and finest level */ MatMult(A, y, ctx->parent->normal); /* y = A Q x */
599: MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][0]); /* y = A^T A Q x */
600: }
601: VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x); /* y = (I - A Q) x */
602: PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0]); /* y = M^-1 (I - A Q) x */
603: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
604: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1]); /* z = A^T y */
605: else {
606: MatMult(A, ctx->v[1][0], ctx->parent->normal);
607: MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][1]); /* z = A^T A y */
608: }
609: PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1]);
610: VecAXPBYPCZ(y, -1.0, 1.0, 1.0, ctx->v[1][1], ctx->v[1][0]); /* y = (I - Q A^T) y + Q x */
611: } else VecAXPY(y, 1.0, ctx->v[1][0]); /* y = Q M^-1 (I - A Q) x + Q x */
612: } else {
614: PCApply(ctx->pc, x, ctx->v[1][0]);
615: VecAXPY(y, 1.0, ctx->v[1][0]); /* y = M^-1 x + Q x */
616: }
617: return 0;
618: }
620: /*
621: PCMatApply_HPDDMShell - Variant of PCApply_HPDDMShell() for blocks of vectors.
623: Input Parameters:
624: + pc - preconditioner context
625: - X - block of input vectors
627: Output Parameter:
628: . Y - block of output vectors
630: Level: advanced
632: .seealso: `PCHPDDM`, `PCApply_HPDDMShell()`, `PCHPDDMCoarseCorrectionType`
633: */
634: static PetscErrorCode PCMatApply_HPDDMShell(PC pc, Mat X, Mat Y)
635: {
636: PC_HPDDM_Level *ctx;
637: Mat A, *ptr;
638: PetscContainer container = NULL;
639: PetscScalar *array;
640: PetscInt m, M, N, prev = 0;
641: PetscBool reset = PETSC_FALSE;
643: PCShellGetContext(pc, &ctx);
645: MatGetSize(X, NULL, &N);
646: KSPGetOperators(ctx->ksp, &A, NULL);
647: PetscObjectQuery((PetscObject)A, "_HPDDM_MatProduct", (PetscObject *)&container);
648: if (container) { /* MatProduct container already attached */
649: PetscContainerGetPointer(container, (void **)&ptr);
650: if (ptr[1] != ctx->V[2]) /* Mat has changed or may have been set first in KSPHPDDM */
651: for (m = 0; m < 2; ++m) {
652: MatDestroy(ctx->V + m + 1);
653: ctx->V[m + 1] = ptr[m];
654: PetscObjectReference((PetscObject)ctx->V[m + 1]);
655: }
656: }
657: if (ctx->V[1]) MatGetSize(ctx->V[1], NULL, &prev);
658: if (N != prev || !ctx->V[0]) {
659: MatDestroy(ctx->V);
660: VecGetLocalSize(ctx->v[0][0], &m);
661: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, ctx->V);
662: if (N != prev) {
663: MatDestroy(ctx->V + 1);
664: MatDestroy(ctx->V + 2);
665: MatGetLocalSize(X, &m, NULL);
666: MatGetSize(X, &M, NULL);
667: if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) MatDenseGetArrayWrite(ctx->V[0], &array);
668: else array = NULL;
669: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, array, ctx->V + 1);
670: if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) MatDenseRestoreArrayWrite(ctx->V[0], &array);
671: MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, NULL, ctx->V + 2);
672: MatProductCreateWithMat(A, Y, NULL, ctx->V[1]);
673: MatProductSetType(ctx->V[1], MATPRODUCT_AB);
674: MatProductSetFromOptions(ctx->V[1]);
675: MatProductSymbolic(ctx->V[1]);
676: if (!container) { /* no MatProduct container attached, create one to be queried in KSPHPDDM or at the next call to PCMatApply() */
677: PetscContainerCreate(PetscObjectComm((PetscObject)A), &container);
678: PetscObjectCompose((PetscObject)A, "_HPDDM_MatProduct", (PetscObject)container);
679: }
680: PetscContainerSetPointer(container, ctx->V + 1); /* need to compose B and D from MatProductCreateWithMath(A, B, NULL, D), which are stored in the contiguous array ctx->V */
681: }
682: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
683: MatProductCreateWithMat(A, ctx->V[1], NULL, ctx->V[2]);
684: MatProductSetType(ctx->V[2], MATPRODUCT_AtB);
685: MatProductSetFromOptions(ctx->V[2]);
686: MatProductSymbolic(ctx->V[2]);
687: }
688: ctx->P->start(N);
689: }
690: if (N == prev || container) { /* when MatProduct container is attached, always need to MatProductReplaceMats() since KSPHPDDM may have replaced the Mat as well */
691: MatProductReplaceMats(NULL, Y, NULL, ctx->V[1]);
692: if (container && ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) {
693: MatDenseGetArrayWrite(ctx->V[0], &array);
694: MatDensePlaceArray(ctx->V[1], array);
695: MatDenseRestoreArrayWrite(ctx->V[0], &array);
696: reset = PETSC_TRUE;
697: }
698: }
699: PCHPDDMDeflate_Private(pc, X, Y);
700: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
701: MatProductNumeric(ctx->V[1]);
702: MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN);
703: MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN);
704: PCMatApply(ctx->pc, ctx->V[2], ctx->V[1]);
705: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
706: MatProductNumeric(ctx->V[2]);
707: PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2]);
708: MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN);
709: }
710: MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN);
711: } else {
713: PCMatApply(ctx->pc, X, ctx->V[1]);
714: MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN);
715: }
716: if (reset) MatDenseResetArray(ctx->V[1]);
717: return 0;
718: }
720: static PetscErrorCode PCDestroy_HPDDMShell(PC pc)
721: {
722: PC_HPDDM_Level *ctx;
723: PetscContainer container;
725: PCShellGetContext(pc, &ctx);
726: HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE);
727: VecDestroyVecs(1, &ctx->v[0]);
728: VecDestroyVecs(2, &ctx->v[1]);
729: PetscObjectQuery((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", (PetscObject *)&container);
730: PetscContainerDestroy(&container);
731: PetscObjectCompose((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", NULL);
732: MatDestroy(ctx->V);
733: MatDestroy(ctx->V + 1);
734: MatDestroy(ctx->V + 2);
735: VecDestroy(&ctx->D);
736: VecScatterDestroy(&ctx->scatter);
737: PCDestroy(&ctx->pc);
738: return 0;
739: }
741: static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu)
742: {
743: Mat B, X;
744: PetscInt n, N, j = 0;
746: KSPGetOperators(ctx->ksp, &B, NULL);
747: MatGetLocalSize(B, &n, NULL);
748: MatGetSize(B, &N, NULL);
749: if (ctx->parent->log_separate) {
750: j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx));
751: PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, 0, 0, 0);
752: }
753: if (mu == 1) {
754: if (!ctx->ksp->vec_rhs) {
755: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, NULL, &ctx->ksp->vec_rhs);
756: VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol);
757: }
758: VecPlaceArray(ctx->ksp->vec_rhs, rhs);
759: KSPSolve(ctx->ksp, NULL, NULL);
760: VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs);
761: VecResetArray(ctx->ksp->vec_rhs);
762: } else {
763: MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B);
764: MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, NULL, &X);
765: KSPMatSolve(ctx->ksp, B, X);
766: MatCopy(X, B, SAME_NONZERO_PATTERN);
767: MatDestroy(&X);
768: MatDestroy(&B);
769: }
770: if (ctx->parent->log_separate) PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, 0, 0, 0);
771: return 0;
772: }
774: static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc)
775: {
776: PC_HPDDM *data = (PC_HPDDM *)pc->data;
778: if (data->setup) {
779: Mat P;
780: Vec x, xt = NULL;
781: PetscReal t = 0.0, s = 0.0;
783: PCGetOperators(pc, NULL, &P);
784: PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x);
785: PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx));
786: }
787: return 0;
788: }
790: static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[])
791: {
792: Mat A;
795: /* previously composed Mat */
796: PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A);
798: if (scall == MAT_INITIAL_MATRIX) {
799: PetscCalloc1(1, submat);
800: MatDuplicate(A, MAT_COPY_VALUES, *submat);
801: } else MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN);
802: return 0;
803: }
805: static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted)
806: {
807: void (*op)(void);
809: /* previously-composed Mat */
810: PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C);
811: MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op);
812: /* trick suggested by Barry https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-January/025491.html */
813: MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (void (*)(void))PCHPDDMCreateSubMatrices_Private);
814: if (sorted) PCASMSetSortIndices(pc, PETSC_FALSE); /* everything is already sorted */
815: PCSetFromOptions(pc); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */
816: PCSetUp(pc);
817: /* reset MatCreateSubMatrices() */
818: MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op);
819: PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", NULL);
820: return 0;
821: }
823: static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p)
824: {
825: IS perm;
826: const PetscInt *ptr;
827: PetscInt *concatenate, size, n;
828: std::map<PetscInt, PetscInt> order;
829: PetscBool sorted;
831: ISSorted(is, &sorted);
832: if (!sorted) {
833: ISGetLocalSize(is, &size);
834: ISGetIndices(is, &ptr);
835: /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */
836: for (n = 0; n < size; ++n) order.insert(std::make_pair(ptr[n], n));
837: ISRestoreIndices(is, &ptr);
838: if (out_C) {
839: PetscMalloc1(size, &concatenate);
840: for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second;
841: concatenate -= size;
842: ISCreateGeneral(PetscObjectComm((PetscObject)in_C), size, concatenate, PETSC_OWN_POINTER, &perm);
843: ISSetPermutation(perm);
844: /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */
845: MatPermute(in_C, perm, perm, out_C);
846: if (p) *p = perm;
847: else ISDestroy(&perm); /* no need to save the permutation */
848: }
849: if (out_is) {
850: PetscMalloc1(size, &concatenate);
851: for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first;
852: concatenate -= size;
853: /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */
854: ISCreateGeneral(PetscObjectComm((PetscObject)in_is), size, concatenate, PETSC_OWN_POINTER, out_is);
855: }
856: } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */
857: if (out_C) MatDuplicate(in_C, MAT_COPY_VALUES, out_C);
858: if (out_is) ISDuplicate(in_is, out_is);
859: }
860: return 0;
861: }
863: static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub)
864: {
865: IS is;
867: if (!flg) {
868: if (algebraic) {
869: PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is);
870: ISDestroy(&is);
871: PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", NULL);
872: PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", NULL);
873: }
874: MatDestroySubMatrices(algebraic ? 2 : 1, &sub);
875: }
876: return 0;
877: }
879: static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block)
880: {
881: IS icol[3], irow[2];
882: Mat *M, Q;
883: PetscReal *ptr;
884: PetscInt *idx, p = 0, n, bs = PetscAbs(P->cmap->bs);
885: PetscBool flg;
887: ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2);
888: ISSetBlockSize(icol[2], bs);
889: ISSetIdentity(icol[2]);
890: PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg);
891: if (flg) {
892: /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */
893: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q);
894: std::swap(P, Q);
895: }
896: MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M);
897: if (flg) {
898: std::swap(P, Q);
899: MatDestroy(&Q);
900: }
901: ISDestroy(icol + 2);
902: ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow);
903: ISSetBlockSize(irow[0], bs);
904: ISSetIdentity(irow[0]);
905: if (!block) {
906: PetscMalloc2(P->cmap->N, &ptr, P->cmap->N, &idx);
907: MatGetColumnNorms(M[0], NORM_INFINITY, ptr);
908: /* check for nonzero columns so that M[0] may be expressed in compact form */
909: for (n = 0; n < P->cmap->N; n += bs)
910: if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) {
911: std::iota(idx + p, idx + p + bs, n);
912: p += bs;
913: }
914: ISCreateGeneral(PETSC_COMM_SELF, p, idx, PETSC_USE_POINTER, icol + 1);
915: ISSetBlockSize(icol[1], bs);
916: ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE);
917: ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2);
918: irow[1] = irow[0];
919: /* first Mat will be used in PCASM (if it is used as a PC on this level) and as the left-hand side of GenEO */
920: icol[0] = is[0];
921: MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub);
922: ISDestroy(icol + 1);
923: PetscFree2(ptr, idx);
924: /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */
925: PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2]);
926: /* Mat used in eq. (3.1) of [2022b] */
927: PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1]);
928: } else {
929: Mat aux;
930: MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
931: /* diagonal block of the overlapping rows */
932: MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub);
933: MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux);
934: MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
935: if (bs == 1) { /* scalar case */
936: Vec sum[2];
937: MatCreateVecs(aux, sum, sum + 1);
938: MatGetRowSum(M[0], sum[0]);
939: MatGetRowSum(aux, sum[1]);
940: /* off-diagonal block row sum (full rows - diagonal block rows) */
941: VecAXPY(sum[0], -1.0, sum[1]);
942: /* subdomain matrix plus off-diagonal block row sum */
943: MatDiagonalSet(aux, sum[0], ADD_VALUES);
944: VecDestroy(sum);
945: VecDestroy(sum + 1);
946: } else { /* vectorial case */
947: /* TODO: missing MatGetValuesBlocked(), so the code below is */
948: /* an extension of the scalar case for when bs > 1, but it could */
949: /* be more efficient by avoiding all these MatMatMult() */
950: Mat sum[2], ones;
951: PetscScalar *ptr;
952: PetscCalloc1(M[0]->cmap->n * bs, &ptr);
953: MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones);
954: for (n = 0; n < M[0]->cmap->n; n += bs) {
955: for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0;
956: }
957: MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum);
958: MatDestroy(&ones);
959: MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones);
960: MatDenseSetLDA(ones, M[0]->cmap->n);
961: MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum + 1);
962: MatDestroy(&ones);
963: PetscFree(ptr);
964: /* off-diagonal block row sum (full rows - diagonal block rows) */
965: MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN);
966: MatDestroy(sum + 1);
967: /* re-order values to be consistent with MatSetValuesBlocked() */
968: /* equivalent to MatTranspose() which does not truly handle */
969: /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */
970: MatDenseGetArrayWrite(sum[0], &ptr);
971: HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs);
972: /* subdomain matrix plus off-diagonal block row sum */
973: for (n = 0; n < aux->cmap->n / bs; ++n) MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES);
974: MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY);
975: MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY);
976: MatDenseRestoreArrayWrite(sum[0], &ptr);
977: MatDestroy(sum);
978: }
979: MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
980: /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */
981: PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux);
982: }
983: ISDestroy(irow);
984: MatDestroySubMatrices(1, &M);
985: return 0;
986: }
988: static PetscErrorCode PCSetUp_HPDDM(PC pc)
989: {
990: PC_HPDDM *data = (PC_HPDDM *)pc->data;
991: PC inner;
992: KSP *ksp;
993: Mat *sub, A, P, N, C = NULL, uaux = NULL, weighted, subA[2];
994: Vec xin, v;
995: std::vector<Vec> initial;
996: IS is[1], loc, uis = data->is;
997: ISLocalToGlobalMapping l2g;
998: char prefix[256];
999: const char *pcpre;
1000: const PetscScalar *const *ev;
1001: PetscInt n, requested = data->N, reused = 0;
1002: MatStructure structure = UNKNOWN_NONZERO_PATTERN;
1003: PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE;
1004: DM dm;
1007: PCGetOptionsPrefix(pc, &pcpre);
1008: PCGetOperators(pc, &A, &P);
1009: if (!data->levels[0]->ksp) {
1010: KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp);
1011: PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse");
1012: KSPSetOptionsPrefix(data->levels[0]->ksp, prefix);
1013: KSPSetType(data->levels[0]->ksp, KSPPREONLY);
1014: } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) {
1015: /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */
1016: /* then just propagate the appropriate flag to the coarser levels */
1017: for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1018: /* the following KSP and PC may be NULL for some processes, hence the check */
1019: if (data->levels[n]->ksp) KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE);
1020: if (data->levels[n]->pc) PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE);
1021: }
1022: /* early bail out because there is nothing to do */
1023: return 0;
1024: } else {
1025: /* reset coarser levels */
1026: for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1027: 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) {
1028: reused = data->N - n;
1029: break;
1030: }
1031: KSPDestroy(&data->levels[n]->ksp);
1032: PCDestroy(&data->levels[n]->pc);
1033: }
1034: /* check if some coarser levels are being reused */
1035: MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc));
1036: const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0;
1038: if (addr != &HPDDM::i__0 && reused != data->N - 1) {
1039: /* reuse previously computed eigenvectors */
1040: ev = data->levels[0]->P->getVectors();
1041: if (ev) {
1042: initial.reserve(*addr);
1043: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin);
1044: for (n = 0; n < *addr; ++n) {
1045: VecDuplicate(xin, &v);
1046: VecPlaceArray(xin, ev[n]);
1047: VecCopy(xin, v);
1048: initial.emplace_back(v);
1049: VecResetArray(xin);
1050: }
1051: VecDestroy(&xin);
1052: }
1053: }
1054: }
1055: data->N -= reused;
1056: KSPSetOperators(data->levels[0]->ksp, A, P);
1058: PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis);
1059: if (!data->is && !ismatis) {
1060: PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode(**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = NULL;
1061: PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = NULL;
1062: void *uctx = NULL;
1064: /* first see if we can get the data from the DM */
1065: MatGetDM(P, &dm);
1066: if (!dm) MatGetDM(A, &dm);
1067: if (!dm) PCGetDM(pc, &dm);
1068: if (dm) { /* this is the hook for DMPLEX and DMDA for which the auxiliary Mat is the local Neumann matrix */
1069: PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create);
1070: if (create) {
1071: (*create)(dm, &uis, &uaux, &usetup, &uctx);
1072: data->Neumann = PETSC_TRUE;
1073: }
1074: }
1075: if (!create) {
1076: if (!uis) {
1077: PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis);
1078: PetscObjectReference((PetscObject)uis);
1079: }
1080: if (!uaux) {
1081: PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux);
1082: PetscObjectReference((PetscObject)uaux);
1083: }
1084: /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */
1085: if (!uis) {
1086: PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis);
1087: PetscObjectReference((PetscObject)uis);
1088: }
1089: if (!uaux) {
1090: PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux);
1091: PetscObjectReference((PetscObject)uaux);
1092: }
1093: }
1094: PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx);
1095: MatDestroy(&uaux);
1096: ISDestroy(&uis);
1097: }
1099: if (!ismatis) {
1100: PCHPDDMSetUpNeumannOverlap_Private(pc);
1101: PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_block_splitting", &block, NULL);
1102: if (data->is && block) {
1103: ISDestroy(&data->is);
1104: MatDestroy(&data->aux);
1105: }
1106: if (!data->is && data->N > 1) {
1107: char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */
1108: PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "");
1109: if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) {
1110: Mat B;
1111: PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre);
1112: if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED;
1113: MatDestroy(&B);
1114: } else {
1115: PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg);
1116: if (flg) {
1117: Mat A00, P00, A01, A10, A11, B, N;
1118: const PetscScalar *array;
1119: PetscReal norm;
1120: MatSchurComplementAinvType type;
1122: MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11);
1123: if (A11) {
1124: MatNorm(A11, NORM_INFINITY, &norm);
1126: }
1127: if (PetscDefined(USE_DEBUG)) {
1128: Mat T, U = NULL;
1129: IS z;
1130: PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg);
1131: if (flg) MatTransposeGetMat(A10, &U);
1132: else {
1133: PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg);
1134: if (flg) MatHermitianTransposeGetMat(A10, &U);
1135: }
1136: if (U) MatDuplicate(U, MAT_COPY_VALUES, &T);
1137: else MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T);
1138: PetscLayoutCompare(T->rmap, A01->rmap, &flg);
1139: if (flg) {
1140: PetscLayoutCompare(T->cmap, A01->cmap, &flg);
1141: if (flg) {
1142: MatFindZeroRows(A01, &z); /* for essential boundary conditions, some implementations will */
1143: if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */
1144: MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
1145: MatZeroRowsIS(T, z, 0.0, NULL, NULL); /* corresponding zero rows from A01 */
1146: ISDestroy(&z);
1147: }
1148: MatMultEqual(A01, T, 10, &flg);
1150: } else PetscInfo(pc, "A01 and A10^T have non-congruent column layouts, cannot test for equality\n");
1151: }
1152: MatDestroy(&T);
1153: }
1154: MatCreateVecs(P00, &v, NULL);
1155: MatSchurComplementGetAinvType(P, &type);
1157: if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
1158: MatGetRowSum(P00, v);
1159: if (A00 == P00) PetscObjectReference((PetscObject)A00);
1160: MatDestroy(&P00);
1161: VecGetArrayRead(v, &array);
1162: MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, NULL, 0, NULL, &P00);
1163: MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
1164: for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES);
1165: MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY);
1166: MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY);
1167: VecRestoreArrayRead(v, &array);
1168: MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11); /* replace P00 by diag(sum of each row of P00) */
1169: MatDestroy(&P00);
1170: } else MatGetDiagonal(P00, v);
1171: VecReciprocal(v); /* inv(diag(P00)) */
1172: VecSqrtAbs(v); /* sqrt(inv(diag(P00))) */
1173: MatDuplicate(A01, MAT_COPY_VALUES, &B);
1174: MatDiagonalScale(B, v, NULL);
1175: VecDestroy(&v);
1176: MatCreateNormalHermitian(B, &N);
1177: PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre);
1178: PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg);
1179: if (!flg) {
1180: MatDestroy(&P);
1181: P = N;
1182: PetscObjectReference((PetscObject)P);
1183: } else MatScale(P, -1.0);
1184: MatScale(N, -1.0);
1185: PCSetOperators(pc, N, P); /* replace P by -A01' inv(diag(P00)) A01 */
1186: MatDestroy(&N);
1187: MatDestroy(&P);
1188: MatDestroy(&B);
1189: return 0;
1190: } else {
1191: PetscOptionsGetString(NULL, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), NULL);
1192: PetscStrcmp(type, PCMAT, &algebraic);
1193: PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_block_splitting", &block, NULL);
1195: if (block) algebraic = PETSC_TRUE;
1196: if (algebraic) {
1197: ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is);
1198: MatIncreaseOverlap(P, 1, &data->is, 1);
1199: ISSort(data->is);
1200: } else PetscInfo(pc, "Cannot assemble a fully-algebraic coarse operator with an assembled Pmat and -%spc_hpddm_levels_1_st_pc_type != mat and -%spc_hpddm_block_splitting != true\n", pcpre ? pcpre : "", pcpre ? pcpre : "");
1201: }
1202: }
1203: }
1204: }
1206: if (data->is || (ismatis && data->N > 1)) {
1207: if (ismatis) {
1208: std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ};
1209: MatISGetLocalMat(P, &N);
1210: std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name);
1211: MatISRestoreLocalMat(P, &N);
1212: switch (std::distance(list.begin(), it)) {
1213: case 0:
1214: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
1215: break;
1216: case 1:
1217: /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */
1218: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C);
1219: MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE);
1220: break;
1221: default:
1222: MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C);
1223: }
1224: MatISGetLocalToGlobalMapping(P, &l2g, NULL);
1225: PetscObjectReference((PetscObject)P);
1226: KSPSetOperators(data->levels[0]->ksp, A, C);
1227: std::swap(C, P);
1228: ISLocalToGlobalMappingGetSize(l2g, &n);
1229: ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc);
1230: ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0]);
1231: ISDestroy(&loc);
1232: /* the auxiliary Mat is _not_ the local Neumann matrix */
1233: /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */
1234: data->Neumann = PETSC_FALSE;
1235: structure = SAME_NONZERO_PATTERN;
1236: if (data->share) {
1237: data->share = PETSC_FALSE;
1238: PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n");
1239: }
1240: } else {
1241: is[0] = data->is;
1242: if (algebraic) subdomains = PETSC_TRUE;
1243: PetscOptionsGetBool(NULL, pcpre, "-pc_hpddm_define_subdomains", &subdomains, NULL);
1244: if (data->share) {
1245: if (!subdomains) {
1246: PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "");
1247: data->share = PETSC_FALSE;
1248: }
1249: if (data->deflation) {
1250: PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n");
1251: data->share = PETSC_FALSE;
1252: }
1253: }
1254: if (data->Neumann) {
1257: }
1258: if (data->Neumann || block) structure = SAME_NONZERO_PATTERN;
1259: ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc);
1260: }
1261: PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "");
1262: PetscOptionsGetEnum(NULL, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg); /* if not user-provided, force its value when possible */
1263: if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */
1264: PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "");
1265: PetscOptionsSetValue(NULL, prefix, MatStructures[structure]);
1266: }
1267: if (data->N > 1 && (data->aux || ismatis || algebraic)) {
1269: MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
1270: if (ismatis) {
1271: /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */
1272: MatIncreaseOverlap(P, 1, is, 1);
1273: ISDestroy(&data->is);
1274: data->is = is[0];
1275: } else {
1276: if (PetscDefined(USE_DEBUG)) {
1277: PetscBool equal;
1278: IS intersect;
1280: ISIntersect(data->is, loc, &intersect);
1281: ISEqualUnsorted(loc, intersect, &equal);
1282: ISDestroy(&intersect);
1284: }
1285: PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private);
1286: if (!data->Neumann && !algebraic) {
1287: PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg);
1288: if (flg) {
1289: /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */
1290: /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */
1291: MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux);
1292: flg = PETSC_FALSE;
1293: }
1294: }
1295: }
1296: if (algebraic) {
1297: PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block));
1298: if (block) {
1299: PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux);
1300: PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", NULL);
1301: }
1302: } else if (!uaux) {
1303: if (data->Neumann) sub = &data->aux;
1304: else MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub);
1305: } else {
1306: MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub);
1307: MatDestroy(&uaux);
1308: MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub);
1309: }
1310: /* Vec holding the partition of unity */
1311: if (!data->levels[0]->D) {
1312: ISGetLocalSize(data->is, &n);
1313: VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D);
1314: }
1315: if (data->share && structure == SAME_NONZERO_PATTERN) { /* share the KSP only when the MatStructure is SAME_NONZERO_PATTERN */
1316: Mat D;
1317: IS perm = NULL;
1318: PetscInt size = -1;
1319: PCHPDDMPermute_Private(*is, data->is, &uis, data->Neumann || block ? sub[0] : data->aux, &C, &perm);
1320: if (!data->Neumann && !block) {
1321: MatPermute(sub[0], perm, perm, &D); /* permute since PCASM will call ISSort() */
1322: MatHeaderReplace(sub[0], &D);
1323: }
1324: if (data->B) { /* see PCHPDDMSetRHSMat() */
1325: MatPermute(data->B, perm, perm, &D);
1326: MatHeaderReplace(data->B, &D);
1327: }
1328: ISDestroy(&perm);
1329: if (!data->levels[0]->pc) {
1330: PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "");
1331: PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc);
1332: PCSetOptionsPrefix(data->levels[0]->pc, prefix);
1333: PCSetOperators(data->levels[0]->pc, A, P);
1334: }
1335: PCSetType(data->levels[0]->pc, PCASM);
1336: if (!data->levels[0]->pc->setupcalled) PCASMSetLocalSubdomains(data->levels[0]->pc, 1, is, &loc);
1337: PCSetFromOptions(data->levels[0]->pc);
1338: if (block) PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic);
1339: else PCSetUp(data->levels[0]->pc);
1340: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, NULL, &ksp));
1341: if (size != 1) {
1342: PCDestroy(&data->levels[0]->pc);
1343: MatDestroy(&C);
1344: ISDestroy(&uis);
1345: data->share = PETSC_FALSE;
1347: PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n");
1348: } else {
1349: const char *matpre;
1350: PetscBool cmp[4];
1351: KSPGetOperators(ksp[0], subA, subA + 1);
1352: MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D);
1353: MatGetOptionsPrefix(subA[1], &matpre);
1354: MatSetOptionsPrefix(D, matpre);
1355: PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp);
1356: PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1);
1357: if (!cmp[0]) PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2);
1358: else cmp[2] = PETSC_FALSE;
1359: if (!cmp[1]) PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3);
1360: else cmp[3] = PETSC_FALSE;
1362: if (!cmp[0] && !cmp[2]) {
1363: if (!block) MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN);
1364: else {
1365: MatMissingDiagonal(D, cmp, NULL);
1366: if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */
1367: MatAXPY(D, 1.0, data->aux, structure);
1368: }
1369: } else {
1370: Mat mat[2];
1371: if (cmp[0]) {
1372: MatNormalGetMat(D, mat);
1373: MatNormalGetMat(C, mat + 1);
1374: } else {
1375: MatNormalHermitianGetMat(D, mat);
1376: MatNormalHermitianGetMat(C, mat + 1);
1377: }
1378: MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN);
1379: }
1380: MatPropagateSymmetryOptions(C, D);
1381: MatDestroy(&C);
1382: C = D;
1383: /* swap pointers so that variables stay consistent throughout PCSetUp() */
1384: std::swap(C, data->aux);
1385: std::swap(uis, data->is);
1386: swap = PETSC_TRUE;
1387: }
1388: } else if (data->share) {
1389: data->share = PETSC_FALSE;
1390: PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_levels_1_st_matstructure %s (!= %s)\n", pcpre ? pcpre : "", MatStructures[structure], MatStructures[SAME_NONZERO_PATTERN]);
1391: }
1392: if (!data->levels[0]->scatter) {
1393: MatCreateVecs(P, &xin, NULL);
1394: if (ismatis) MatDestroy(&P);
1395: VecScatterCreate(xin, data->is, data->levels[0]->D, NULL, &data->levels[0]->scatter);
1396: VecDestroy(&xin);
1397: }
1398: if (data->levels[0]->P) {
1399: /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */
1400: HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE);
1401: }
1402: if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>();
1403: if (data->log_separate) PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, 0, 0, 0);
1404: else PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, 0, 0, 0);
1405: /* HPDDM internal data structure */
1406: data->levels[0]->P->structure(loc, data->is, sub[0], ismatis ? C : data->aux, data->levels);
1407: if (!data->log_separate) PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, 0, 0, 0);
1408: /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */
1409: if (data->deflation) weighted = data->aux;
1410: else if (!data->B) {
1411: PetscBool cmp[2];
1412: MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted);
1413: PetscObjectTypeCompare((PetscObject)weighted, MATNORMAL, cmp);
1414: if (!cmp[0]) PetscObjectTypeCompare((PetscObject)weighted, MATNORMALHERMITIAN, cmp + 1);
1415: else cmp[1] = PETSC_FALSE;
1416: if (!cmp[0] && !cmp[1]) MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D);
1417: else { /* MATNORMAL applies MatDiagonalScale() in a matrix-free fashion, not what is needed since this won't be passed to SLEPc during the eigensolve */
1418: if (cmp[0]) MatNormalGetMat(weighted, &data->B);
1419: else MatNormalHermitianGetMat(weighted, &data->B);
1420: MatDiagonalScale(data->B, NULL, data->levels[0]->D);
1421: data->B = NULL;
1422: flg = PETSC_FALSE;
1423: }
1424: /* neither MatDuplicate() nor MatDiagonaleScale() handles the symmetry options, so propagate the options explicitly */
1425: /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */
1426: MatPropagateSymmetryOptions(sub[0], weighted);
1427: } else weighted = data->B;
1428: /* SLEPc is used inside the loaded symbol */
1429: (*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : (algebraic && !block ? sub[0] : data->aux), weighted, data->B, initial, data->levels);
1430: if (data->share) {
1431: Mat st[2];
1432: KSPGetOperators(ksp[0], st, st + 1);
1433: MatCopy(subA[0], st[0], structure);
1434: if (subA[1] != subA[0] || st[1] != st[0]) MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN);
1435: }
1436: if (data->log_separate) PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, 0, 0, 0);
1437: if (ismatis) MatISGetLocalMat(C, &N);
1438: else N = data->aux;
1439: P = sub[0];
1440: /* going through the grid hierarchy */
1441: for (n = 1; n < data->N; ++n) {
1442: if (data->log_separate) PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, 0, 0, 0);
1443: /* method composed in the loaded symbol since there, SLEPc is used as well */
1444: PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels));
1445: if (data->log_separate) PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, 0, 0, 0);
1446: }
1447: /* reset to NULL to avoid any faulty use */
1448: PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", NULL);
1449: if (!ismatis) PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", NULL);
1450: else PetscObjectDereference((PetscObject)C); /* matching PetscObjectReference() above */
1451: for (n = 0; n < data->N - 1; ++n)
1452: if (data->levels[n]->P) {
1453: /* HPDDM internal work buffers */
1454: data->levels[n]->P->setBuffer();
1455: data->levels[n]->P->super::start();
1456: }
1457: if (ismatis || !subdomains) PCHPDDMDestroySubMatrices_Private(data->Neumann, PetscBool(algebraic && !block), sub);
1458: if (ismatis) data->is = NULL;
1459: for (n = 0; n < data->N - 1 + (reused > 0); ++n) {
1460: if (data->levels[n]->P) {
1461: PC spc;
1463: /* force the PC to be PCSHELL to do the coarse grid corrections */
1464: KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE);
1465: KSPGetPC(data->levels[n]->ksp, &spc);
1466: PCSetType(spc, PCSHELL);
1467: PCShellSetContext(spc, data->levels[n]);
1468: PCShellSetSetUp(spc, PCSetUp_HPDDMShell);
1469: PCShellSetApply(spc, PCApply_HPDDMShell);
1470: PCShellSetMatApply(spc, PCMatApply_HPDDMShell);
1471: PCShellSetDestroy(spc, PCDestroy_HPDDMShell);
1472: if (!data->levels[n]->pc) PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc);
1473: if (n < reused) {
1474: PCSetReusePreconditioner(spc, PETSC_TRUE);
1475: PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE);
1476: }
1477: PCSetUp(spc);
1478: }
1479: }
1480: PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", NULL);
1481: } else flg = reused ? PETSC_FALSE : PETSC_TRUE;
1482: if (!ismatis && subdomains) {
1483: if (flg) KSPGetPC(data->levels[0]->ksp, &inner);
1484: else inner = data->levels[0]->pc;
1485: if (inner) {
1486: PCSetType(inner, PCASM); /* inner is the fine-level PC for which one must ensure */
1487: /* PCASMSetLocalSubdomains() has been called when -pc_hpddm_define_subdomains */
1488: if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */
1489: PCASMSetLocalSubdomains(inner, 1, is, &loc);
1490: if (!data->Neumann && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */
1491: PCHPDDMPermute_Private(*is, NULL, NULL, sub[0], &C, NULL);
1492: PCHPDDMCommunicationAvoidingPCASM_Private(inner, C, algebraic);
1493: MatDestroy(&C);
1494: }
1495: }
1496: }
1497: if (data->N > 1) PCHPDDMDestroySubMatrices_Private(data->Neumann, PetscBool(algebraic && !block), sub);
1498: }
1499: ISDestroy(&loc);
1500: } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */
1501: if (requested != data->N + reused) {
1502: PetscCall(PetscInfo(pc, "%" PetscInt_FMT " levels requested, only %" PetscInt_FMT " built + %" PetscInt_FMT " reused. Options for level(s) > %" PetscInt_FMT ", including -%spc_hpddm_coarse_ will not be taken into account\n", requested, data->N, reused,
1503: data->N, pcpre ? pcpre : ""));
1504: PetscInfo(pc, "It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%" PetscInt_FMT "_eps_threshold so that at least one local deflation vector will be selected\n", pcpre ? pcpre : "", data->N);
1505: /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */
1506: for (n = data->N - 1; n < requested - 1; ++n) {
1507: if (data->levels[n]->P) {
1508: HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE);
1509: VecDestroyVecs(1, &data->levels[n]->v[0]);
1510: VecDestroyVecs(2, &data->levels[n]->v[1]);
1511: MatDestroy(data->levels[n]->V);
1512: MatDestroy(data->levels[n]->V + 1);
1513: MatDestroy(data->levels[n]->V + 2);
1514: VecDestroy(&data->levels[n]->D);
1515: VecScatterDestroy(&data->levels[n]->scatter);
1516: }
1517: }
1518: if (reused) {
1519: for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1520: KSPDestroy(&data->levels[n]->ksp);
1521: PCDestroy(&data->levels[n]->pc);
1522: }
1523: }
1525: data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N);
1526: }
1527: /* these solvers are created after PCSetFromOptions() is called */
1528: if (pc->setfromoptionscalled) {
1529: for (n = 0; n < data->N; ++n) {
1530: if (data->levels[n]->ksp) KSPSetFromOptions(data->levels[n]->ksp);
1531: if (data->levels[n]->pc) PCSetFromOptions(data->levels[n]->pc);
1532: }
1533: pc->setfromoptionscalled = 0;
1534: }
1535: data->N += reused;
1536: if (data->share && swap) {
1537: /* swap back pointers so that variables follow the user-provided numbering */
1538: std::swap(C, data->aux);
1539: std::swap(uis, data->is);
1540: MatDestroy(&C);
1541: ISDestroy(&uis);
1542: }
1543: return 0;
1544: }
1546: /*@
1547: PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type.
1549: Collective
1551: Input Parameters:
1552: + pc - preconditioner context
1553: - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, or `PC_HPDDM_COARSE_CORRECTION_BALANCED`
1555: Options Database Key:
1556: . -pc_hpddm_coarse_correction <deflated, additive, balanced> - type of coarse correction to apply
1558: Level: intermediate
1560: .seealso: `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
1561: @*/
1562: PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type)
1563: {
1566: PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type));
1567: return 0;
1568: }
1570: /*@
1571: PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type.
1573: Input Parameter:
1574: . pc - preconditioner context
1576: Output Parameter:
1577: . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, or `PC_HPDDM_COARSE_CORRECTION_BALANCED`
1579: Level: intermediate
1581: .seealso: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
1582: @*/
1583: PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type)
1584: {
1586: if (type) {
1588: PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type));
1589: }
1590: return 0;
1591: }
1593: static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type)
1594: {
1595: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1598: data->correction = type;
1599: return 0;
1600: }
1602: static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type)
1603: {
1604: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1606: *type = data->correction;
1607: return 0;
1608: }
1610: /*@
1611: PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared.
1613: Input Parameter:
1614: . pc - preconditioner context
1616: Output Parameter:
1617: . share - whether the `KSP` is shared or not
1619: Note:
1620: 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
1621: when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`.
1623: Level: advanced
1625: .seealso: `PCHPDDM`
1626: @*/
1627: PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share)
1628: {
1630: if (share) {
1632: PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share));
1633: }
1634: return 0;
1635: }
1637: static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share)
1638: {
1639: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1641: *share = data->share;
1642: return 0;
1643: }
1645: /*@
1646: PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator.
1648: Input Parameters:
1649: + pc - preconditioner context
1650: . is - index set of the local deflation matrix
1651: - U - deflation sequential matrix stored as a `MATSEQDENSE`
1653: Level: advanced
1655: .seealso: `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()`
1656: @*/
1657: PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U)
1658: {
1662: PetscUseMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U));
1663: return 0;
1664: }
1666: static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U)
1667: {
1668: PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE);
1669: return 0;
1670: }
1672: PetscErrorCode HPDDMLoadDL_Private(PetscBool *found)
1673: {
1674: char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN];
1677: PetscStrcpy(dir, "${PETSC_LIB_DIR}");
1678: PetscOptionsGetString(NULL, NULL, "-hpddm_dir", dir, sizeof(dir), NULL);
1679: PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir);
1680: PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found);
1681: #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure since */
1682: if (!*found) { /* slepcconf.h is not yet built (and thus can't be included) */
1683: PetscStrcpy(dir, HPDDM_STR(SLEPC_LIB_DIR));
1684: PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir);
1685: PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found);
1686: }
1687: #endif
1689: PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib);
1690: return 0;
1691: }
1693: /*MC
1694: PCHPDDM - Interface with the HPDDM library.
1696: 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
1697: AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in [2021] (see references below)
1699: The matrix to be preconditioned (Pmat) may be unassembled (`MATIS`), assembled (`MATAIJ`, `MATBAIJ`, or `MATSBAIJ`), hierarchical (`MATHTOOL`), or `MATNORMAL`.
1701: For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using
1702: `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`.
1704: Options Database Keys:
1705: + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()`
1706: (not relevant with an unassembled Pmat)
1707: . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()`
1708: - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()`
1710: Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the options database prefixes
1711: .vb
1712: -pc_hpddm_levels_%d_pc_
1713: -pc_hpddm_levels_%d_ksp_
1714: -pc_hpddm_levels_%d_eps_
1715: -pc_hpddm_levels_%d_p
1716: -pc_hpddm_levels_%d_mat_type_
1717: -pc_hpddm_coarse_
1718: -pc_hpddm_coarse_p
1719: -pc_hpddm_coarse_mat_type_
1720: .ve
1722: 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
1723: -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1",
1724: aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2",
1725: and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`).
1727: 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.
1729: This preconditioner requires that you build PETSc with SLEPc (``--download-slepc``). By default, the underlying concurrent eigenproblems are solved using
1730: SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. [2011, 2013]. As stated above, SLEPc options
1731: 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.
1733: References:
1734: + 2011 - A robust two-level domain decomposition preconditioner for systems of PDEs. Spillane, Dolean, Hauret, Nataf, Pechstein, and Scheichl. Comptes Rendus Mathematique.
1735: . 2013 - Scalable domain decomposition preconditioners for heterogeneous elliptic problems. Jolivet, Hecht, Nataf, and Prud'homme. SC13.
1736: . 2015 - An introduction to domain decomposition methods: algorithms, theory, and parallel implementation. Dolean, Jolivet, and Nataf. SIAM.
1737: . 2019 - A multilevel Schwarz preconditioner based on a hierarchy of robust coarse spaces. Al Daas, Grigori, Jolivet, and Tournier. SIAM Journal on Scientific Computing.
1738: . 2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini. Computer & Mathematics with Applications.
1739: . 2022a - A robust algebraic domain decomposition preconditioner for sparse normal equations. Al Daas, Jolivet, and Scott. SIAM Journal on Scientific Computing.
1740: - 2022b - A robust algebraic multilevel domain decomposition preconditioner for sparse symmetric positive definite matrices. Al Daas and Jolivet.
1742: Level: intermediate
1744: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`,
1745: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMGetSTShareSubKSP()`,
1746: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMGetComplexities()`
1747: M*/
1748: PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc)
1749: {
1750: PC_HPDDM *data;
1751: PetscBool found;
1753: if (!loadedSym) {
1754: HPDDMLoadDL_Private(&found);
1755: if (found) PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, NULL, "PCHPDDM_Internal", (void **)&loadedSym);
1756: }
1758: PetscNew(&data);
1759: pc->data = data;
1760: pc->ops->reset = PCReset_HPDDM;
1761: pc->ops->destroy = PCDestroy_HPDDM;
1762: pc->ops->setfromoptions = PCSetFromOptions_HPDDM;
1763: pc->ops->setup = PCSetUp_HPDDM;
1764: pc->ops->apply = PCApply_HPDDM;
1765: pc->ops->matapply = PCMatApply_HPDDM;
1766: pc->ops->view = PCView_HPDDM;
1767: pc->ops->presolve = PCPreSolve_HPDDM;
1769: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM);
1770: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM);
1771: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM);
1772: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM);
1773: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM);
1774: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM);
1775: PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM);
1776: return 0;
1777: }
1779: /*@C
1780: PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`.
1782: Level: developer
1784: .seealso: `PetscInitialize()`
1785: @*/
1786: PetscErrorCode PCHPDDMInitializePackage(void)
1787: {
1788: char ename[32];
1789: PetscInt i;
1791: if (PCHPDDMPackageInitialized) return 0;
1792: PCHPDDMPackageInitialized = PETSC_TRUE;
1793: PetscRegisterFinalize(PCHPDDMFinalizePackage);
1794: /* general events registered once during package initialization */
1795: /* some of these events are not triggered in libpetsc, */
1796: /* but rather directly in libhpddm_petsc, */
1797: /* which is in charge of performing the following operations */
1799: /* domain decomposition structure from Pmat sparsity pattern */
1800: PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc);
1801: /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */
1802: PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP);
1803: /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */
1804: PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP);
1805: /* next level construction using PtAP and PtBP (not triggered in libpetsc) */
1806: PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next);
1807: static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high");
1808: for (i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) {
1809: PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i);
1810: /* events during a PCSetUp() at level #i _except_ the assembly */
1811: /* of the Galerkin operator of the coarser level #(i + 1) */
1812: PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1]);
1813: PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i);
1814: /* events during a PCApply() at level #i _except_ */
1815: /* the KSPSolve() of the coarser level #(i + 1) */
1816: PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1]);
1817: }
1818: return 0;
1819: }
1821: /*@C
1822: PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`.
1824: Level: developer
1826: .seealso: `PetscFinalize()`
1827: @*/
1828: PetscErrorCode PCHPDDMFinalizePackage(void)
1829: {
1830: PCHPDDMPackageInitialized = PETSC_FALSE;
1831: return 0;
1832: }