Actual source code: pchpddm.cxx
1: #include <petscsf.h>
2: #include <petsc/private/vecimpl.h>
3: #include <petsc/private/matimpl.h>
4: #include <petsc/private/petschpddm.h>
5: #include <petsc/private/pcimpl.h>
6: #include <petsc/private/dmimpl.h>
7: /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */
9: static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar> *const, IS, Mat, Mat, Mat, std::vector<Vec>, PC_HPDDM_Level **const) = nullptr;
11: static PetscBool PCHPDDMPackageInitialized = PETSC_FALSE;
12: static constexpr int i_0 = 0;
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", "NONE", "DEFLATED_REVERSED", "PCHPDDMCoarseCorrectionType", "PC_HPDDM_COARSE_CORRECTION_", nullptr};
22: const char *const PCHPDDMSchurPreTypes[] = {"LEAST_SQUARES", "GENEO", "PCHPDDMSchurPreType", "PC_HPDDM_SCHUR_PRE", nullptr};
24: static PetscErrorCode PCHPDDMInitializeLevels_Private(PC_HPDDM *data)
25: {
26: PetscFunctionBegin;
27: if (!data->levels) { /* usually allocated in PCSetFromOptions_HPDDM(), but PCSetUp_HPDDM() may be called without a prior PCSetFromOptions() */
28: PetscCall(PetscCalloc1(PETSC_PCHPDDM_MAXLEVELS, &data->levels));
29: PetscCall(PetscNew(data->levels));
30: data->levels[0]->parent = data;
31: data->N = 1;
32: }
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode PCReset_HPDDM(PC pc)
37: {
38: PC_HPDDM *data = (PC_HPDDM *)pc->data;
40: PetscFunctionBegin;
41: if (data->levels) {
42: for (PetscInt i = 0; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) {
43: PetscCall(KSPDestroy(&data->levels[i]->ksp));
44: PetscCall(PCDestroy(&data->levels[i]->pc));
45: PetscCall(PetscFree(data->levels[i]));
46: }
47: PetscCall(PetscFree(data->levels));
48: data->N = 0;
49: }
50: PetscCall(ISDestroy(&data->is));
51: PetscCall(MatDestroy(&data->aux));
52: PetscCall(MatDestroy(&data->B));
53: PetscCall(VecDestroy(&data->normal));
54: data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED;
55: data->Neumann = PETSC_BOOL3_UNKNOWN;
56: data->deflation = PETSC_FALSE;
57: data->setup = nullptr;
58: data->setup_ctx = nullptr;
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode PCDestroy_HPDDM(PC pc)
63: {
64: PC_HPDDM *data = (PC_HPDDM *)pc->data;
66: PetscFunctionBegin;
67: PetscCall(PCReset_HPDDM(pc));
68: PetscCall(PetscFree(data));
69: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, nullptr));
70: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", nullptr));
71: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", nullptr));
72: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", nullptr));
73: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", nullptr));
74: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", nullptr));
75: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", nullptr));
76: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", nullptr));
77: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", nullptr));
78: PetscCall(PetscObjectCompose((PetscObject)pc, "_PCHPDDM_Schur", nullptr));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static inline PetscErrorCode PCHPDDMSetAuxiliaryMat_Private(PC pc, IS is, Mat A, PetscBool deflation)
83: {
84: PC_HPDDM *data = (PC_HPDDM *)pc->data;
85: PCHPDDMCoarseCorrectionType type = data->correction;
87: PetscFunctionBegin;
89: if (is && A) {
90: PetscInt m[2];
92: PetscCall(ISGetLocalSize(is, m));
93: PetscCall(MatGetLocalSize(A, m + 1, nullptr));
94: PetscCheck(m[0] == m[1], PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Inconsistent IS and Mat sizes (%" PetscInt_FMT " v. %" PetscInt_FMT ")", m[0], m[1]);
95: }
96: if (is) {
97: PetscCall(PetscObjectReference((PetscObject)is));
98: if (data->is) { /* new overlap definition resets the PC */
99: PetscCall(PCReset_HPDDM(pc));
100: pc->setfromoptionscalled = 0;
101: pc->setupcalled = PETSC_FALSE;
102: data->correction = type;
103: }
104: PetscCall(ISDestroy(&data->is));
105: data->is = is;
106: }
107: if (A) {
108: PetscCall(PetscObjectReference((PetscObject)A));
109: PetscCall(MatDestroy(&data->aux));
110: data->aux = A;
111: }
112: data->deflation = deflation;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static inline PetscErrorCode PCHPDDMSplittingMatNormal_Private(Mat A, IS *is, Mat *splitting[])
117: {
118: Mat *sub;
119: IS zero;
121: PetscFunctionBegin;
122: PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
123: PetscCall(MatCreateSubMatrices(A, 1, is + 2, is, MAT_INITIAL_MATRIX, splitting));
124: PetscCall(MatCreateSubMatrices(**splitting, 1, is + 2, is + 1, MAT_INITIAL_MATRIX, &sub));
125: PetscCall(MatFindZeroRows(*sub, &zero));
126: PetscCall(MatDestroySubMatrices(1, &sub));
127: PetscCall(MatSetOption(**splitting, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
128: PetscCall(MatZeroRowsIS(**splitting, zero, 0.0, nullptr, nullptr));
129: PetscCall(ISDestroy(&zero));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static inline PetscErrorCode PCHPDDMSetAuxiliaryMatNormal_Private(PC pc, Mat A, Mat N, Mat *B, const char *pcpre, Vec *diagonal = nullptr, Mat B01 = nullptr)
134: {
135: PC_HPDDM *data = (PC_HPDDM *)pc->data;
136: Mat *splitting[2] = {}, aux;
137: Vec d;
138: IS is[3];
139: PetscReal norm;
140: PetscBool flg;
141: char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */
143: PetscFunctionBegin;
144: if (!B01) PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, B));
145: else PetscCall(MatTransposeMatMult(B01, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, B));
146: PetscCall(MatEliminateZeros(*B, PETSC_TRUE));
147: PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, is));
148: PetscCall(MatIncreaseOverlap(*B, 1, is, 1));
149: PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, is + 2));
150: PetscCall(ISEmbed(is[0], is[2], PETSC_TRUE, is + 1));
151: PetscCall(ISDestroy(is + 2));
152: PetscCall(ISCreateStride(PETSC_COMM_SELF, A->rmap->N, 0, 1, is + 2));
153: PetscCall(PCHPDDMSplittingMatNormal_Private(A, is, &splitting[0]));
154: if (B01) {
155: PetscCall(PCHPDDMSplittingMatNormal_Private(B01, is, &splitting[1]));
156: PetscCall(MatDestroy(&B01));
157: }
158: PetscCall(ISDestroy(is + 2));
159: PetscCall(ISDestroy(is + 1));
160: PetscCall(PetscOptionsGetString(((PetscObject)pc)->options, pcpre, "-pc_hpddm_levels_1_sub_pc_type", type, sizeof(type), nullptr));
161: PetscCall(PetscStrcmp(type, PCQR, &flg));
162: if (!flg) {
163: Mat conjugate = *splitting[splitting[1] ? 1 : 0];
165: if (PetscDefined(USE_COMPLEX) && !splitting[1]) {
166: PetscCall(MatDuplicate(*splitting[0], MAT_COPY_VALUES, &conjugate));
167: PetscCall(MatConjugate(conjugate));
168: }
169: PetscCall(MatTransposeMatMult(conjugate, *splitting[0], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &aux));
170: if (PetscDefined(USE_COMPLEX) && !splitting[1]) PetscCall(MatDestroy(&conjugate));
171: else if (splitting[1]) PetscCall(MatDestroySubMatrices(1, &splitting[1]));
172: PetscCall(MatNorm(aux, NORM_FROBENIUS, &norm));
173: PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
174: if (diagonal) {
175: PetscReal norm;
177: PetscCall(VecScale(*diagonal, -1.0));
178: PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm));
179: if (norm > PETSC_SMALL) {
180: PetscSF scatter;
181: PetscInt n;
183: PetscCall(ISGetLocalSize(*is, &n));
184: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), n, PETSC_DECIDE, &d));
185: PetscCall(VecScatterCreate(*diagonal, *is, d, nullptr, &scatter));
186: PetscCall(VecScatterBegin(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD));
187: PetscCall(VecScatterEnd(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD));
188: PetscCall(PetscSFDestroy(&scatter));
189: PetscCall(MatDiagonalSet(aux, d, ADD_VALUES));
190: PetscCall(VecDestroy(&d));
191: } else PetscCall(VecDestroy(diagonal));
192: }
193: if (!diagonal) PetscCall(MatShift(aux, PETSC_SMALL * norm));
194: PetscCall(MatEliminateZeros(aux, PETSC_TRUE));
195: } else {
196: PetscBool flg;
198: PetscCheck(!splitting[1], PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use PCQR when A01 != A10^T");
199: if (diagonal) {
200: PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm));
201: PetscCheck(norm < PETSC_SMALL, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Nonzero diagonal A11 block");
202: PetscCall(VecDestroy(diagonal));
203: }
204: PetscCall(PetscObjectTypeCompare((PetscObject)N, MATNORMAL, &flg));
205: if (flg) PetscCall(MatCreateNormal(*splitting[0], &aux));
206: else PetscCall(MatCreateNormalHermitian(*splitting[0], &aux));
207: }
208: PetscCall(MatDestroySubMatrices(1, &splitting[0]));
209: PetscCall(PCHPDDMSetAuxiliaryMat(pc, *is, aux, nullptr, nullptr));
210: data->Neumann = PETSC_BOOL3_TRUE;
211: PetscCall(ISDestroy(is));
212: PetscCall(MatDestroy(&aux));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx)
217: {
218: PC_HPDDM *data = (PC_HPDDM *)pc->data;
220: PetscFunctionBegin;
221: PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, A, PETSC_FALSE));
222: if (setup) {
223: data->setup = setup;
224: data->setup_ctx = setup_ctx;
225: }
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level.
232: Input Parameters:
233: + pc - preconditioner context
234: . is - index set of the local auxiliary, e.g., Neumann, matrix
235: . A - auxiliary sequential matrix
236: . setup - function for generating the auxiliary matrix entries, may be `NULL`
237: - ctx - context for `setup`, may be `NULL`
239: Calling sequence of `setup`:
240: + J - matrix whose values are to be set
241: . t - time
242: . X - linearization point
243: . X_t - time-derivative of the linearization point
244: . s - step
245: . ovl - index set of the local auxiliary, e.g., Neumann, matrix
246: - ctx - context for `setup`, may be `NULL`
248: Level: intermediate
250: Note:
251: As an example, in a finite element context with nonoverlapping subdomains plus (overlapping) ghost elements, this could be the unassembled (Neumann)
252: local overlapping operator. As opposed to the assembled (Dirichlet) local overlapping operator obtained by summing neighborhood contributions
253: at the interface of ghost elements.
255: Fortran Notes:
256: Only `PETSC_NULL_FUNCTION` is supported for `setup` and `ctx` is never accessed
258: .seealso: [](ch_ksp), `PCHPDDM`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetRHSMat()`, `MATIS`
259: @*/
260: PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, PetscCtx ctx), PetscCtx ctx)
261: {
262: PetscFunctionBegin;
266: PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode (*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *), (pc, is, A, setup, ctx));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has)
271: {
272: PC_HPDDM *data = (PC_HPDDM *)pc->data;
274: PetscFunctionBegin;
275: data->Neumann = PetscBoolToBool3(has);
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: PCHPDDMHasNeumannMat - Informs `PCHPDDM` that the `Mat` passed to `PCHPDDMSetAuxiliaryMat()` is the local Neumann matrix.
282: Input Parameters:
283: + pc - preconditioner context
284: - has - Boolean value
286: Level: intermediate
288: Notes:
289: This may be used to bypass a call to `MatCreateSubMatrices()` and to `MatConvert()` for `MATSBAIJ` matrices.
291: If a function is composed with DMCreateNeumannOverlap_C 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`.
293: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()`
294: @*/
295: PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has)
296: {
297: PetscFunctionBegin;
299: PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B)
304: {
305: PC_HPDDM *data = (PC_HPDDM *)pc->data;
307: PetscFunctionBegin;
308: PetscCall(PetscObjectReference((PetscObject)B));
309: PetscCall(MatDestroy(&data->B));
310: data->B = B;
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@
315: PCHPDDMSetRHSMat - Sets the right-hand side matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level.
317: Input Parameters:
318: + pc - preconditioner context
319: - B - right-hand side sequential matrix
321: Level: advanced
323: Note:
324: Must be used in conjunction with `PCHPDDMSetAuxiliaryMat`(N), so that Nv = lambda Bv is solved using `EPSSetOperators`(N, B).
325: 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.
327: .seealso: [](ch_ksp), `PCHPDDMSetAuxiliaryMat()`, `PCHPDDM`
328: @*/
329: PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B)
330: {
331: PetscFunctionBegin;
333: if (B) {
335: PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B));
336: }
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode PCSetFromOptions_HPDDM(PC pc, PetscOptionItems PetscOptionsObject)
341: {
342: PC_HPDDM *data = (PC_HPDDM *)pc->data;
343: char prefix[256], deprecated[256];
344: int i = 1;
345: PetscMPIInt size, previous;
346: PetscInt n, overlap = 1;
347: PCHPDDMCoarseCorrectionType type;
348: PetscBool flg = PETSC_TRUE, set;
350: PetscFunctionBegin;
351: PetscCall(PCHPDDMInitializeLevels_Private(data));
352: PetscOptionsHeadBegin(PetscOptionsObject, "PCHPDDM options");
353: PetscCall(PetscOptionsBoundedInt("-pc_hpddm_harmonic_overlap", "Overlap prior to computing local harmonic extensions", "PCHPDDM", overlap, &overlap, &set, 1));
354: if (!set) overlap = -1;
355: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
356: previous = size;
357: while (i < PETSC_PCHPDDM_MAXLEVELS) {
358: PetscInt p = 1;
360: if (!data->levels[i - 1]) PetscCall(PetscNew(data->levels + i - 1));
361: data->levels[i - 1]->parent = data;
362: /* if the previous level has a single process, it is not possible to coarsen further */
363: if (previous == 1 || !flg) break;
364: data->levels[i - 1]->nu = 0;
365: data->levels[i - 1]->threshold = -1.0;
366: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i));
367: PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "EPSSetDimensions", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, nullptr, 0));
368: PetscCall(PetscSNPrintf(deprecated, sizeof(deprecated), "-pc_hpddm_levels_%d_eps_threshold", i));
369: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold_absolute", i));
370: PetscCall(PetscOptionsDeprecated(deprecated, prefix, "3.24", nullptr));
371: PetscCall(PetscOptionsReal(prefix, "Local absolute threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, nullptr));
372: if (i == 1) {
373: PetscCheck(overlap == -1 || PetscAbsReal(data->levels[i - 1]->threshold + static_cast<PetscReal>(1.0)) < PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply both -pc_hpddm_levels_1_eps_threshold_absolute and -pc_hpddm_harmonic_overlap");
374: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_nsv", i));
375: if (overlap != -1) {
376: PetscInt nsv = 0;
377: PetscBool set[2] = {PETSC_FALSE, PETSC_FALSE};
379: PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "SVDSetDimensions", nsv, &nsv, nullptr, 0));
380: PetscCheck(data->levels[0]->nu == 0 || nsv == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply both -pc_hpddm_levels_1_eps_nev and -pc_hpddm_levels_1_svd_nsv");
381: if (data->levels[0]->nu == 0) { /* -eps_nev has not been used, so nu is 0 */
382: data->levels[0]->nu = nsv; /* nu may still be 0 if -svd_nsv has not been used */
383: PetscCall(PetscSNPrintf(deprecated, sizeof(deprecated), "-pc_hpddm_levels_%d_svd_relative_threshold", i));
384: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_threshold_relative", i));
385: PetscCall(PetscOptionsDeprecated(deprecated, prefix, "3.24", nullptr));
386: PetscCall(PetscOptionsReal(prefix, "Local relative threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[0]->threshold, &data->levels[0]->threshold, set)); /* cache whether this option has been used or not to error out in case of exclusive options being used simultaneously later on */
387: }
388: if (data->levels[0]->nu == 0 || nsv == 0) { /* if neither -eps_nev nor -svd_nsv has been used */
389: PetscCall(PetscSNPrintf(deprecated, sizeof(deprecated), "-pc_hpddm_levels_%d_eps_relative_threshold", i));
390: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold_relative", i));
391: PetscCall(PetscOptionsDeprecated(deprecated, prefix, "3.24", nullptr));
392: PetscCall(PetscOptionsReal(prefix, "Local relative threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[0]->threshold, &data->levels[0]->threshold, set + 1));
393: PetscCheck(!set[0] || !set[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply both -pc_hpddm_levels_1_eps_threshold_relative and -pc_hpddm_levels_1_svd_threshold_relative");
394: }
395: PetscCheck(data->levels[0]->nu || PetscAbsReal(data->levels[i - 1]->threshold + static_cast<PetscReal>(1.0)) > PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need to supply at least one of 1) -pc_hpddm_levels_1_eps_nev, 2) -pc_hpddm_levels_1_svd_nsv, 3) -pc_hpddm_levels_1_eps_threshold_relative, 4) -pc_hpddm_levels_1_svd_threshold_relative (for nonsymmetric matrices, only option 2 and option 4 are appropriate)");
396: } else if (PetscDefined(USE_DEBUG)) {
397: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg));
398: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply -%spc_hpddm_levels_%d_svd_nsv without -%spc_hpddm_harmonic_overlap", PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "", i,
399: PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "");
400: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_threshold_relative", i));
401: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg));
402: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply -%spc_hpddm_levels_%d_svd_threshold_relative without -%spc_hpddm_harmonic_overlap", PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "", i,
403: PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "");
404: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold_relative", i));
405: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg));
406: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply -%spc_hpddm_levels_%d_eps_threshold_relative without -%spc_hpddm_harmonic_overlap, maybe you meant to use -%spc_hpddm_levels_%d_eps_threshold_absolute?",
407: PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "", i, PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "", PetscOptionsObject->prefix ? PetscOptionsObject->prefix : "", i);
408: }
409: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp"));
410: PetscCall(PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "PCHPDDMSetSTShareSubKSP", PETSC_FALSE, &data->share, nullptr));
411: }
412: /* if there is no prescribed coarsening, just break out of the loop */
413: if (data->levels[i - 1]->threshold <= PetscReal() && data->levels[i - 1]->nu <= 0 && !(data->deflation && i == 1)) break;
414: else {
415: ++i;
416: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i));
417: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg));
418: if (!flg) {
419: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold_absolute", i));
420: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg));
421: }
422: if (flg) {
423: /* if there are coarsening options for the next level, then register it */
424: /* otherwise, don't to avoid having both options levels_N_p and coarse_p */
425: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i));
426: PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "PCHPDDM", p, &p, &flg, 1, PetscMax(1, previous / 2)));
427: previous = p;
428: }
429: }
430: }
431: data->N = i;
432: n = 1;
433: if (i > 1) {
434: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p"));
435: PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "PCHPDDM", n, &n, nullptr, 1, PetscMax(1, previous / 2)));
436: #if PetscDefined(HAVE_MUMPS)
437: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "pc_hpddm_coarse_"));
438: PetscCall(PetscOptionsHasName(PetscOptionsObject->options, prefix, "-mat_mumps_use_omp_threads", &flg));
439: if (flg) {
440: char type[64]; /* same size as in src/ksp/pc/impls/factor/factimpl.c */
442: PetscCall(PetscStrncpy(type, n > 1 && PetscDefined(HAVE_MUMPS) ? MATSOLVERMUMPS : MATSOLVERPETSC, sizeof(type))); /* default solver for a MatMPIAIJ or a MatSeqAIJ */
443: PetscCall(PetscOptionsGetString(PetscOptionsObject->options, prefix, "-pc_factor_mat_solver_type", type, sizeof(type), nullptr));
444: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
445: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "-%smat_mumps_use_omp_threads and -%spc_factor_mat_solver_type != %s", prefix, prefix, MATSOLVERMUMPS);
446: size = n;
447: n = -1;
448: PetscCall(PetscOptionsGetInt(PetscOptionsObject->options, prefix, "-mat_mumps_use_omp_threads", &n, nullptr));
449: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need to specify a positive integer for -%smat_mumps_use_omp_threads", prefix);
450: PetscCheck(n * size <= previous, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%d MPI process%s x %d OpenMP thread%s greater than %d available MPI process%s for the coarsest operator", (int)size, size > 1 ? "es" : "", (int)n, n > 1 ? "s" : "", (int)previous, previous > 1 ? "es" : "");
451: }
452: #endif
453: PetscCall(PetscOptionsEnum("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, (PetscEnum)data->correction, (PetscEnum *)&type, &flg));
454: if (flg) PetscCall(PCHPDDMSetCoarseCorrectionType(pc, type));
455: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann"));
456: PetscCall(PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", PetscBool3ToBool(data->Neumann), &flg, &set));
457: if (set) data->Neumann = PetscBoolToBool3(flg);
458: data->log_separate = PETSC_FALSE;
459: if (PetscDefined(USE_LOG)) {
460: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate"));
461: PetscCall(PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", nullptr, data->log_separate, &data->log_separate, nullptr));
462: }
463: }
464: PetscOptionsHeadEnd();
465: for (; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) {
466: PetscCall(KSPDestroy(&data->levels[i]->ksp));
467: PetscCall(PCDestroy(&data->levels[i]->pc));
468: PetscCall(PetscFree(data->levels[i]));
469: }
470: if (data->levels[0]->ksp) { /* PCSetUp_HPDDM() may have created this KSP initially as a single-level solver before PCSetFromOptions() enabled multiple levels */
471: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", data->N > 1 ? "levels_1" : "coarse"));
472: PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix));
473: }
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: template <bool transpose>
478: static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y)
479: {
480: PC_HPDDM *data = (PC_HPDDM *)pc->data;
482: PetscFunctionBegin;
483: PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
484: PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM");
485: if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); /* coarser-level events are directly triggered in HPDDM */
486: if (!transpose) PetscCall(KSPSolve(data->levels[0]->ksp, x, y));
487: else PetscCall(KSPSolveTranspose(data->levels[0]->ksp, x, y));
488: PetscCall(KSPCheckSolve(data->levels[0]->ksp, pc, y));
489: if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: template <bool transpose>
494: static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y)
495: {
496: PC_HPDDM *data = (PC_HPDDM *)pc->data;
498: PetscFunctionBegin;
499: PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
500: PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM");
501: if (!transpose) PetscCall(KSPMatSolve(data->levels[0]->ksp, X, Y));
502: else PetscCall(KSPMatSolveTranspose(data->levels[0]->ksp, X, Y));
503: PetscCall(KSPCheckMatSolve(data->levels[0]->ksp, pc, Y));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@
508: PCHPDDMGetComplexities - Computes the grid and operator complexities.
510: Collective
512: Input Parameter:
513: . pc - preconditioner context
515: Output Parameters:
516: + gc - grid complexity $ \sum_i m_i / m_1 $
517: - oc - operator complexity $ \sum_i nnz_i / nnz_1 $
519: Level: advanced
521: .seealso: [](ch_ksp), `PCMGGetGridComplexity()`, `PCHPDDM`, `PCHYPRE`, `PCGAMG`
522: @*/
523: PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc)
524: {
525: PC_HPDDM *data = (PC_HPDDM *)pc->data;
526: MatInfo info;
527: PetscLogDouble accumulate[2]{}, nnz1 = 1.0, m1 = 1.0;
529: PetscFunctionBegin;
530: if (gc) {
531: PetscAssertPointer(gc, 2);
532: *gc = 0;
533: }
534: if (oc) {
535: PetscAssertPointer(oc, 3);
536: *oc = 0;
537: }
538: for (PetscInt n = 0; n < data->N; ++n) {
539: if (data->levels[n]->ksp) {
540: Mat P, A = nullptr;
541: PetscInt m;
542: PetscBool flg = PETSC_FALSE;
544: PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &P));
545: PetscCall(MatGetSize(P, &m, nullptr));
546: accumulate[0] += m;
547: if (n == 0) {
548: PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, ""));
549: if (flg) {
550: PetscCall(MatConvert(P, MATAIJ, MAT_INITIAL_MATRIX, &A));
551: P = A;
552: } else {
553: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg));
554: PetscCall(PetscObjectReference((PetscObject)P));
555: }
556: }
557: if (!A && flg) accumulate[1] += m * m; /* assumption that a MATSCHURCOMPLEMENT is dense if stored explicitly */
558: else if (P->ops->getinfo) {
559: PetscCall(MatGetInfo(P, MAT_GLOBAL_SUM, &info));
560: accumulate[1] += info.nz_used;
561: }
562: if (n == 0) {
563: m1 = m;
564: if (!A && flg) nnz1 = m * m;
565: else if (P->ops->getinfo) nnz1 = info.nz_used;
566: PetscCall(MatDestroy(&P));
567: }
568: }
569: }
570: /* only process #0 has access to the full hierarchy by construction, so broadcast to ensure consistent outputs */
571: PetscCallMPI(MPI_Bcast(accumulate, 2, MPIU_PETSCLOGDOUBLE, 0, PetscObjectComm((PetscObject)pc)));
572: if (gc) *gc = static_cast<PetscReal>(accumulate[0] / m1);
573: if (oc) *oc = static_cast<PetscReal>(accumulate[1] / nnz1);
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer)
578: {
579: PC_HPDDM *data = (PC_HPDDM *)pc->data;
580: PetscViewer subviewer;
581: PetscViewerFormat format;
582: PetscSubcomm subcomm;
583: PetscReal oc, gc;
584: PetscInt tabs;
585: PetscMPIInt size, color, rank;
586: PetscBool flg;
587: const char *name;
589: PetscFunctionBegin;
590: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg));
591: if (flg) {
592: PetscCall(PetscViewerASCIIPrintf(viewer, "level%s: %" PetscInt_FMT "\n", data->N > 1 ? "s" : "", data->N));
593: PetscCall(PCHPDDMGetComplexities(pc, &gc, &oc));
594: if (data->N > 1) {
595: if (!data->deflation) {
596: PetscCall(PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[PetscBool3ToBool(data->Neumann)]));
597: PetscCall(PetscViewerASCIIPrintf(viewer, "shared subdomain KSP between SLEPc and PETSc? %s\n", PetscBools[data->share]));
598: } else PetscCall(PetscViewerASCIIPrintf(viewer, "user-supplied deflation matrix\n"));
599: PetscCall(PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction]));
600: PetscCall(PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : ""));
601: PetscCall(PetscViewerASCIIGetTab(viewer, &tabs));
602: PetscCall(PetscViewerASCIISetTab(viewer, 0));
603: for (PetscInt i = 1; i < data->N; ++i) {
604: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, data->levels[i - 1]->nu));
605: if (data->levels[i - 1]->threshold > static_cast<PetscReal>(-0.1)) PetscCall(PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold));
606: }
607: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
608: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
609: }
610: PetscCall(PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc));
611: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
612: if (data->levels && data->levels[0]->ksp) {
613: PetscCall(KSPView(data->levels[0]->ksp, viewer));
614: if (data->levels[0]->pc) PetscCall(PCView(data->levels[0]->pc, viewer));
615: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
616: for (PetscInt i = 1; i < data->N; ++i) {
617: if (data->levels[i]->ksp) color = 1;
618: else color = 0;
619: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm));
620: PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2)));
621: PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
622: PetscCall(PetscViewerASCIIPushTab(viewer));
623: PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
624: if (color == 1) {
625: PetscCall(KSPView(data->levels[i]->ksp, subviewer));
626: if (data->levels[i]->pc) PetscCall(PCView(data->levels[i]->pc, subviewer));
627: PetscCall(PetscViewerFlush(subviewer));
628: }
629: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
630: PetscCall(PetscViewerASCIIPopTab(viewer));
631: PetscCall(PetscSubcommDestroy(&subcomm));
632: }
633: }
634: PetscCall(PetscViewerGetFormat(viewer, &format));
635: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
636: PetscCall(PetscViewerFileGetName(viewer, &name));
637: if (name) {
638: Mat aux[2];
639: IS is;
640: const PetscInt *indices;
641: PetscInt m, n, sizes[5] = {pc->mat->rmap->n, pc->mat->cmap->n, pc->mat->rmap->N, pc->mat->cmap->N, 0};
642: char *tmp;
643: std::string prefix, suffix;
644: size_t pos;
646: PetscCall(PetscStrstr(name, ".", &tmp));
647: if (tmp) {
648: pos = std::distance(const_cast<char *>(name), tmp);
649: prefix = std::string(name, pos);
650: suffix = std::string(name + pos + 1);
651: } else prefix = name;
652: if (data->aux) {
653: PetscCall(MatGetSize(data->aux, &m, &n));
654: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), aux));
655: PetscCall(MatSetSizes(aux[0], m, n, PETSC_DETERMINE, PETSC_DETERMINE));
656: PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg));
657: if (flg) PetscCall(MatSetType(aux[0], MATMPIAIJ));
658: else {
659: PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQBAIJ, &flg));
660: if (flg) PetscCall(MatSetType(aux[0], MATMPIBAIJ));
661: else {
662: PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQSBAIJ, &flg));
663: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "MatType of auxiliary Mat (%s) is not any of the following: MATSEQAIJ, MATSEQBAIJ, or MATSEQSBAIJ", ((PetscObject)data->aux)->type_name);
664: PetscCall(MatSetType(aux[0], MATMPISBAIJ));
665: }
666: }
667: PetscCall(MatSetBlockSizesFromMats(aux[0], data->aux, data->aux));
668: PetscCall(MatAssemblyBegin(aux[0], MAT_FINAL_ASSEMBLY));
669: PetscCall(MatAssemblyEnd(aux[0], MAT_FINAL_ASSEMBLY));
670: PetscCall(MatGetDiagonalBlock(aux[0], aux + 1));
671: PetscCall(MatCopy(data->aux, aux[1], DIFFERENT_NONZERO_PATTERN));
672: PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_aux_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer));
673: PetscCall(MatView(aux[0], subviewer));
674: PetscCall(PetscViewerDestroy(&subviewer));
675: PetscCall(MatDestroy(aux));
676: }
677: if (data->is) {
678: PetscCall(ISGetIndices(data->is, &indices));
679: PetscCall(ISGetSize(data->is, sizes + 4));
680: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), sizes[4], indices, PETSC_USE_POINTER, &is));
681: PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_is_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer));
682: PetscCall(ISView(is, subviewer));
683: PetscCall(PetscViewerDestroy(&subviewer));
684: PetscCall(ISDestroy(&is));
685: PetscCall(ISRestoreIndices(data->is, &indices));
686: }
687: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), PETSC_STATIC_ARRAY_LENGTH(sizes), sizes, PETSC_USE_POINTER, &is));
688: PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_sizes_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer));
689: PetscCall(ISView(is, subviewer));
690: PetscCall(PetscViewerDestroy(&subviewer));
691: PetscCall(ISDestroy(&is));
692: }
693: }
694: }
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec)
699: {
700: PC_HPDDM *data = (PC_HPDDM *)pc->data;
701: Mat A;
702: PetscBool flg;
704: PetscFunctionBegin;
705: if (ksp) {
706: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg));
707: if (flg && !data->normal) {
708: PetscCall(KSPGetOperators(ksp, &A, nullptr));
709: PetscCall(MatCreateVecs(A, nullptr, &data->normal)); /* temporary Vec used in PCApply_HPDDMShell() for coarse grid corrections */
710: } else if (!flg) {
711: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &flg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPECGRR, KSPPIPELCG, KSPPIPEPRCG, KSPPIPECG2, KSPSTCG, KSPFCG, KSPPIPEFCG, KSPMINRES, KSPNASH, KSPSYMMLQ, ""));
712: if (!flg) {
713: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPHPDDM, &flg));
714: if (flg) {
715: KSPHPDDMType type;
717: PetscCall(KSPHPDDMGetType(ksp, &type));
718: flg = (type == KSP_HPDDM_TYPE_CG || type == KSP_HPDDM_TYPE_BCG || type == KSP_HPDDM_TYPE_BFBCG ? PETSC_TRUE : PETSC_FALSE);
719: }
720: }
721: }
722: if (flg) {
723: if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
724: PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_hpddm_coarse_correction", &flg));
725: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "PCHPDDMCoarseCorrectionType %s is known to be not symmetric, but KSPType %s requires a symmetric PC, if you insist on using this configuration, use the additional option -%spc_hpddm_coarse_correction %s, or alternatively, switch to a symmetric PCHPDDMCoarseCorrectionType such as %s",
726: PCHPDDMCoarseCorrectionTypes[data->correction], ((PetscObject)ksp)->type_name, ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", PCHPDDMCoarseCorrectionTypes[data->correction], PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_BALANCED]);
727: }
728: for (PetscInt n = 0; n < data->N; ++n) {
729: if (data->levels[n]->pc) {
730: PetscCall(PetscObjectTypeCompare((PetscObject)data->levels[n]->pc, PCASM, &flg));
731: if (flg) {
732: PCASMType type;
734: PetscCall(PCASMGetType(data->levels[n]->pc, &type));
735: if (type == PC_ASM_RESTRICT || type == PC_ASM_INTERPOLATE) {
736: PetscCall(PetscOptionsHasName(((PetscObject)data->levels[n]->pc)->options, ((PetscObject)data->levels[n]->pc)->prefix, "-pc_asm_type", &flg));
737: PetscCheck(flg, PetscObjectComm((PetscObject)data->levels[n]->pc), PETSC_ERR_ARG_INCOMP, "PCASMType %s is known to be not symmetric, but KSPType %s requires a symmetric PC, if you insist on using this configuration, use the additional option -%spc_asm_type %s, or alternatively, switch to a symmetric PCASMType such as %s", PCASMTypes[type],
738: ((PetscObject)ksp)->type_name, ((PetscObject)data->levels[n]->pc)->prefix, PCASMTypes[type], PCASMTypes[PC_ASM_BASIC]);
739: }
740: }
741: }
742: }
743: }
744: }
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: static PetscErrorCode PCSetUp_HPDDMShell(PC pc)
749: {
750: PC_HPDDM_Level *ctx;
751: Mat A, P;
752: Vec x;
753: const char *pcpre;
755: PetscFunctionBegin;
756: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
757: PetscCall(KSPGetOptionsPrefix(ctx->ksp, &pcpre));
758: PetscCall(KSPGetOperators(ctx->ksp, &A, &P));
759: /* smoother */
760: PetscCall(PCSetOptionsPrefix(ctx->pc, pcpre));
761: PetscCall(PCSetOperators(ctx->pc, A, P));
762: if (!ctx->v[0]) {
763: PetscCall(VecDuplicateVecs(ctx->D, 1, &ctx->v[0]));
764: if (!std::is_same<PetscScalar, PetscReal>::value) PetscCall(VecDestroy(&ctx->D));
765: PetscCall(MatCreateVecs(A, &x, nullptr));
766: PetscCall(VecDuplicateVecs(x, 2, &ctx->v[1]));
767: PetscCall(VecDestroy(&x));
768: }
769: std::fill_n(ctx->V, 3, nullptr);
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: template <bool transpose = false, class Type = Vec, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
774: static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y)
775: {
776: PC_HPDDM_Level *ctx;
778: PetscFunctionBegin;
779: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
780: /* going from PETSc to HPDDM numbering */
781: PetscCall(VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD));
782: PetscCall(VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD));
783: PetscCall(ctx->P->deflation<false, transpose>(ctx->v[0][0], ctx->D)); /* y = Q x */
784: /* going from HPDDM to PETSc numbering */
785: PetscCall(VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE));
786: PetscCall(VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: template <bool transpose = false, class Type = Mat, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
791: static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y)
792: {
793: PC_HPDDM_Level *ctx;
794: PetscInt N, ld[2];
796: PetscFunctionBegin;
797: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
798: PetscCall(MatGetSize(X, nullptr, &N));
799: PetscCall(MatDenseGetLDA(X, ld));
800: PetscCall(MatDenseGetLDA(Y, ld + 1));
801: PetscCheck(ld[0] == ld[1], PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Leading dimension of input Mat different than the one of output Mat");
802: /* going from PETSc to HPDDM numbering */
803: PetscCall(MatDenseScatter_Private(ctx->scatter, X, ctx->V[0], INSERT_VALUES, SCATTER_FORWARD));
804: PetscCall(ctx->P->deflation<false, transpose>(ctx->V[0], ctx->D)); /* Y = Q X */
805: /* going from HPDDM to PETSc numbering */
806: PetscCall(MatDenseScatter_Private(ctx->scatter, ctx->V[0], Y, INSERT_VALUES, SCATTER_REVERSE));
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: static PetscErrorCode PCApply_HPDDMShell(PC pc, Vec x, Vec y)
811: {
812: PC_HPDDM_Level *ctx;
813: Mat A;
815: PetscFunctionBegin;
816: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
817: PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
818: PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr));
819: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCApply(ctx->pc, x, y)); /* y = M^-1 x */
820: else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
821: PetscCall(PCApply(ctx->pc, x, y)); /* y = M^-1 x */
822: PetscCall(MatMult(A, y, ctx->v[1][0]));
823: PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* z = (I - A M^-1) x */
824: PetscCall(PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][0])); /* z = Q (I - A M^-1) x */
825: PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-1 x + Q (I - A M^-1) x */
826: } else {
827: PetscCall(PCHPDDMDeflate_Private(pc, x, y)); /* y = Q x */
828: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
829: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMult(A, y, ctx->v[1][0])); /* y = A Q x */
830: else {
831: /* KSPLSQR and finest level */
832: PetscCall(MatMult(A, y, ctx->parent->normal)); /* y = A Q x */
833: PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][0])); /* y = A^T A Q x */
834: }
835: PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* y = (I - A Q) x */
836: PetscCall(PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0])); /* y = M^-1 (I - A Q) x */
837: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
838: if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1])); /* z = A^T y */
839: else {
840: PetscCall(MatMult(A, ctx->v[1][0], ctx->parent->normal));
841: PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][1])); /* z = A^T A y */
842: }
843: PetscCall(PCHPDDMDeflate_Private<true>(pc, ctx->v[1][1], ctx->v[1][1])); /* z = Q^T z */
844: PetscCall(VecAXPBYPCZ(y, -1.0, 1.0, 1.0, ctx->v[1][1], ctx->v[1][0])); /* y = (I - Q^T A^T) y + Q x */
845: } else PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = Q M^-1 (I - A Q) x + Q x */
846: } else {
847: PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
848: PetscCall(PCApply(ctx->pc, x, ctx->v[1][0]));
849: PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-1 x + Q x */
850: }
851: }
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: template <bool transpose>
856: static PetscErrorCode PCHPDDMMatApply_Private(PC_HPDDM_Level *ctx, Mat Y, PetscBool *reset)
857: {
858: Mat A, *ptr;
859: PetscScalar *array;
860: PetscInt m, M, N, prev = 0;
861: PetscContainer container = nullptr;
863: PetscFunctionBegin;
864: PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr));
865: PetscCall(MatGetSize(Y, nullptr, &N));
866: PetscCall(PetscObjectQuery((PetscObject)A, "_HPDDM_MatProduct", (PetscObject *)&container));
867: if (container) { /* MatProduct container already attached */
868: PetscCall(PetscContainerGetPointer(container, static_cast<void *>(&ptr)));
869: if (ptr[1] != ctx->V[2]) /* Mat has changed or may have been set first in KSPHPDDM */
870: for (m = 0; m < 2; ++m) {
871: PetscCall(MatDestroy(ctx->V + m + 1));
872: ctx->V[m + 1] = ptr[m];
873: PetscCall(PetscObjectReference((PetscObject)ctx->V[m + 1]));
874: }
875: }
876: if (ctx->V[1]) PetscCall(MatGetSize(ctx->V[1], nullptr, &prev));
877: if (N != prev || !ctx->V[0]) {
878: PetscCall(MatDestroy(ctx->V));
879: PetscCall(VecGetLocalSize(ctx->v[0][0], &m));
880: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)Y), A->defaultvectype, m, PETSC_DECIDE, PETSC_DECIDE, N, PETSC_DECIDE, nullptr, ctx->V));
881: if (N != prev) {
882: PetscMemType mtype;
884: PetscCall(MatDestroy(ctx->V + 1));
885: PetscCall(MatDestroy(ctx->V + 2));
886: PetscCall(MatGetLocalSize(Y, &m, nullptr));
887: PetscCall(MatGetSize(Y, &M, nullptr));
888: PetscCall(MatDenseGetArrayWriteAndMemType(ctx->V[0], &array, &mtype));
889: PetscCall(MatCreateDenseWithMemType(PetscObjectComm((PetscObject)Y), mtype, m, PETSC_DECIDE, M, N, PETSC_DECIDE, ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED ? array : nullptr, ctx->V + 1));
890: PetscCall(MatDenseRestoreArrayWriteAndMemType(ctx->V[0], &array));
891: PetscCall(MatDuplicate(ctx->V[1], MAT_DO_NOT_COPY_VALUES, ctx->V + 2));
892: PetscCall(MatProductCreateWithMat(A, !transpose ? Y : ctx->V[2], nullptr, ctx->V[1]));
893: PetscCall(MatProductSetType(ctx->V[1], !transpose ? MATPRODUCT_AB : MATPRODUCT_AtB));
894: PetscCall(MatProductSetFromOptions(ctx->V[1]));
895: PetscCall(MatProductSymbolic(ctx->V[1]));
896: if (!container)
897: PetscCall(PetscObjectContainerCompose((PetscObject)A, "_HPDDM_MatProduct", static_cast<void *>(ctx->V + 1), nullptr)); /* no MatProduct container attached, create one to be queried in KSPHPDDM or at the next call to PCMatApply() */
898: else PetscCall(PetscContainerSetPointer(container, static_cast<void *>(ctx->V + 1))); /* need to compose B and D from MatProductCreateWithMat(A, B, NULL, D), which are stored in the contiguous array ctx->V */
899: }
900: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
901: PetscCall(MatProductCreateWithMat(A, !transpose ? ctx->V[1] : Y, nullptr, ctx->V[2]));
902: PetscCall(MatProductSetType(ctx->V[2], !transpose ? MATPRODUCT_AtB : MATPRODUCT_AB));
903: PetscCall(MatProductSetFromOptions(ctx->V[2]));
904: PetscCall(MatProductSymbolic(ctx->V[2]));
905: }
906: PetscCallCXX(ctx->P->start(N));
907: }
908: if (N == prev || container) { /* when MatProduct container is attached, always need to MatProductReplaceMats() since KSPHPDDM may have replaced the Mat as well */
909: PetscCall(MatProductReplaceMats(nullptr, !transpose ? Y : ctx->V[2], nullptr, ctx->V[1]));
910: if (container && ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) {
911: PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array));
912: PetscCall(MatDensePlaceArray(ctx->V[1], array));
913: PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array));
914: *reset = PETSC_TRUE;
915: }
916: }
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*
921: PCMatApply_HPDDMShell - Variant of PCApply_HPDDMShell() for blocks of vectors.
923: Input Parameters:
924: + pc - preconditioner context
925: - X - block of input vectors
927: Output Parameter:
928: . Y - block of output vectors
930: Level: advanced
932: .seealso: [](ch_ksp), `PCHPDDM`, `PCApply_HPDDMShell()`, `PCHPDDMCoarseCorrectionType`
933: */
934: static PetscErrorCode PCMatApply_HPDDMShell(PC pc, Mat X, Mat Y)
935: {
936: PC_HPDDM_Level *ctx;
937: PetscBool reset = PETSC_FALSE;
939: PetscFunctionBegin;
940: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
941: PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
942: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCMatApply(ctx->pc, X, Y));
943: else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
944: PetscCall(PCMatApply(ctx->pc, X, Y));
945: PetscCall(PCHPDDMMatApply_Private<false>(ctx, Y, &reset));
946: PetscCall(MatProductNumeric(ctx->V[1]));
947: PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN));
948: PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN));
949: PetscCall(PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2]));
950: PetscCall(MatAXPY(Y, -1.0, ctx->V[2], SAME_NONZERO_PATTERN));
951: } else {
952: PetscCall(PCHPDDMMatApply_Private<false>(ctx, Y, &reset));
953: PetscCall(PCHPDDMDeflate_Private(pc, X, Y));
954: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
955: PetscCall(MatProductNumeric(ctx->V[1]));
956: PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN));
957: PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN));
958: PetscCall(PCMatApply(ctx->pc, ctx->V[2], ctx->V[1]));
959: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
960: PetscCall(MatProductNumeric(ctx->V[2]));
961: PetscCall(PCHPDDMDeflate_Private<true>(pc, ctx->V[2], ctx->V[2]));
962: PetscCall(MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN));
963: }
964: PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN));
965: } else {
966: PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
967: PetscCall(PCMatApply(ctx->pc, X, ctx->V[1]));
968: PetscCall(MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN));
969: }
970: }
971: if (reset) PetscCall(MatDenseResetArray(ctx->V[1]));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: static PetscErrorCode PCApplyTranspose_HPDDMShell(PC pc, Vec x, Vec y)
976: {
977: PC_HPDDM_Level *ctx;
978: Mat A;
980: PetscFunctionBegin;
981: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
982: PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
983: PetscCheck(!ctx->parent->normal, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not implemented for the normal equations");
984: PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr));
985: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCApplyTranspose(ctx->pc, x, y)); /* y = M^-T x */
986: else {
987: PetscCall(PCHPDDMDeflate_Private<true>(pc, x, y)); /* y = Q^T x */
988: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
989: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
990: /* TODO: checking whether Q^T = Q would make it possible to skip this coarse correction */
991: PetscCall(PCHPDDMDeflate_Private(pc, x, ctx->v[1][1])); /* y = Q x */
992: PetscCall(MatMult(A, ctx->v[1][1], ctx->v[1][0])); /* y = A Q x */
993: PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* y = (I - A Q) x */
994: PetscCall(PCApplyTranspose(ctx->pc, ctx->v[1][1], ctx->v[1][0])); /* y = M^-T (I - A Q) x */
995: } else PetscCall(PCApplyTranspose(ctx->pc, x, ctx->v[1][0])); /* y = M^-T x */
996: PetscCall(MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1])); /* z = A^T y */
997: PetscCall(PCHPDDMDeflate_Private<true>(pc, ctx->v[1][1], ctx->v[1][1])); /* z = Q^T z */
998: PetscCall(VecAXPBYPCZ(y, -1.0, 1.0, 1.0, ctx->v[1][1], ctx->v[1][0])); /* y = (I - Q^T A^T) y + Q^T x */
999: } else {
1000: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
1001: PetscCall(MatMultHermitianTranspose(A, y, ctx->v[1][0]));
1002: PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x));
1003: } else PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
1004: PetscCall(PCApplyTranspose(ctx->pc, ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE ? x : ctx->v[1][1], ctx->v[1][0]));
1005: PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-T x + Q^T x or M^-T (I - A^T Q^T) x + Q^T x */
1006: }
1007: }
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*
1012: PCMatApplyTranspose_HPDDMShell - Variant of PCApplyTranspose_HPDDMShell() for blocks of vectors.
1014: Input Parameters:
1015: + pc - preconditioner context
1016: - X - block of input vectors
1018: Output Parameter:
1019: . Y - block of output vectors
1021: Level: advanced
1023: .seealso: [](ch_ksp), `PCHPDDM`, `PCApplyTranspose_HPDDMShell()`, `PCHPDDMCoarseCorrectionType`
1024: */
1025: static PetscErrorCode PCMatApplyTranspose_HPDDMShell(PC pc, Mat X, Mat Y)
1026: {
1027: PC_HPDDM_Level *ctx;
1028: PetscBool reset = PETSC_FALSE;
1030: PetscFunctionBegin;
1031: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
1032: PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object");
1033: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCMatApplyTranspose(ctx->pc, X, Y));
1034: else if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
1035: /* similar code as in PCMatApply_HPDDMShell() with an extra call to PCHPDDMDeflate_Private<true>() */
1036: PetscCall(PCHPDDMMatApply_Private<false>(ctx, Y, &reset));
1037: PetscCall(PCHPDDMDeflate_Private(pc, X, Y));
1038: PetscCall(MatProductNumeric(ctx->V[1]));
1039: PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN));
1040: PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN));
1041: PetscCall(PCMatApplyTranspose(ctx->pc, ctx->V[2], ctx->V[1]));
1042: PetscCall(MatProductNumeric(ctx->V[2]));
1043: PetscCall(PCHPDDMDeflate_Private<true>(pc, ctx->V[2], ctx->V[2]));
1044: PetscCall(MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN));
1045: PetscCall(PCHPDDMDeflate_Private<true>(pc, X, Y)); /* TODO: checking whether Q^T = Q would make it possible to skip this coarse correction */
1046: PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN));
1047: } else {
1048: PetscCall(PCHPDDMMatApply_Private<true>(ctx, Y, &reset));
1049: PetscCall(PCHPDDMDeflate_Private<true>(pc, X, Y));
1050: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) {
1051: PetscCall(PCMatApplyTranspose(ctx->pc, X, ctx->V[2]));
1052: PetscCall(MatAXPY(Y, 1.0, ctx->V[2], SAME_NONZERO_PATTERN));
1053: PetscCall(MatProductNumeric(ctx->V[1]));
1054: /* ctx->V[0] and ctx->V[1] memory regions overlap, so need to copy to ctx->V[2] and switch array */
1055: PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN));
1056: if (reset) PetscCall(MatDenseResetArray(ctx->V[1]));
1057: PetscCall(PCHPDDMDeflate_Private<true>(pc, ctx->V[2], ctx->V[2]));
1058: PetscCall(MatAXPY(Y, -1.0, ctx->V[2], SAME_NONZERO_PATTERN));
1059: } else {
1060: if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
1061: PetscCall(MatCopy(Y, ctx->V[2], SAME_NONZERO_PATTERN));
1062: PetscCall(MatProductNumeric(ctx->V[1]));
1063: PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN));
1064: PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN));
1065: PetscCall(PCMatApplyTranspose(ctx->pc, ctx->V[2], ctx->V[1]));
1066: PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN));
1067: } else {
1068: PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction);
1069: PetscCall(PCMatApplyTranspose(ctx->pc, X, ctx->V[1]));
1070: PetscCall(MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN));
1071: }
1072: if (reset) PetscCall(MatDenseResetArray(ctx->V[1]));
1073: }
1074: }
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: static PetscErrorCode PCDestroy_HPDDMShell(PC pc)
1079: {
1080: PC_HPDDM_Level *ctx;
1082: PetscFunctionBegin;
1083: PetscCall(PCShellGetContext(pc, static_cast<void *>(&ctx)));
1084: PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE));
1085: PetscCall(VecDestroyVecs(1, &ctx->v[0]));
1086: PetscCall(VecDestroyVecs(2, &ctx->v[1]));
1087: PetscCall(PetscObjectCompose((PetscObject)ctx->pc->mat, "_HPDDM_MatProduct", nullptr));
1088: PetscCall(MatDestroy(ctx->V));
1089: PetscCall(MatDestroy(ctx->V + 1));
1090: PetscCall(MatDestroy(ctx->V + 2));
1091: PetscCall(VecDestroy(&ctx->D));
1092: PetscCall(PetscSFDestroy(&ctx->scatter));
1093: PetscCall(PCDestroy(&ctx->pc));
1094: PetscFunctionReturn(PETSC_SUCCESS);
1095: }
1097: template <class Type, bool T = false, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
1098: static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type x, Type y)
1099: {
1100: PetscFunctionBegin;
1101: PetscCall(VecISCopy(std::get<2>(*p)[0], std::get<1>(*p), SCATTER_FORWARD, x));
1102: if (!T) PetscCall(PCApply(factor, std::get<2>(*p)[0], std::get<2>(*p)[1]));
1103: else PetscCall(PCApplyTranspose(factor, std::get<2>(*p)[0], std::get<2>(*p)[1]));
1104: PetscCall(VecISCopy(std::get<2>(*p)[1], std::get<1>(*p), SCATTER_REVERSE, y));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: template <class Type, bool = false, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
1109: static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type X, Type Y)
1110: {
1111: Mat B[2];
1112: Vec x, y;
1114: PetscFunctionBegin;
1115: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B));
1116: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B + 1));
1117: for (PetscInt i = 0; i < X->cmap->n; ++i) {
1118: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
1119: PetscCall(MatDenseGetColumnVecWrite(B[0], i, &y));
1120: PetscCall(VecISCopy(y, std::get<1>(*p), SCATTER_FORWARD, x));
1121: PetscCall(MatDenseRestoreColumnVecWrite(B[0], i, &y));
1122: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
1123: }
1124: PetscCall(PCMatApply(factor, B[0], B[1]));
1125: PetscCall(MatDestroy(B));
1126: for (PetscInt i = 0; i < X->cmap->n; ++i) {
1127: PetscCall(MatDenseGetColumnVecRead(B[1], i, &x));
1128: PetscCall(MatDenseGetColumnVecWrite(Y, i, &y));
1129: PetscCall(VecISCopy(x, std::get<1>(*p), SCATTER_REVERSE, y));
1130: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y));
1131: PetscCall(MatDenseRestoreColumnVecRead(B[1], i, &x));
1132: }
1133: PetscCall(MatDestroy(B + 1));
1134: PetscFunctionReturn(PETSC_SUCCESS);
1135: }
1137: template <class Type = Vec, bool T = false>
1138: static PetscErrorCode PCApply_Schur(PC pc, Type x, Type y)
1139: {
1140: PC factor;
1141: Mat A;
1142: MatSolverType type;
1143: PetscBool flg;
1144: std::tuple<KSP, IS, Vec[2]> *p;
1146: PetscFunctionBegin;
1147: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1148: PetscCall(KSPGetPC(std::get<0>(*p), &factor));
1149: PetscCall(PCFactorGetMatSolverType(factor, &type));
1150: PetscCall(PCFactorGetMatrix(factor, &A));
1151: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
1152: if (flg) {
1153: PetscCheck(PetscDefined(HAVE_MUMPS), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType");
1154: PetscCall(MatMumpsSetIcntl(A, 26, 0));
1155: } else {
1156: PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &flg));
1157: PetscCheck(flg && PetscDefined(HAVE_MKL_PARDISO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType");
1158: flg = PETSC_FALSE;
1159: #if PetscDefined(HAVE_MKL_PARDISO)
1160: PetscCall(MatMkl_PardisoSetCntl(A, 70, 1));
1161: #endif
1162: }
1163: PetscCall(PCApply_Schur_Private<Type, T>(p, factor, x, y));
1164: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1));
1165: else {
1166: #if PetscDefined(HAVE_MKL_PARDISO)
1167: PetscCall(MatMkl_PardisoSetCntl(A, 70, 0));
1168: #endif
1169: }
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: static PetscErrorCode PCDestroy_Schur(PC pc)
1174: {
1175: std::tuple<KSP, IS, Vec[2]> *p;
1177: PetscFunctionBegin;
1178: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1179: PetscCall(ISDestroy(&std::get<1>(*p)));
1180: PetscCall(VecDestroy(std::get<2>(*p)));
1181: PetscCall(VecDestroy(std::get<2>(*p) + 1));
1182: PetscCall(PetscFree(p));
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: template <bool transpose>
1187: static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu)
1188: {
1189: Mat B, X;
1190: PetscInt n, N, j = 0;
1192: PetscFunctionBegin;
1193: PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr));
1194: PetscCall(MatGetLocalSize(B, &n, nullptr));
1195: PetscCall(MatGetSize(B, &N, nullptr));
1196: if (ctx->parent->log_separate) {
1197: j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx));
1198: PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr));
1199: }
1200: if (mu == 1) {
1201: if (!ctx->ksp->vec_rhs) {
1202: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs));
1203: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol));
1204: }
1205: PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs));
1206: if (!transpose) PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr));
1207: else {
1208: PetscCall(VecConjugate(ctx->ksp->vec_rhs));
1209: PetscCall(KSPSolveTranspose(ctx->ksp, nullptr, nullptr)); /* TODO: missing KSPSolveHermitianTranspose() */
1210: PetscCall(VecConjugate(ctx->ksp->vec_sol));
1211: }
1212: PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs));
1213: PetscCall(VecResetArray(ctx->ksp->vec_rhs));
1214: } else {
1215: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B));
1216: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X));
1217: if (!transpose) PetscCall(KSPMatSolve(ctx->ksp, B, X));
1218: else {
1219: PetscCall(MatConjugate(B));
1220: PetscCall(KSPMatSolveTranspose(ctx->ksp, B, X)); /* TODO: missing KSPMatSolveHermitianTranspose() */
1221: PetscCall(MatConjugate(X));
1222: }
1223: PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN));
1224: PetscCall(MatDestroy(&X));
1225: PetscCall(MatDestroy(&B));
1226: }
1227: if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr));
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc)
1232: {
1233: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1235: PetscFunctionBegin;
1236: if (data->setup) {
1237: Mat P;
1238: Vec x, xt = nullptr;
1239: PetscReal t = 0.0, s = 0.0;
1241: PetscCall(PCGetOperators(pc, nullptr, &P));
1242: PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x));
1243: PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx));
1244: }
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[])
1249: {
1250: Mat A;
1251: PetscBool flg;
1253: PetscFunctionBegin;
1254: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n);
1255: /* previously composed Mat */
1256: PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A));
1257: PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat");
1258: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSCHURCOMPLEMENT, &flg)); /* MATSCHURCOMPLEMENT has neither a MatDuplicate() nor a MatCopy() implementation */
1259: if (scall == MAT_INITIAL_MATRIX) {
1260: PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */
1261: if (!flg) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat));
1262: } else if (!flg) PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN));
1263: if (flg) {
1264: PetscCall(MatDestroy(*submat)); /* previously created Mat has to be destroyed */
1265: (*submat)[0] = A;
1266: PetscCall(PetscObjectReference((PetscObject)A));
1267: }
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted)
1272: {
1273: PetscErrorCodeFn *op;
1275: PetscFunctionBegin;
1276: /* previously-composed Mat */
1277: PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C));
1278: PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op));
1279: /* see https://mailman.cels.anl.gov/archives/list/petsc-dev@lists.mcs.anl.gov/message/22HXNMER6OU7N7CXWA2LJAY6RZPGQYXT/ */
1280: PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (PetscErrorCodeFn *)PCHPDDMCreateSubMatrices_Private));
1281: if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */
1282: PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */
1283: PetscCall(PCSetUp(pc));
1284: /* reset MatCreateSubMatrices() */
1285: PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op));
1286: PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p)
1291: {
1292: IS perm;
1293: const PetscInt *ptr;
1294: PetscInt *compressed, size, bs;
1295: std::map<PetscInt, PetscInt> order;
1296: PetscBool flg;
1298: PetscFunctionBegin;
1301: PetscCall(ISGetLocalSize(is, &size));
1302: PetscCall(ISGetBlockSize(is, &bs));
1303: PetscCall(ISSorted(is, &flg));
1304: if (!flg) {
1305: PetscCall(ISGetIndices(is, &ptr));
1306: /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */
1307: for (PetscInt n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs));
1308: PetscCall(ISRestoreIndices(is, &ptr));
1309: size /= bs;
1310: if (out_C) {
1311: PetscCall(PetscMalloc1(size, &compressed));
1312: for (const std::pair<const PetscInt, PetscInt> &i : order) *compressed++ = i.second;
1313: compressed -= size;
1314: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, compressed, PETSC_OWN_POINTER, &perm));
1315: PetscCall(ISSetPermutation(perm));
1316: /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */
1317: PetscCall(MatPermute(in_C, perm, perm, out_C));
1318: if (p) *p = perm;
1319: else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */
1320: }
1321: if (out_is) {
1322: PetscCall(PetscMalloc1(size, &compressed));
1323: for (const std::pair<const PetscInt, PetscInt> &i : order) *compressed++ = i.first;
1324: compressed -= size;
1325: /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */
1326: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, compressed, PETSC_OWN_POINTER, out_is));
1327: }
1328: } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */
1329: if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C));
1330: if (out_is) {
1331: PetscCall(PetscObjectTypeCompare((PetscObject)in_is, ISBLOCK, &flg));
1332: if (flg) PetscCall(ISDuplicate(in_is, out_is));
1333: else {
1334: PetscCall(ISGetIndices(is, &ptr));
1335: if (bs > 1) {
1336: size /= bs;
1337: PetscCall(PetscMalloc1(size, &compressed));
1338: for (PetscInt n = 0; n < size; ++n) compressed[n] = ptr[n * bs] / bs;
1339: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, compressed, PETSC_OWN_POINTER, out_is));
1340: } else PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), 1, size, ptr, PETSC_COPY_VALUES, out_is));
1341: PetscCall(ISRestoreIndices(is, &ptr));
1342: PetscCall(ISSetInfo(*out_is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1343: }
1344: }
1345: }
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10, Mat *B01 = nullptr)
1350: {
1351: Mat T, U = nullptr, B = nullptr;
1352: IS z;
1353: PetscBool flg, conjugate = PETSC_FALSE;
1355: PetscFunctionBegin;
1356: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
1357: if (B01) *B01 = nullptr;
1358: if (flg) {
1359: PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1360: PetscCall(MatTransposeGetMat(A10, &U));
1361: } else {
1362: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1363: if (flg) {
1364: PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1365: PetscCall(MatHermitianTransposeGetMat(A10, &U));
1366: conjugate = PETSC_TRUE;
1367: }
1368: }
1369: if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T));
1370: else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T));
1371: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
1372: if (flg) {
1373: PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1374: PetscCall(MatTransposeGetMat(A01, &A01));
1375: PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B));
1376: A01 = B;
1377: } else {
1378: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1379: if (flg) {
1380: PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1381: PetscCall(MatHermitianTransposeGetMat(A01, &A01));
1382: PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B));
1383: A01 = B;
1384: }
1385: }
1386: PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, &flg));
1387: if (flg) {
1388: PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, &flg));
1389: if (flg) {
1390: PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */
1391: if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */
1392: if (B01) PetscCall(MatDuplicate(T, MAT_COPY_VALUES, B01));
1393: PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE));
1394: PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */
1395: }
1396: PetscCall(MatMultEqual(A01, T, 20, &flg));
1397: if (!B01) PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T");
1398: else {
1399: PetscCall(PetscInfo(pc, "A01 and A10^T are equal? %s\n", PetscBools[flg]));
1400: if (!flg) {
1401: if (z) PetscCall(MatDestroy(&T));
1402: else *B01 = T;
1403: flg = PETSC_TRUE;
1404: } else PetscCall(MatDestroy(B01));
1405: }
1406: PetscCall(ISDestroy(&z));
1407: }
1408: }
1409: if (!flg) PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent layouts, cannot test for equality\n"));
1410: if (!B01 || !*B01) PetscCall(MatDestroy(&T));
1411: else if (conjugate) PetscCall(MatConjugate(T));
1412: PetscCall(MatDestroy(&B));
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: static PetscErrorCode PCHPDDMCheckInclusion_Private(PC pc, IS is, IS is_local, PetscBool check)
1417: {
1418: IS intersect;
1419: const char *str = "IS of the auxiliary Mat does not include all local rows of A";
1420: PetscBool equal;
1422: PetscFunctionBegin;
1423: PetscCall(ISIntersect(is, is_local, &intersect));
1424: PetscCall(ISEqualUnsorted(is_local, intersect, &equal));
1425: PetscCall(ISDestroy(&intersect));
1426: if (check) PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", str);
1427: else if (!equal) PetscCall(PetscInfo(pc, "%s\n", str));
1428: PetscFunctionReturn(PETSC_SUCCESS);
1429: }
1431: static PetscErrorCode PCHPDDMCheckMatStructure_Private(PC pc, Mat A, Mat B)
1432: {
1433: Mat X, Y;
1434: const PetscInt *i[2], *j[2];
1435: PetscBool flg = PETSC_TRUE;
1437: PetscFunctionBegin;
1438: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &X)); /* no common way to compare sparsity pattern, so just convert to MATSEQAIJ */
1439: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &Y)); /* the second Mat (B = Neumann) should have a SUBSET_NONZERO_PATTERN MatStructure of the first one (A = Dirichlet) */
1440: PetscCall(MatSeqAIJGetCSRAndMemType(X, &i[0], &j[0], nullptr, nullptr));
1441: PetscCall(MatSeqAIJGetCSRAndMemType(Y, &i[1], &j[1], nullptr, nullptr));
1442: for (PetscInt row = 0; (row < X->rmap->n) && flg; ++row) {
1443: const PetscInt n = i[0][row + 1] - i[0][row];
1445: for (PetscInt k = i[1][row], location; k < i[1][row + 1]; ++k) {
1446: PetscCall(PetscFindInt(j[1][k], n, j[0] + i[0][row], &location));
1447: if (location < 0) {
1448: flg = PETSC_FALSE;
1449: break;
1450: }
1451: }
1452: }
1453: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)pc)));
1454: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER_INPUT, "Auxiliary Mat is supposedly the local Neumann matrix but it has a sparsity pattern which is not a subset of the one of the local assembled matrix");
1455: PetscCall(MatDestroy(&Y));
1456: PetscCall(MatDestroy(&X));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub)
1461: {
1462: IS is;
1464: PetscFunctionBegin;
1465: if (!flg) {
1466: if (algebraic) {
1467: PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is));
1468: PetscCall(ISDestroy(&is));
1469: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr));
1470: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr));
1471: }
1472: PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub));
1473: }
1474: PetscFunctionReturn(PETSC_SUCCESS);
1475: }
1477: static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat Q, IS *is, Mat *sub[], PetscBool block)
1478: {
1479: IS icol[3], irow[2];
1480: Mat *M;
1481: Mat P = Q;
1482: PetscReal *ptr;
1483: PetscInt *idx, p = 0, bs = P->cmap->bs;
1484: PetscBool flg;
1486: PetscFunctionBegin;
1487: /* MatCreateSubMatrices_MPISBAIJ() may return a rectangular MATSEQSBAIJ containing only the explicitly stored upper-triangular entries */
1488: /* of the selected rows. But the missing lower-triangular entries are needed by MatGetColumnNorms(), MatGetRowSum(), and MatMatMult() */
1489: PetscCall(PetscObjectTypeCompare((PetscObject)Q, MATMPISBAIJ, &flg));
1490: if (flg) PetscCall(MatConvert(Q, MATBAIJ, MAT_INITIAL_MATRIX, &P));
1491: PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2));
1492: PetscCall(ISSetBlockSize(icol[2], bs));
1493: PetscCall(ISSetIdentity(icol[2]));
1494: PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M));
1495: if (flg) {
1496: PetscCall(MatDestroy(&P));
1497: P = Q; /* continue using the caller-owned (MATMPISBAIJ) matrix */
1498: }
1499: PetscCall(ISDestroy(icol + 2));
1500: PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow));
1501: PetscCall(ISSetBlockSize(irow[0], bs));
1502: PetscCall(ISSetIdentity(irow[0]));
1503: if (!block) {
1504: PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx));
1505: PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr));
1506: /* check for nonzero columns so that M[0] may be expressed in compact form */
1507: for (PetscInt n = 0; n < P->cmap->N; n += bs) {
1508: if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs;
1509: }
1510: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1));
1511: PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1512: PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2));
1513: irow[1] = irow[0];
1514: /* 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 */
1515: icol[0] = is[0];
1516: PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub));
1517: PetscCall(ISDestroy(icol + 1));
1518: PetscCall(PetscFree2(ptr, idx));
1519: PetscCall(MatPropagateSymmetryOptions(P, (*sub)[0]));
1520: if (flg) PetscCall(MatConvert((*sub)[0], MATSBAIJ, MAT_INPLACE_MATRIX, sub[0]));
1521: /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */
1522: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2]));
1523: /* Mat used in eq. (3.1) of [2022b] */
1524: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1]));
1525: } else {
1526: Mat aux;
1528: PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
1529: /* diagonal block of the overlapping rows */
1530: PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub));
1531: PetscCall(MatPropagateSymmetryOptions(P, (*sub)[0]));
1532: if (flg && bs == 1) PetscCall(MatConvert((*sub)[0], MATSBAIJ, MAT_INPLACE_MATRIX, sub[0]));
1533: PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux));
1534: aux->spd = PETSC_BOOL3_UNKNOWN; /* the auxiliary Mat need not be SPD */
1535: PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1536: if (bs == 1) { /* scalar case */
1537: Vec sum[2];
1539: PetscCall(MatCreateVecs(aux, sum, sum + 1));
1540: PetscCall(MatGetRowSum(M[0], sum[0]));
1541: PetscCall(MatGetRowSum(aux, sum[1]));
1542: /* off-diagonal block row sum (full rows - diagonal block rows) */
1543: PetscCall(VecAXPY(sum[0], -1.0, sum[1]));
1544: /* subdomain matrix plus off-diagonal block row sum */
1545: PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES));
1546: PetscCall(VecDestroy(sum));
1547: PetscCall(VecDestroy(sum + 1));
1548: } else { /* vectorial case */
1549: /* TODO: missing MatGetValuesBlocked(), so the code below is */
1550: /* an extension of the scalar case for when bs > 1, but it could */
1551: /* be more efficient by avoiding all these MatMatMult() */
1552: Mat sum[2], ones;
1553: PetscScalar *ptr;
1555: aux->symmetry_eternal = PETSC_FALSE;
1556: aux->symmetric = PETSC_BOOL3_UNKNOWN;
1557: aux->hermitian = PETSC_BOOL3_UNKNOWN;
1558: PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr));
1559: PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones));
1560: for (PetscInt n = 0; n < M[0]->cmap->n; n += bs) {
1561: for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0;
1562: }
1563: PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum));
1564: PetscCall(MatDestroy(&ones));
1565: PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones));
1566: PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n));
1567: PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum + 1));
1568: PetscCall(MatDestroy(&ones));
1569: PetscCall(PetscFree(ptr));
1570: /* off-diagonal block row sum (full rows - diagonal block rows) */
1571: PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN));
1572: PetscCall(MatDestroy(sum + 1));
1573: /* re-order values to be consistent with MatSetValuesBlocked() */
1574: /* equivalent to MatTranspose() which does not truly handle */
1575: /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */
1576: PetscCall(MatDenseGetArrayWrite(sum[0], &ptr));
1577: HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs);
1578: /* subdomain matrix plus off-diagonal block row sum */
1579: for (PetscInt n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES));
1580: PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY));
1581: PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY));
1582: PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr));
1583: PetscCall(MatDestroy(sum));
1584: }
1585: PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1586: /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */
1587: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux));
1588: }
1589: PetscCall(ISDestroy(irow));
1590: PetscCall(MatDestroySubMatrices(1, &M));
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y)
1595: {
1596: Mat A;
1597: MatSolverType type;
1598: IS is[2];
1599: PetscBool flg;
1600: std::pair<PC, Vec[2]> *p;
1602: PetscFunctionBegin;
1603: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1604: if (p->second[0]) { /* in case of a centralized Schur complement, some processes may have no local operator */
1605: PetscCall(PCGetOperators(p->first, &A, nullptr));
1606: PetscCall(MatNestGetISs(A, is, nullptr));
1607: PetscCall(PetscObjectTypeCompareAny((PetscObject)p->first, &flg, PCLU, PCCHOLESKY, ""));
1608: if (flg) { /* partial solve currently only makes sense with exact factorizations */
1609: PetscCall(PCFactorGetMatSolverType(p->first, &type));
1610: PetscCall(PCFactorGetMatrix(p->first, &A));
1611: if (A->schur) {
1612: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
1613: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */
1614: } else flg = PETSC_FALSE;
1615: }
1616: PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */
1617: PetscCall(PCApply(p->first, p->second[0], p->second[1]));
1618: PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */
1619: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */
1620: }
1621: PetscFunctionReturn(PETSC_SUCCESS);
1622: }
1624: static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer)
1625: {
1626: std::pair<PC, Vec[2]> *p;
1628: PetscFunctionBegin;
1629: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1630: PetscCall(PCView(p->first, viewer));
1631: PetscFunctionReturn(PETSC_SUCCESS);
1632: }
1634: static PetscErrorCode PCDestroy_Nest(PC pc)
1635: {
1636: std::pair<PC, Vec[2]> *p;
1638: PetscFunctionBegin;
1639: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1640: PetscCall(VecDestroy(p->second));
1641: PetscCall(VecDestroy(p->second + 1));
1642: PetscCall(PCDestroy(&p->first));
1643: PetscCall(PetscFree(p));
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: template <bool T = false>
1648: static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y)
1649: {
1650: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
1652: PetscFunctionBegin;
1653: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1654: PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */
1655: PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD));
1656: if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */
1657: else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1]));
1658: PetscCall(VecSet(y, 0.0));
1659: PetscCall(VecScatterBegin(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); /* global Vec with summed up contributions on the overlap */
1660: PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE));
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: static PetscErrorCode MatDestroy_Schur(Mat A)
1665: {
1666: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
1668: PetscFunctionBegin;
1669: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1670: PetscCall(VecDestroy(std::get<2>(*ctx)));
1671: PetscCall(VecDestroy(std::get<2>(*ctx) + 1));
1672: PetscCall(PetscFree(ctx));
1673: PetscFunctionReturn(PETSC_SUCCESS);
1674: }
1676: static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y)
1677: {
1678: PC pc;
1679: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1681: PetscFunctionBegin;
1682: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1683: pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc;
1684: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { /* Q_0 is the coarse correction associated to the A00 block from PCFIELDSPLIT */
1685: PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */
1686: PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */
1687: PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */
1688: PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */
1689: } else {
1690: PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */
1691: PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */
1692: PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */
1693: PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */
1694: }
1695: PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */
1696: PetscFunctionReturn(PETSC_SUCCESS);
1697: }
1699: static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer)
1700: {
1701: PetscBool ascii;
1702: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1704: PetscFunctionBegin;
1705: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
1706: if (ascii) {
1707: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1708: PetscCall(PetscViewerASCIIPrintf(viewer, "action of %s\n", std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT ? "(I - M_S^-1 A_10 Q_0 A_01)" : "(I - A_10 Q_0 A_01 M_S^-1)"));
1709: PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */
1710: }
1711: PetscFunctionReturn(PETSC_SUCCESS);
1712: }
1714: static PetscErrorCode MatDestroy_SchurCorrection(Mat A)
1715: {
1716: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1718: PetscFunctionBegin;
1719: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1720: PetscCall(VecDestroy(std::get<3>(*ctx)));
1721: PetscCall(VecDestroy(std::get<3>(*ctx) + 1));
1722: PetscCall(VecDestroy(std::get<3>(*ctx) + 2));
1723: PetscCall(PCDestroy(std::get<0>(*ctx) + 1));
1724: PetscCall(PetscFree(ctx));
1725: PetscFunctionReturn(PETSC_SUCCESS);
1726: }
1728: static PetscErrorCode PCPostSolve_SchurPreLeastSquares(PC, KSP, Vec, Vec x)
1729: {
1730: PetscFunctionBegin;
1731: PetscCall(VecScale(x, -1.0));
1732: PetscFunctionReturn(PETSC_SUCCESS);
1733: }
1735: static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context)
1736: {
1737: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context);
1739: PetscFunctionBegin;
1740: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) {
1741: PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2]));
1742: std::swap(*b, *std::get<3>(*ctx)[2]); /* replace b by M^-1 b, but need to keep a copy of the original RHS, so swap it with the work Vec */
1743: }
1744: PetscFunctionReturn(PETSC_SUCCESS);
1745: }
1747: static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context)
1748: {
1749: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context);
1751: PetscFunctionBegin;
1752: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) std::swap(*b, *std::get<3>(*ctx)[2]); /* put back the original RHS where it belongs */
1753: else {
1754: PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2]));
1755: PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */
1756: }
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec);
1761: static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec);
1762: static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *);
1763: static PetscErrorCode MatProduct_AtB_Harmonic(Mat, Mat, Mat, void *);
1764: static PetscErrorCode MatDestroy_Harmonic(Mat);
1766: static PetscErrorCode PCSetUp_HPDDM(PC pc)
1767: {
1768: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1769: PC inner;
1770: KSP *ksp = nullptr;
1771: Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S;
1772: Vec xin, v;
1773: std::vector<Vec> initial;
1774: IS is[1], loc, uis = data->is, unsorted = nullptr;
1775: ISLocalToGlobalMapping l2g;
1776: char prefix[256];
1777: const char *pcpre;
1778: Mat ev;
1779: PetscInt n, requested, reused = 0, overlap = -1;
1780: MatStructure structure = UNKNOWN_NONZERO_PATTERN;
1781: PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE;
1782: DM dm;
1783: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr;
1784: IS dis = nullptr;
1785: Mat daux = nullptr;
1787: PetscFunctionBegin;
1788: if (!data->levels) PetscCall(PetscInfo(pc, "No level allocated, defaulting to a single level, PCSetFromOptions() should be called before PCSetUp() to avoid this\n"));
1789: PetscCall(PCHPDDMInitializeLevels_Private(data));
1790: requested = data->N;
1791: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
1792: PetscCall(PCGetOperators(pc, &A, &P));
1793: if (!data->levels[0]->ksp) {
1794: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp));
1795: PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel));
1796: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse"));
1797: PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix));
1798: PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY));
1799: } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled && data->levels[0]->ksp->pc->reusepreconditioner) {
1800: /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */
1801: /* then just propagate the appropriate flag to the coarser levels */
1802: for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1803: /* the following KSP and PC may be NULL for some processes, hence the check */
1804: if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE));
1805: if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE));
1806: }
1807: /* early bail out because there is nothing to do */
1808: PetscFunctionReturn(PETSC_SUCCESS);
1809: } else {
1810: /* reset coarser levels */
1811: for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1812: if (data->levels[n]->ksp && data->levels[n]->ksp->pc && data->levels[n]->ksp->pc->setupcalled && data->levels[n]->ksp->pc->reusepreconditioner && n < data->N) {
1813: reused = data->N - n;
1814: break;
1815: }
1816: PetscCall(KSPDestroy(&data->levels[n]->ksp));
1817: PetscCall(PCDestroy(&data->levels[n]->pc));
1818: }
1819: /* check if some coarser levels are being reused */
1820: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1821: const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &i_0;
1823: if (*addr != 0 && reused != data->N - 1) {
1824: /* reuse previously computed eigenvectors */
1825: ev = data->levels[0]->P->getMat();
1826: if (ev) {
1827: initial.reserve(*addr);
1828: for (n = 0; n < *addr; ++n) {
1829: PetscCall(MatDenseGetColumnVecRead(ev, n, &xin));
1830: PetscCall(VecDuplicate(xin, &v));
1831: PetscCall(VecCopy(xin, v));
1832: initial.emplace_back(v);
1833: PetscCall(MatDenseRestoreColumnVecRead(ev, n, &xin));
1834: }
1835: }
1836: }
1837: }
1838: data->N -= reused;
1839: PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P));
1841: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis));
1842: if (!data->is && !ismatis) {
1843: PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr;
1844: PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr;
1845: void *uctx = nullptr;
1847: /* first see if we can get the data from the DM */
1848: PetscCall(MatGetDM(P, &dm));
1849: if (!dm) PetscCall(MatGetDM(A, &dm));
1850: if (!dm) PetscCall(PCGetDM(pc, &dm));
1851: if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */
1852: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create));
1853: if (create) {
1854: PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx));
1855: if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */
1856: }
1857: }
1858: if (!create) {
1859: if (!uis) {
1860: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis));
1861: PetscCall(PetscObjectReference((PetscObject)uis));
1862: }
1863: if (!uaux) {
1864: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux));
1865: PetscCall(PetscObjectReference((PetscObject)uaux));
1866: }
1867: /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */
1868: if (!uis) {
1869: PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis));
1870: PetscCall(PetscObjectReference((PetscObject)uis));
1871: }
1872: if (!uaux) {
1873: PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux));
1874: PetscCall(PetscObjectReference((PetscObject)uaux));
1875: }
1876: }
1877: PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx));
1878: PetscCall(MatDestroy(&uaux));
1879: PetscCall(ISDestroy(&uis));
1880: }
1882: if (!ismatis) {
1883: PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc));
1884: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_block_splitting", &block, nullptr));
1885: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr));
1886: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg));
1887: if (data->is || flg) {
1888: if (block || overlap != -1) {
1889: PetscCall(ISDestroy(&data->is));
1890: PetscCall(MatDestroy(&data->aux));
1891: } else if (flg) {
1892: PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO;
1894: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg));
1895: if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) {
1896: PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */
1897: PetscCall(MatDestroy(&data->aux));
1898: } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) {
1899: PetscContainer container = nullptr;
1901: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container));
1902: if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */
1903: PC_HPDDM *data_00;
1904: KSP ksp, inner_ksp;
1905: PC pc_00;
1906: Mat A11;
1907: Vec d = nullptr;
1908: PetscReal norm;
1909: const PetscInt *ranges;
1910: PetscMPIInt size;
1911: char *prefix;
1913: PetscCall(MatSchurComplementGetKSP(P, &ksp));
1914: PetscCall(KSPGetPC(ksp, &pc_00));
1915: PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg));
1916: PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)pc_00)->prefix ? ((PetscObject)pc_00)->prefix : "",
1917: ((PetscObject)pc_00)->type_name, PCHPDDM);
1918: data_00 = (PC_HPDDM *)pc_00->data;
1919: PetscCheck(data_00->N == 2, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and %" PetscInt_FMT " level%s instead of 2 for the A00 block -%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type],
1920: data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix);
1921: PetscCheck(data_00->levels[0]->pc, PetscObjectComm((PetscObject)P), PETSC_ERR_ORDER, "PC of the first block%s not setup yet", ((PetscObject)pc_00)->prefix ? std::string(std::string(" (") + ((PetscObject)pc_00)->prefix + std::string(")")).c_str() : "");
1922: PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg));
1923: PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)data_00->levels[0]->pc)->prefix,
1924: ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM);
1925: PetscCall(PetscNew(&ctx)); /* context to pass data around for the inner-most PC, which will be a proper PCHPDDM (or a dummy variable if the Schur complement is centralized on a single process) */
1926: PetscCall(MatSchurComplementGetSubMatrices(P, nullptr, nullptr, nullptr, nullptr, &A11));
1927: PetscCall(MatGetOwnershipRanges(P, &ranges));
1928: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)P), &size));
1929: flg = PetscBool(std::find_if(ranges, ranges + size + 1, [&](PetscInt v) { return v != ranges[0] && v != ranges[size]; }) == ranges + size + 1); /* are all local matrices but one of dimension 0 (centralized Schur complement)? */
1930: if (!flg) {
1931: if (PetscDefined(USE_DEBUG) || !data->is) {
1932: Mat A01, A10, B = nullptr, C = nullptr, *sub;
1934: PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr));
1935: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
1936: if (flg) {
1937: PetscCall(MatTransposeGetMat(A10, &C));
1938: PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B));
1939: } else {
1940: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1941: if (flg) {
1942: PetscCall(MatHermitianTransposeGetMat(A10, &C));
1943: PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B));
1944: }
1945: }
1946: if (flg)
1947: PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1948: if (!B) {
1949: B = A10;
1950: PetscCall(PetscObjectReference((PetscObject)B));
1951: } else if (!data->is) {
1952: PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
1953: if (!flg) C = A01;
1954: else
1955: PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
1956: }
1957: PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis));
1958: PetscCall(ISSetIdentity(uis));
1959: if (!data->is) {
1960: if (!C) PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C));
1961: else PetscCall(PetscObjectReference((PetscObject)C));
1962: PetscCall(ISDuplicate(data_00->is, is));
1963: PetscCall(MatIncreaseOverlap(A, 1, is, 1));
1964: PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
1965: PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub));
1966: PetscCall(MatDestroy(&C));
1967: PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C));
1968: PetscCall(MatDestroySubMatrices(1, &sub));
1969: PetscCall(MatFindNonzeroRows(C, &data->is));
1970: PetscCheck(data->is, PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "No empty row, which likely means that some rows of A_10 are dense");
1971: PetscCall(MatDestroy(&C));
1972: PetscCall(ISDestroy(is));
1973: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc));
1974: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_FALSE));
1975: PetscCall(ISExpand(data->is, loc, is));
1976: PetscCall(ISDestroy(&loc));
1977: PetscCall(ISDestroy(&data->is));
1978: data->is = is[0];
1979: is[0] = nullptr;
1980: }
1981: if (PetscDefined(USE_DEBUG)) {
1982: PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10));
1983: PetscCall(MatCreateSubMatrices(B, 1, &uis, &data_00->is, MAT_INITIAL_MATRIX, &sub)); /* expensive check since all processes fetch all rows (but only some columns) of the constraint matrix */
1984: PetscCall(ISDestroy(&uis));
1985: PetscCall(ISDuplicate(data->is, &uis));
1986: PetscCall(ISSort(uis));
1987: PetscCall(ISComplement(uis, 0, B->rmap->N, is));
1988: PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C));
1989: PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr));
1990: PetscCall(ISDestroy(is));
1991: PetscCall(MatMultEqual(sub[0], C, 20, &flg));
1992: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The image of A_10 (R_i^p)^T from the local primal (e.g., velocity) space to the full dual (e.g., pressure) space is not restricted to the local dual space: A_10 (R_i^p)^T != R_i^d (R_i^d)^T A_10 (R_i^p)^T"); /* cf. eq. (9) of https://hal.science/hal-02343808v6/document */
1993: PetscCall(MatDestroy(&C));
1994: PetscCall(MatDestroySubMatrices(1, &sub));
1995: }
1996: PetscCall(ISDestroy(&uis));
1997: PetscCall(MatDestroy(&B));
1998: }
1999: flg = PETSC_FALSE;
2000: if (!data->aux && A11) {
2001: Mat D;
2003: PetscCall(MatCreateVecs(A11, &d, nullptr));
2004: PetscCall(MatGetDiagonal(A11, d));
2005: PetscCall(PetscObjectTypeCompareAny((PetscObject)A11, &flg, MATDIAGONAL, MATCONSTANTDIAGONAL, ""));
2006: if (!flg) {
2007: PetscCall(MatCreateDiagonal(d, &D));
2008: PetscCall(MatMultEqual(A11, D, 20, &flg));
2009: PetscCall(MatDestroy(&D));
2010: }
2011: if (flg) PetscCall(PetscInfo(pc, "A11 block is likely diagonal so the PC will build an auxiliary Mat (which was not initially provided by the user)\n"));
2012: }
2013: if (PetscDefined(USE_DEBUG) || (data->Neumann != PETSC_BOOL3_TRUE && !flg)) {
2014: if (A11) PetscCall(MatNorm(A11, NORM_INFINITY, &norm));
2015: else norm = 0.0;
2016: if (data->Neumann != PETSC_BOOL3_TRUE && !flg) {
2017: PetscCheck(norm < PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0), PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_has_neumann != true with a nonzero or non-diagonal A11 block", pcpre ? pcpre : "", pcpre ? pcpre : "");
2018: PetscCall(PetscInfo(pc, "A11 block is likely zero so the PC will build an auxiliary Mat (which was%s initially provided by the user)\n", data->aux ? "" : " not"));
2019: PetscCall(MatDestroy(&data->aux));
2020: flg = PETSC_TRUE;
2021: }
2022: }
2023: if (!data->aux) { /* if A11 is near zero, e.g., Stokes equation, or diagonal, build an auxiliary (Neumann) Mat which is a (possibly slightly shifted) diagonal weighted by the inverse of the multiplicity */
2024: PetscSF scatter;
2025: const PetscScalar *read;
2026: PetscScalar *write, *diagonal = nullptr;
2028: PetscCall(MatDestroy(&data->aux));
2029: PetscCall(ISGetLocalSize(data->is, &n));
2030: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &xin));
2031: PetscCall(VecDuplicate(xin, &v));
2032: PetscCall(VecScatterCreate(xin, data->is, v, nullptr, &scatter));
2033: PetscCall(VecSet(v, 1.0));
2034: PetscCall(VecSet(xin, 1.0));
2035: PetscCall(VecScatterBegin(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE));
2036: PetscCall(VecScatterEnd(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */
2037: PetscCall(PetscSFDestroy(&scatter));
2038: if (d) {
2039: PetscCall(VecScatterCreate(d, data->is, v, nullptr, &scatter));
2040: PetscCall(VecScatterBegin(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD));
2041: PetscCall(VecScatterEnd(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD));
2042: PetscCall(PetscSFDestroy(&scatter));
2043: PetscCall(VecDestroy(&d));
2044: PetscCall(PetscMalloc1(n, &diagonal));
2045: PetscCall(VecGetArrayRead(v, &read));
2046: PetscCallCXX(std::copy_n(read, n, diagonal));
2047: PetscCall(VecRestoreArrayRead(v, &read));
2048: }
2049: PetscCall(VecDestroy(&v));
2050: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v));
2051: PetscCall(VecGetArrayRead(xin, &read));
2052: PetscCall(VecGetArrayWrite(v, &write));
2053: for (PetscInt i = 0; i < n; ++i) write[i] = (!diagonal || std::abs(diagonal[i]) < PETSC_MACHINE_EPSILON) ? PETSC_SMALL / (static_cast<PetscReal>(1000.0) * read[i]) : diagonal[i] / read[i];
2054: PetscCall(PetscFree(diagonal));
2055: PetscCall(VecRestoreArrayRead(xin, &read));
2056: PetscCall(VecRestoreArrayWrite(v, &write));
2057: PetscCall(VecDestroy(&xin));
2058: PetscCall(MatCreateDiagonal(v, &data->aux));
2059: PetscCall(VecDestroy(&v));
2060: }
2061: uis = data->is;
2062: uaux = data->aux;
2063: PetscCall(PetscObjectReference((PetscObject)uis));
2064: PetscCall(PetscObjectReference((PetscObject)uaux));
2065: PetscCall(PetscStrallocpy(pcpre, &prefix));
2066: PetscCall(PCSetOptionsPrefix(pc, nullptr));
2067: PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */
2068: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */
2069: PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n));
2070: PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2));
2071: PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat));
2072: PetscCall(KSPSetOptionsPrefix(inner_ksp, prefix));
2073: PetscCall(KSPAppendOptionsPrefix(inner_ksp, "pc_hpddm_"));
2074: PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE));
2075: PetscCall(KSPSetFromOptions(inner_ksp));
2076: PetscCall(KSPGetPC(inner_ksp, &inner));
2077: PetscCall(PCSetOptionsPrefix(inner, nullptr));
2078: PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */
2079: PetscCall(PCKSPSetKSP(pc, inner_ksp));
2080: std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */
2081: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1]));
2082: PetscCall(PCSetOptionsPrefix(pc, prefix)); /* both PC share the same prefix so that the outer PC can be reset with PCSetFromOptions() */
2083: PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix));
2084: PetscCall(PetscFree(prefix));
2085: PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat));
2086: PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM));
2087: PetscCall(PCHPDDMSetAuxiliaryMat(std::get<0>(*ctx)[1], uis, uaux, nullptr, nullptr)); /* transfer ownership of the auxiliary inputs from the inner (PCKSP) to the inner-most (PCHPDDM) PC */
2088: if (flg) static_cast<PC_HPDDM *>(std::get<0>(*ctx)[1]->data)->Neumann = PETSC_BOOL3_TRUE;
2089: else if (PetscDefined(USE_DEBUG) && norm > PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0)) {
2090: /* no check when A11 is near zero */
2091: PetscCall(MatCreateSubMatrices(A11, 1, &uis, &uis, MAT_INITIAL_MATRIX, &sub));
2092: PetscCall(PCHPDDMCheckMatStructure_Private(pc, sub[0], uaux));
2093: PetscCall(MatDestroySubMatrices(1, &sub));
2094: }
2095: PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1]));
2096: PetscCall(PetscObjectDereference((PetscObject)uis));
2097: PetscCall(PetscObjectDereference((PetscObject)uaux));
2098: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), inner->mat->rmap->n, inner->mat->cmap->n, inner->mat->rmap->N, inner->mat->cmap->N, ctx, &S)); /* MatShell computing the action of M^-1 A or A M^-1 */
2099: PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_SchurCorrection));
2100: PetscCall(MatShellSetOperation(S, MATOP_VIEW, (PetscErrorCodeFn *)MatView_SchurCorrection));
2101: PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_SchurCorrection));
2102: PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx))));
2103: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) {
2104: PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx));
2105: } else { /* no support for PC_SYMMETRIC */
2106: PetscCheck(std::get<2>(*ctx) == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCSide %s (!= %s or %s or %s)", PCSides[std::get<2>(*ctx)], PCSides[PC_SIDE_DEFAULT], PCSides[PC_LEFT], PCSides[PC_RIGHT]);
2107: }
2108: PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx));
2109: PetscCall(PetscObjectContainerCompose((PetscObject)std::get<0>(*ctx)[1], "_PCHPDDM_Schur", ctx, nullptr));
2110: PetscCall(PCSetUp(std::get<0>(*ctx)[1]));
2111: PetscCall(KSPSetOperators(inner_ksp, S, S));
2112: PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1));
2113: PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2));
2114: PetscCall(PetscObjectDereference((PetscObject)inner_ksp));
2115: PetscCall(PetscObjectDereference((PetscObject)S));
2116: } else {
2117: std::get<0>(*ctx)[0] = pc_00;
2118: PetscCall(PetscObjectContainerCompose((PetscObject)pc, "_PCHPDDM_Schur", ctx, nullptr));
2119: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data_00->is), P->rmap->n, P->rmap->rstart, 1, &data->is)); /* dummy variables in the case of a centralized Schur complement */
2120: if (A11) {
2121: PetscCall(MatGetDiagonalBlock(A11, &data->aux));
2122: PetscCall(PetscObjectReference((PetscObject)data->aux));
2123: } else {
2124: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, P->rmap->n, P->rmap->n, 0, nullptr, &data->aux));
2125: PetscCall(MatAssemblyBegin(data->aux, MAT_FINAL_ASSEMBLY));
2126: PetscCall(MatAssemblyEnd(data->aux, MAT_FINAL_ASSEMBLY));
2127: }
2128: PetscCall(PCSetUp(pc));
2129: }
2130: for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it));
2131: PetscFunctionReturn(PETSC_SUCCESS);
2132: } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */
2133: PetscCall(PetscContainerGetPointer(container, static_cast<void *>(&ctx)));
2134: }
2135: }
2136: }
2137: }
2138: if (!data->is && data->N > 1) {
2139: char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */
2141: PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, ""));
2142: if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) {
2143: Mat B;
2145: PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre));
2146: if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED;
2147: PetscCall(MatDestroy(&B));
2148: } else {
2149: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg));
2150: if (flg) {
2151: Mat A00, P00, A01, A10, A11, B, N;
2152: PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES;
2154: PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11));
2155: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg));
2156: if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) {
2157: Mat B01;
2158: Vec diagonal = nullptr;
2159: const PetscScalar *array;
2160: MatSchurComplementAinvType type;
2162: PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10, &B01));
2163: if (A11) {
2164: PetscCall(MatCreateVecs(A11, &diagonal, nullptr));
2165: PetscCall(MatGetDiagonal(A11, diagonal));
2166: }
2167: PetscCall(MatCreateVecs(P00, &v, nullptr));
2168: PetscCall(MatSchurComplementGetAinvType(P, &type));
2169: PetscCheck(type == MAT_SCHUR_COMPLEMENT_AINV_DIAG || type == MAT_SCHUR_COMPLEMENT_AINV_LUMP || type == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "-%smat_schur_complement_ainv_type %s",
2170: ((PetscObject)P)->prefix ? ((PetscObject)P)->prefix : "", MatSchurComplementAinvTypes[type]);
2171: if (type != MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
2172: if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
2173: PetscCall(MatGetRowSum(P00, v));
2174: if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00));
2175: PetscCall(MatDestroy(&P00));
2176: PetscCall(VecGetArrayRead(v, &array));
2177: PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00));
2178: PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2179: for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES));
2180: PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY));
2181: PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY));
2182: PetscCall(VecRestoreArrayRead(v, &array));
2183: PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */
2184: PetscCall(MatDestroy(&P00));
2185: } else PetscCall(MatGetDiagonal(P00, v));
2186: PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */
2187: PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */
2188: PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B));
2189: PetscCall(MatDiagonalScale(B, v, nullptr));
2190: if (B01) PetscCall(MatDiagonalScale(B01, v, nullptr));
2191: } else {
2192: Mat D00;
2193: MatType type;
2195: PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), &D00));
2196: PetscCall(MatSetType(D00, MATAIJ));
2197: PetscCall(MatSetOptionsPrefix(D00, ((PetscObject)A00)->prefix));
2198: PetscCall(MatAppendOptionsPrefix(D00, "block_diagonal_"));
2199: PetscCall(MatSetFromOptions(D00)); /* for setting -mat_block_size dynamically */
2200: PetscCall(MatConvert(A00, MATAIJ, MAT_INITIAL_MATRIX, &B)); /* not all MatTypes have a MatInvertBlockDiagonal() implementation, plus one may want to use a different block size than the one of A00 */
2201: PetscCall(MatSetBlockSizesFromMats(B, D00, D00));
2202: PetscCall(MatInvertBlockDiagonalMat(B, D00));
2203: PetscCall(MatDestroy(&B));
2204: PetscCall(MatGetType(A01, &type)); /* cache MatType */
2205: PetscCall(MatConvert(A01, MATAIJ, MAT_INPLACE_MATRIX, &A01)); /* MatProduct is not versatile enough to fallback gracefully if no implementation found, so MatConvert() */
2206: PetscCall(MatMatMult(D00, A01, MAT_INITIAL_MATRIX, PETSC_CURRENT, &B));
2207: PetscCall(MatDestroy(&D00));
2208: PetscCall(MatConvert(A01, type, MAT_INPLACE_MATRIX, &A01)); /* reset to previous MatType */
2209: PetscCall(MatConvert(B, type, MAT_INPLACE_MATRIX, &B));
2210: if (!B01) { /* symmetric case */
2211: B01 = A01;
2212: PetscCall(PetscObjectReference((PetscObject)B01));
2213: }
2214: }
2215: if (B01 && B01 != A01) PetscCall(MatSetBlockSizesFromMats(B01, A01, A01)); /* TODO: remove this line once Firedrake is fixed */
2216: PetscCall(VecDestroy(&v));
2217: PetscCall(MatCreateNormalHermitian(B, &N));
2218: PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal, B01));
2219: PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg));
2220: if (!flg) {
2221: PetscCall(MatDestroy(&P));
2222: P = N;
2223: PetscCall(PetscObjectReference((PetscObject)P));
2224: }
2225: if (diagonal) {
2226: PetscCall(MatSetOption(P, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); /* may have missing diagonal entries */
2227: PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES));
2228: PetscCall(PCSetOperators(pc, P, P)); /* replace P by A01^T inv(diag(P00)) A01 - diag(P11) */
2229: PetscCall(VecDestroy(&diagonal));
2230: } else PetscCall(PCSetOperators(pc, B01 ? P : N, P)); /* replace P by A01^T inv(diag(P00)) A01 */
2231: pc->ops->postsolve = PCPostSolve_SchurPreLeastSquares; /* PCFIELDSPLIT expect a KSP for (P11 - A10 inv(diag(P00)) A01) */
2232: PetscCall(MatDestroy(&N)); /* but a PC for (A10 inv(diag(P00)) A10 - P11) is setup instead */
2233: PetscCall(MatDestroy(&P)); /* so the sign of the solution must be flipped */
2234: PetscCall(MatDestroy(&B));
2235: } else
2236: PetscCheck(type != PC_HPDDM_SCHUR_PRE_GENEO, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s without a prior call to PCHPDDMSetAuxiliaryMat() on the A11 block%s%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], pcpre ? " -" : "", pcpre ? pcpre : "");
2237: for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it));
2238: PetscFunctionReturn(PETSC_SUCCESS);
2239: } else {
2240: PetscCall(PetscOptionsGetString(((PetscObject)pc)->options, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr));
2241: PetscCall(PetscStrcmp(type, PCMAT, &algebraic));
2242: PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_levels_1_st_pc_type mat and -%spc_hpddm_block_splitting", pcpre ? pcpre : "", pcpre ? pcpre : "");
2243: if (overlap != -1) {
2244: PetscCheck(!block && !algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_%s and -%spc_hpddm_harmonic_overlap", pcpre ? pcpre : "", block ? "block_splitting" : "levels_1_st_pc_type mat", pcpre ? pcpre : "");
2245: PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-%spc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", pcpre ? pcpre : "", overlap);
2246: }
2247: if (block || overlap != -1) algebraic = PETSC_TRUE;
2248: if (algebraic) {
2249: PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is));
2250: PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1));
2251: PetscCall(ISSort(data->is));
2252: } else
2253: PetscCall(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 and -%spc_hpddm_harmonic_overlap < 1\n", pcpre ? pcpre : "", pcpre ? pcpre : "", pcpre ? pcpre : ""));
2254: }
2255: }
2256: }
2257: }
2258: if (PetscDefined(USE_DEBUG)) {
2259: if (data->is) PetscCall(ISDuplicate(data->is, &dis));
2260: if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux));
2261: }
2262: if (data->is || (ismatis && data->N > 1)) {
2263: if (ismatis) {
2264: PetscCall(MatISGetLocalMat(P, &N));
2265: PetscCall(PetscObjectTypeCompareAny((PetscObject)N, &flg, MATSEQBAIJ, MATSEQSBAIJ, ""));
2266: PetscCall(MatISRestoreLocalMat(P, &N));
2267: PetscCall(MatConvert(P, flg ? MATMPIBAIJ : MATMPIAIJ, MAT_INITIAL_MATRIX, &C));
2268: PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr));
2269: PetscCall(PetscObjectReference((PetscObject)P));
2270: PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C));
2271: std::swap(C, P);
2272: PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n));
2273: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc));
2274: PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0]));
2275: PetscCall(ISDestroy(&loc));
2276: /* the auxiliary Mat is _not_ the local Neumann matrix */
2277: /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */
2278: data->Neumann = PETSC_BOOL3_FALSE;
2279: structure = SAME_NONZERO_PATTERN;
2280: } else {
2281: is[0] = data->is;
2282: if (algebraic || ctx) subdomains = PETSC_TRUE;
2283: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr));
2284: if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre);
2285: if (PetscBool3ToBool(data->Neumann)) {
2286: PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_block_splitting and -%spc_hpddm_has_neumann", pcpre ? pcpre : "", pcpre ? pcpre : "");
2287: PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_harmonic_overlap %" PetscInt_FMT " and -%spc_hpddm_has_neumann", pcpre ? pcpre : "", overlap, pcpre ? pcpre : "");
2288: PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_levels_1_st_pc_type mat and -%spc_hpddm_has_neumann", pcpre ? pcpre : "", pcpre ? pcpre : "");
2289: }
2290: if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN;
2291: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc));
2292: }
2293: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2294: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */
2295: if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */
2296: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : ""));
2297: PetscCall(PetscOptionsSetValue(((PetscObject)pc)->options, prefix, MatStructures[structure]));
2298: }
2299: flg = PETSC_FALSE;
2300: if (data->share) {
2301: data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */
2302: if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : ""));
2303: else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n"));
2304: else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n"));
2305: else if (!algebraic && structure != SAME_NONZERO_PATTERN)
2306: PetscCall(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]));
2307: else {
2308: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHTOOL, &flg));
2309: if (flg) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since Pmat is of type MATHTOOL\n"));
2310: else data->share = PETSC_TRUE;
2311: }
2312: if (!data->share) {
2313: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_share_sub_ksp", pcpre ? pcpre : ""));
2314: PetscCall(PetscOptionsClearValue(((PetscObject)pc)->options, prefix));
2315: }
2316: }
2317: if (!ismatis) {
2318: if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted));
2319: else unsorted = is[0];
2320: }
2321: if ((ctx || data->N > 1) && (data->aux || ismatis || algebraic)) {
2322: PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level");
2323: PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2324: if (ismatis) {
2325: /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */
2326: PetscCall(MatIncreaseOverlap(P, 1, is, 1));
2327: PetscCall(ISDestroy(&data->is));
2328: data->is = is[0];
2329: } else {
2330: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_TRUE));
2331: if (!ctx && overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private));
2332: }
2333: if (algebraic && overlap == -1) {
2334: PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block));
2335: if (block) {
2336: PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux));
2337: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr));
2338: }
2339: } else if (!ctx) {
2340: if (PetscBool3ToBool(data->Neumann)) sub = &data->aux;
2341: else {
2342: PetscBool flg;
2344: if (overlap != -1) {
2345: Harmonic h;
2346: Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */
2347: IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */
2348: const PetscInt *i[2], bs = P->cmap->bs; /* with a GEVP: [ A_00 A_01 ] */
2349: PetscInt n[2], location; /* [ A_10 A_11 A_12 ] */
2350: std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */
2352: do {
2353: PetscCall(ISDuplicate(data->is, ov));
2354: if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1));
2355: PetscCall(ISDuplicate(ov[0], ov + 1));
2356: PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1));
2357: PetscCall(ISGetLocalSize(ov[0], n));
2358: PetscCall(ISGetLocalSize(ov[1], n + 1));
2359: flg = PetscBool(n[0] == n[1] && n[0] != P->rmap->n);
2360: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
2361: if (flg) {
2362: PetscCall(ISDestroy(ov));
2363: PetscCall(ISDestroy(ov + 1));
2364: PetscCheck(--overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No oversampling possible");
2365: PetscCall(PetscInfo(pc, "Supplied -%spc_hpddm_harmonic_overlap parameter is too large, it has been decreased to %" PetscInt_FMT "\n", pcpre ? pcpre : "", overlap));
2366: } else break;
2367: } while (1);
2368: PetscCall(PetscNew(&h));
2369: h->ksp = nullptr;
2370: PetscCall(PetscCalloc1(2, &h->A));
2371: PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-eps_nev", &flg));
2372: if (!flg) {
2373: PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-svd_nsv", &flg));
2374: if (!flg) PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-svd_threshold_relative", &flg));
2375: } else flg = PETSC_FALSE;
2376: PetscCall(ISSort(ov[0]));
2377: if (!flg) PetscCall(ISSort(ov[1]));
2378: PetscCall(PetscCalloc1(5, &h->is));
2379: PetscCall(MatCreateSubMatrices(P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */
2380: for (PetscInt j = 0; j < 2; ++j) PetscCall(ISGetIndices(ov[j], i + j));
2381: v[1].reserve((n[1] - n[0]) / bs);
2382: for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */
2383: PetscCall(ISLocate(ov[0], i[1][j], &location));
2384: if (location < 0) v[1].emplace_back(j / bs);
2385: }
2386: if (!flg) {
2387: h->A[1] = a[0];
2388: PetscCall(PetscObjectReference((PetscObject)h->A[1]));
2389: v[0].reserve((n[0] - P->rmap->n) / bs);
2390: for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */
2391: PetscCall(ISLocate(loc, i[1][j], &location));
2392: if (location < 0) {
2393: PetscCall(ISLocate(ov[0], i[1][j], &location));
2394: if (location >= 0) v[0].emplace_back(j / bs);
2395: }
2396: }
2397: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows));
2398: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4));
2399: PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */
2400: PetscCall(ISDestroy(&rows));
2401: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows));
2402: PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */
2403: PetscCall(ISDestroy(&rows));
2404: v[0].clear();
2405: PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3));
2406: PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2));
2407: }
2408: v[0].reserve((n[0] - P->rmap->n) / bs);
2409: for (PetscInt j = 0; j < n[0]; j += bs) {
2410: PetscCall(ISLocate(loc, i[0][j], &location));
2411: if (location < 0) v[0].emplace_back(j / bs);
2412: }
2413: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows));
2414: for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j));
2415: if (flg) {
2416: PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &stride));
2417: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols));
2418: PetscCall(MatCreateSubMatrix(a[0], stride, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */
2419: PetscCall(ISDestroy(&cols));
2420: PetscCall(ISDestroy(&stride));
2421: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg));
2422: if (flg) { /* initial Pmat was MATSBAIJ, convert back to the same format since this submatrix is square */
2423: PetscCall(MatSetOption(A0, MAT_SYMMETRIC, PETSC_TRUE));
2424: PetscCall(MatConvert(A0, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &A0));
2425: }
2426: flg = PETSC_TRUE;
2427: PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2));
2428: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols));
2429: PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */
2430: PetscCall(ISDestroy(&cols));
2431: }
2432: PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride));
2433: PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is));
2434: PetscCall(ISDestroy(&stride));
2435: PetscCall(ISDestroy(&rows));
2436: PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1));
2437: if (subdomains) {
2438: if (!data->levels[0]->pc) {
2439: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc));
2440: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2441: PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix));
2442: PetscCall(PCSetOperators(data->levels[0]->pc, A, P));
2443: }
2444: PetscCall(PCSetType(data->levels[0]->pc, PCASM));
2445: if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc));
2446: PetscCall(PCSetModifySubMatrices(data->levels[0]->pc, pc->modifysubmatrices, pc->modifysubmatricesP));
2447: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE));
2448: if (!flg) ++overlap;
2449: if (data->share) {
2450: PetscInt n = -1;
2451: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp));
2452: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2453: if (flg) {
2454: h->ksp = ksp[0];
2455: PetscCall(PetscObjectReference((PetscObject)h->ksp));
2456: }
2457: }
2458: }
2459: if (!h->ksp) {
2460: PetscBool share = data->share;
2462: PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp));
2463: PetscCall(KSPSetType(h->ksp, KSPPREONLY));
2464: PetscCall(KSPSetOperators(h->ksp, A0, A0));
2465: do {
2466: if (!data->share) {
2467: share = PETSC_FALSE;
2468: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_"));
2469: PetscCall(KSPSetOptionsPrefix(h->ksp, prefix));
2470: PetscCall(KSPSetFromOptions(h->ksp));
2471: } else {
2472: MatSolverType type;
2474: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp[0]->pc, &data->share, PCLU, PCCHOLESKY, ""));
2475: if (data->share) {
2476: PetscCall(PCFactorGetMatSolverType(ksp[0]->pc, &type));
2477: if (!type) {
2478: if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(ksp[0]->pc, MATSOLVERMUMPS));
2479: else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(ksp[0]->pc, MATSOLVERMKL_PARDISO));
2480: else data->share = PETSC_FALSE;
2481: if (data->share) PetscCall(PCSetFromOptions(ksp[0]->pc));
2482: } else {
2483: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share));
2484: if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share));
2485: }
2486: if (data->share) {
2487: std::tuple<KSP, IS, Vec[2]> *p;
2489: PetscCall(PCFactorGetMatrix(ksp[0]->pc, &A));
2490: PetscCall(MatFactorSetSchurIS(A, h->is[4]));
2491: PetscCall(KSPSetUp(ksp[0]));
2492: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : ""));
2493: PetscCall(KSPSetOptionsPrefix(h->ksp, prefix));
2494: PetscCall(KSPSetFromOptions(h->ksp));
2495: PetscCall(PCSetType(h->ksp->pc, PCSHELL));
2496: PetscCall(PetscNew(&p));
2497: std::get<0>(*p) = ksp[0];
2498: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p)));
2499: PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1));
2500: PetscCall(PCShellSetContext(h->ksp->pc, p));
2501: PetscCall(PCShellSetApply(h->ksp->pc, PCApply_Schur));
2502: PetscCall(PCShellSetApplyTranspose(h->ksp->pc, PCApply_Schur<Vec, true>));
2503: PetscCall(PCShellSetMatApply(h->ksp->pc, PCApply_Schur<Mat>));
2504: PetscCall(PCShellSetDestroy(h->ksp->pc, PCDestroy_Schur));
2505: }
2506: }
2507: if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n"));
2508: }
2509: } while (!share != !data->share); /* if data->share is initially PETSC_TRUE, but then reset to PETSC_FALSE, then go back to the beginning of the do loop */
2510: }
2511: PetscCall(ISDestroy(ov));
2512: PetscCall(ISDestroy(ov + 1));
2513: if (overlap == 1 && subdomains && flg) {
2514: *subA = A0;
2515: sub = subA;
2516: } else PetscCall(MatDestroy(&A0));
2517: PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux));
2518: PetscCall(MatSetVecType(data->aux, h->A[0]->defaultvectype));
2519: PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */
2520: PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr));
2521: PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Harmonic));
2522: PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Harmonic));
2523: PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE));
2524: PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AtB, nullptr, MatProduct_AtB_Harmonic, nullptr, MATDENSE, MATDENSE));
2525: PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Harmonic));
2526: PetscCall(MatDestroySubMatrices(1, &a));
2527: }
2528: if (overlap != 1 || !subdomains) {
2529: PetscCall(MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub));
2530: if (ismatis) {
2531: PetscCall(MatISGetLocalMat(C, &N));
2532: PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg));
2533: if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub));
2534: PetscCall(MatISRestoreLocalMat(C, &N));
2535: }
2536: }
2537: }
2538: }
2539: if (data->N > 1) {
2540: /* Vec holding the partition of unity */
2541: if (!data->levels[0]->D) {
2542: PetscCall(ISGetLocalSize(data->is, &n));
2543: PetscCall(VecCreate(PETSC_COMM_SELF, &data->levels[0]->D));
2544: PetscCall(VecSetSizes(data->levels[0]->D, n, n));
2545: PetscCall(VecSetType(data->levels[0]->D, A->defaultvectype));
2546: }
2547: if (data->share && overlap == -1) {
2548: Mat D;
2549: IS perm = nullptr;
2550: PetscInt size = -1;
2552: if (!data->levels[0]->pc) {
2553: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2554: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc));
2555: PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix));
2556: PetscCall(PCSetOperators(data->levels[0]->pc, A, P));
2557: }
2558: PetscCall(PCSetType(data->levels[0]->pc, PCASM));
2559: if (!ctx) {
2560: if (!data->levels[0]->pc->setupcalled) {
2561: IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */
2563: PetscCall(ISDuplicate(is[0], &sorted));
2564: PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc));
2565: PetscCall(PetscObjectDereference((PetscObject)sorted));
2566: }
2567: PetscCall(PCSetFromOptions(data->levels[0]->pc));
2568: PetscCall(PCSetModifySubMatrices(data->levels[0]->pc, pc->modifysubmatrices, pc->modifysubmatricesP));
2569: if (block) {
2570: PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm));
2571: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic));
2572: } else PetscCall(PCSetUp(data->levels[0]->pc));
2573: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp));
2574: if (size != 1) {
2575: data->share = PETSC_FALSE;
2576: PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size);
2577: PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n"));
2578: PetscCall(ISDestroy(&unsorted));
2579: unsorted = is[0];
2580: } else {
2581: const char *matpre;
2582: PetscBool cmp[4];
2584: if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm));
2585: if (perm) { /* unsorted input IS */
2586: if (!PetscBool3ToBool(data->Neumann) && !block) {
2587: PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */
2588: PetscCall(MatHeaderReplace(sub[0], &D));
2589: }
2590: if (data->B) { /* see PCHPDDMSetRHSMat() */
2591: PetscCall(MatPermute(data->B, perm, perm, &D));
2592: PetscCall(MatHeaderReplace(data->B, &D));
2593: }
2594: PetscCall(ISDestroy(&perm));
2595: }
2596: PetscCall(KSPGetOperators(ksp[0], subA, subA + 1));
2597: PetscCall(PetscObjectReference((PetscObject)subA[0]));
2598: PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D));
2599: PetscCall(MatGetOptionsPrefix(subA[1], &matpre));
2600: PetscCall(MatSetOptionsPrefix(D, matpre));
2601: PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp));
2602: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1));
2603: if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2));
2604: else cmp[2] = PETSC_FALSE;
2605: if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3));
2606: else cmp[3] = PETSC_FALSE;
2607: PetscCheck(cmp[0] == cmp[1] && cmp[2] == cmp[3], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_levels_1_pc_asm_sub_mat_type %s and auxiliary Mat of type %s", pcpre ? pcpre : "", ((PetscObject)D)->type_name, ((PetscObject)C)->type_name);
2608: if (!cmp[0] && !cmp[2]) {
2609: if (!block) {
2610: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckMatStructure_Private(pc, D, C));
2611: PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN));
2612: } else {
2613: structure = DIFFERENT_NONZERO_PATTERN;
2614: PetscCall(MatAXPY(D, 1.0, data->aux, structure));
2615: }
2616: } else {
2617: Mat mat[2];
2619: if (cmp[0]) {
2620: PetscCall(MatNormalGetMat(D, mat));
2621: PetscCall(MatNormalGetMat(C, mat + 1));
2622: } else {
2623: PetscCall(MatNormalHermitianGetMat(D, mat));
2624: PetscCall(MatNormalHermitianGetMat(C, mat + 1));
2625: }
2626: PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN));
2627: }
2628: PetscCall(MatPropagateSymmetryOptions(C, D));
2629: PetscCall(MatDestroy(&C));
2630: C = D;
2631: /* swap pointers so that variables stay consistent throughout PCSetUp() */
2632: std::swap(C, data->aux);
2633: std::swap(uis, data->is);
2634: swap = PETSC_TRUE;
2635: }
2636: }
2637: }
2638: }
2639: if (ctx) {
2640: PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data;
2641: PC s;
2642: Mat A00, P00, A01 = nullptr, A10, N, b[4];
2643: IS sorted, is[2], *is_00;
2644: MatSolverType type;
2645: std::pair<PC, Vec[2]> *p;
2647: n = -1;
2648: PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp));
2649: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2650: PetscCall(KSPGetOperators(ksp[0], subA, subA + 1));
2651: PetscCall(ISGetLocalSize(data_00->is, &n));
2652: if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) {
2653: PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr));
2654: PetscCall(ISGetLocalSize(*is_00, &n));
2655: PetscCheck(n == subA[0]->rmap->n && n == subA[0]->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre ? pcpre : "", ((PetscObject)pc)->prefix);
2656: } else is_00 = &data_00->is;
2657: PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */
2658: std::swap(C, data->aux);
2659: std::swap(uis, data->is);
2660: swap = PETSC_TRUE;
2661: PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, nullptr));
2662: std::get<1>(*ctx)[1] = A10;
2663: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
2664: if (flg) PetscCall(MatTransposeGetMat(A10, &A01));
2665: else {
2666: PetscBool flg;
2668: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
2669: if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01));
2670: }
2671: PetscCall(ISDuplicate(*is_00, &sorted)); /* during setup of the PC associated to the A00 block, this IS has already been sorted, but it's put back to its original state at the end of PCSetUp_HPDDM(), which may be unsorted */
2672: PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */
2673: if (!A01) {
2674: PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2675: PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub));
2676: b[2] = sub[0];
2677: PetscCall(PetscObjectReference((PetscObject)sub[0]));
2678: PetscCall(MatDestroySubMatrices(1, &sub));
2679: PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg));
2680: A10 = nullptr;
2681: if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10));
2682: else {
2683: PetscBool flg;
2685: PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg));
2686: if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10));
2687: }
2688: if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub));
2689: else {
2690: if (flg) PetscCall(MatCreateTranspose(b[2], b + 1));
2691: else PetscCall(MatCreateHermitianTranspose(b[2], b + 1));
2692: }
2693: } else {
2694: PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2695: PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub));
2696: if (flg) PetscCall(MatCreateTranspose(*sub, b + 2));
2697: else PetscCall(MatCreateHermitianTranspose(*sub, b + 2));
2698: }
2699: if (A01 || !A10) {
2700: b[1] = sub[0];
2701: PetscCall(PetscObjectReference((PetscObject)sub[0]));
2702: }
2703: PetscCall(MatDestroySubMatrices(1, &sub));
2704: PetscCall(ISDestroy(&sorted));
2705: b[3] = data->aux;
2706: PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S));
2707: PetscCall(MatSchurComplementSetKSP(S, ksp[0]));
2708: if (data->N != 1) {
2709: PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */
2710: PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc));
2711: PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */
2712: s = data->levels[0]->pc;
2713: } else {
2714: is[0] = data->is;
2715: PetscCall(PetscObjectReference((PetscObject)is[0]));
2716: PetscCall(PetscObjectReference((PetscObject)b[3]));
2717: PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */
2718: data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */
2719: PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */
2720: PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc));
2721: PetscCall(ISDestroy(is));
2722: PetscCall(ISDestroy(&loc));
2723: s = pc;
2724: }
2725: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */
2726: PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp));
2727: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2728: PetscCall(KSPGetPC(ksp[0], &inner));
2729: PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */
2730: b[0] = subA[0];
2731: PetscCall(MatCreateNest(PETSC_COMM_SELF, 2, nullptr, 2, nullptr, b, &N)); /* instead of computing inv(A11 - A10 inv(A00) A01), compute inv([A00, A01; A10, A11]) followed by a partial solution associated to the A11 block */
2732: if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3]));
2733: PetscCall(PetscObjectDereference((PetscObject)b[1]));
2734: PetscCall(PetscObjectDereference((PetscObject)b[2]));
2735: PetscCall(PCCreate(PETSC_COMM_SELF, &s));
2736: PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix));
2737: PetscCall(PCSetOptionsPrefix(inner, nullptr));
2738: PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE));
2739: PetscCall(PCSetType(s, PCLU));
2740: if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(s, MATSOLVERMUMPS)); /* only MATSOLVERMUMPS handles MATNEST, so for the others, e.g., MATSOLVERPETSC or MATSOLVERMKL_PARDISO, convert to plain MATAIJ */
2741: PetscCall(PCSetFromOptions(s));
2742: PetscCall(PCFactorGetMatSolverType(s, &type));
2743: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
2744: PetscCall(MatGetLocalSize(P, &n, nullptr));
2745: if (flg || n == 0) {
2746: PetscCall(PCSetOperators(s, N, N));
2747: if (n) {
2748: PetscCall(PCFactorGetMatrix(s, b));
2749: PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix));
2750: n = -1;
2751: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr));
2752: if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */
2753: PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */
2754: PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */
2755: }
2756: } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */
2757: } else {
2758: MatInfo info;
2760: PetscCall(MatGetInfo(b[3], MAT_LOCAL, &info));
2761: if (info.nz_used == 0.0) {
2762: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, b[3]->rmap->n, b[3]->cmap->n, 1, nullptr, b + 3));
2763: for (PetscInt i = 0; i < b[3]->rmap->n; ++i) PetscCall(MatSetValue(b[3], i, i, 0.0, INSERT_VALUES));
2764: PetscCall(MatAssemblyBegin(b[3], MAT_FINAL_ASSEMBLY));
2765: PetscCall(MatAssemblyEnd(b[3], MAT_FINAL_ASSEMBLY));
2766: PetscCall(MatNestSetSubMat(N, 1, 1, b[3]));
2767: PetscCall(PetscObjectDereference((PetscObject)b[3]));
2768: }
2769: PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b));
2770: PetscCall(PCSetOperators(s, N, *b));
2771: PetscCall(PetscObjectDereference((PetscObject)*b));
2772: PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, ""));
2773: if (flg) {
2774: PetscCall(PCFactorGetMatrix(s, b)); /* MATSOLVERMKL_PARDISO cannot compute in PETSc (yet) a partial solution associated to the A11 block, only partial solution associated to the A00 block or full solution */
2775: if (info.nz_used == 0.0) {
2776: PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, ""));
2777: if (flg) {
2778: PetscCall(MatFactorGetSolverType(*b, &type));
2779: PetscCall(PetscStrcmp(type, MATSOLVERPETSC, &flg));
2780: if (flg) {
2781: PetscCall(PetscOptionsHasName(((PetscObject)s)->options, ((PetscObject)s)->prefix, "-pc_factor_mat_ordering_type", &flg));
2782: if (!flg) PetscCall(PCFactorSetMatOrderingType(s, MATORDERINGNATURAL));
2783: }
2784: }
2785: flg = PETSC_TRUE;
2786: }
2787: }
2788: }
2789: PetscCall(PetscNew(&p));
2790: p->first = s;
2791: if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1));
2792: else p->second[0] = p->second[1] = nullptr;
2793: PetscCall(PCShellSetContext(inner, p));
2794: PetscCall(PCShellSetApply(inner, PCApply_Nest));
2795: PetscCall(PCShellSetView(inner, PCView_Nest));
2796: PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest));
2797: PetscCall(PetscObjectDereference((PetscObject)N));
2798: if (!data) {
2799: PetscCall(MatDestroy(&S));
2800: PetscCall(ISDestroy(&unsorted));
2801: PetscCall(MatDestroy(&C));
2802: PetscCall(ISDestroy(&uis));
2803: PetscCall(PetscFree(ctx));
2804: if (PetscDefined(USE_DEBUG)) {
2805: PetscCall(ISDestroy(&dis));
2806: PetscCall(MatDestroy(&daux));
2807: }
2808: PetscFunctionReturn(PETSC_SUCCESS);
2809: }
2810: }
2811: if (!data->levels[0]->scatter) {
2812: PetscCall(MatCreateVecs(P, &xin, nullptr));
2813: if (ismatis) PetscCall(MatDestroy(&P));
2814: PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter));
2815: PetscCall(VecDestroy(&xin));
2816: }
2817: if (data->levels[0]->P) {
2818: /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */
2819: PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], !pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE));
2820: }
2821: if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>();
2822: if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr));
2823: else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr));
2824: /* HPDDM internal data structure */
2825: PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels));
2826: if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr));
2827: /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */
2828: if (!ctx) {
2829: if (data->deflation || overlap != -1) weighted = data->aux;
2830: else if (!data->B) {
2831: PetscBool cmp;
2833: PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted));
2834: PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, ""));
2835: if (cmp) flg = PETSC_FALSE;
2836: PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D));
2837: /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */
2838: /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */
2839: PetscCall(MatPropagateSymmetryOptions(sub[0], weighted));
2840: if (PetscDefined(USE_DEBUG) && PetscBool3ToBool(data->Neumann)) {
2841: Mat *sub, A[2];
2842: PetscReal norm[2];
2844: PetscCall(MatCreateSubMatrices(P, 1, &data->is, &data->is, MAT_INITIAL_MATRIX, &sub));
2845: PetscCall(MatDiagonalScale(sub[0], data->levels[0]->D, data->levels[0]->D));
2846: PetscCall(MatConvert(sub[0], MATSEQAIJ, MAT_INITIAL_MATRIX, A)); /* too many corner cases to handle (MATNORMAL, MATNORMALHERMITIAN, MATBAIJ with different block sizes...), so just MatConvert() to MATSEQAIJ since this is just for debugging */
2847: PetscCall(MatConvert(weighted, MATSEQAIJ, MAT_INITIAL_MATRIX, A + 1));
2848: PetscCall(MatAXPY(A[0], -1.0, A[1], UNKNOWN_NONZERO_PATTERN));
2849: PetscCall(MatNorm(A[0], NORM_FROBENIUS, norm));
2850: if (norm[0]) {
2851: PetscCall(MatNorm(A[1], NORM_FROBENIUS, norm + 1));
2852: PetscCheck(PetscAbsReal(norm[0] / norm[1]) < PetscSqrtReal(PETSC_SMALL), PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Auxiliary Mat is different from the (assembled) subdomain Mat for the interior unknowns, so it cannot be the Neumann matrix, remove -%spc_hpddm_has_neumann", pcpre ? pcpre : "");
2853: }
2854: PetscCall(MatDestroySubMatrices(1, &sub));
2855: for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(A + i));
2856: }
2857: } else weighted = data->B;
2858: } else weighted = nullptr;
2859: /* SLEPc is used inside the loaded symbol */
2860: PetscCall((*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : (algebraic && !block && overlap == -1 ? sub[0] : (!ctx ? data->aux : S)), weighted, data->B, initial, data->levels));
2861: if (!ctx && data->share && overlap == -1) {
2862: Mat st[2];
2864: PetscCheck(ksp && ksp[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing shared subdomain KSP");
2865: PetscCall(KSPGetOperators(ksp[0], st, st + 1));
2866: PetscCall(MatCopy(subA[0], st[0], structure));
2867: if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN));
2868: PetscCall(PetscObjectDereference((PetscObject)subA[0]));
2869: }
2870: if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr));
2871: if (ismatis) PetscCall(MatISGetLocalMat(C, &N));
2872: else N = data->aux;
2873: if (!ctx) P = sub[0];
2874: else P = S;
2875: /* going through the grid hierarchy */
2876: for (n = 1; n < data->N; ++n) {
2877: if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr));
2878: /* method composed in the loaded symbol since there, SLEPc is used as well */
2879: PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels));
2880: if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr));
2881: }
2882: /* reset to NULL to avoid any faulty use */
2883: PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr));
2884: if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr));
2885: else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */
2886: for (n = 0; n < data->N - 1; ++n)
2887: if (data->levels[n]->P) {
2888: /* HPDDM internal work buffers */
2889: PetscCallCXX(data->levels[n]->P->setBuffer());
2890: PetscCallCXX(data->levels[n]->P->super::start());
2891: }
2892: if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub));
2893: if (ismatis) data->is = nullptr;
2894: for (n = 0; n < data->N - 1 + (reused > 0); ++n) {
2895: if (data->levels[n]->P) {
2896: PC spc;
2898: /* force the PC to be PCSHELL to do the coarse grid corrections */
2899: PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE));
2900: PetscCall(KSPGetPC(data->levels[n]->ksp, &spc));
2901: PetscCall(PCSetType(spc, PCSHELL));
2902: PetscCall(PCShellSetContext(spc, data->levels[n]));
2903: PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell));
2904: PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell));
2905: PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell));
2906: PetscCall(PCShellSetApplyTranspose(spc, PCApplyTranspose_HPDDMShell));
2907: PetscCall(PCShellSetMatApplyTranspose(spc, PCMatApplyTranspose_HPDDMShell));
2908: if (ctx && n == 0) {
2909: Mat Amat, Pmat;
2910: PetscInt m, M;
2911: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
2913: PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat));
2914: PetscCall(MatGetLocalSize(Pmat, &m, nullptr));
2915: PetscCall(MatGetSize(Pmat, &M, nullptr));
2916: PetscCall(PetscNew(&ctx));
2917: std::get<0>(*ctx) = S;
2918: std::get<1>(*ctx) = data->levels[n]->scatter;
2919: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat));
2920: PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Schur<false>));
2921: PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMult_Schur<true>));
2922: PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Schur));
2923: PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1));
2924: PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat));
2925: PetscCall(PetscObjectDereference((PetscObject)Amat));
2926: }
2927: PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell));
2928: if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc));
2929: if (n < reused) {
2930: PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE));
2931: PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE));
2932: }
2933: PetscCall(PCSetUp(spc));
2934: }
2935: }
2936: if (ctx) PetscCall(MatDestroy(&S));
2937: if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr));
2938: } else flg = reused ? PETSC_FALSE : PETSC_TRUE;
2939: if (!ismatis && subdomains) {
2940: if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner));
2941: else inner = data->levels[0]->pc;
2942: if (inner) {
2943: if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM));
2944: PetscCall(PCSetFromOptions(inner));
2945: PetscCall(PCSetModifySubMatrices(inner, pc->modifysubmatrices, pc->modifysubmatricesP));
2946: PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg));
2947: if (flg) {
2948: if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */
2949: IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */
2951: PetscCall(ISDuplicate(is[0], &sorted));
2952: PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc));
2953: PetscCall(PetscObjectDereference((PetscObject)sorted));
2954: }
2955: if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */
2956: PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr));
2957: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic));
2958: PetscCall(PetscObjectDereference((PetscObject)P));
2959: }
2960: }
2961: }
2962: if (data->N > 1) {
2963: if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub));
2964: if (overlap == 1) PetscCall(MatDestroy(subA));
2965: }
2966: }
2967: PetscCall(ISDestroy(&loc));
2968: } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */
2969: if (requested != data->N + reused) {
2970: 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,
2971: data->N, pcpre ? pcpre : ""));
2972: PetscCall(PetscInfo(pc, "It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%" PetscInt_FMT "_eps_threshold_absolute or a lower value for -%spc_hpddm_levels_%" PetscInt_FMT "_svd_threshold_relative, so that at least one local deflation vector will be selected\n", pcpre ? pcpre : "",
2973: data->N, pcpre ? pcpre : "", data->N));
2974: /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */
2975: for (n = data->N - 1; n < requested - 1; ++n) {
2976: if (data->levels[n]->P) {
2977: PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE));
2978: PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0]));
2979: PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1]));
2980: PetscCall(MatDestroy(data->levels[n]->V));
2981: PetscCall(MatDestroy(data->levels[n]->V + 1));
2982: PetscCall(MatDestroy(data->levels[n]->V + 2));
2983: PetscCall(VecDestroy(&data->levels[n]->D));
2984: PetscCall(PetscSFDestroy(&data->levels[n]->scatter));
2985: }
2986: }
2987: if (reused) {
2988: for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
2989: PetscCall(KSPDestroy(&data->levels[n]->ksp));
2990: PetscCall(PCDestroy(&data->levels[n]->pc));
2991: }
2992: }
2993: PetscCheck(!PetscDefined(USE_DEBUG), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "%" 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. It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%" PetscInt_FMT "_eps_threshold or a lower value for -%spc_hpddm_levels_%" PetscInt_FMT "_svd_threshold_relative, 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,
2994: data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N, pcpre ? pcpre : "", data->N);
2995: }
2996: /* these solvers are created after PCSetFromOptions() is called */
2997: if (pc->setfromoptionscalled) {
2998: for (n = 0; n < data->N; ++n) {
2999: if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp));
3000: if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc));
3001: }
3002: pc->setfromoptionscalled = 0;
3003: }
3004: data->N += reused;
3005: if (data->share && swap) {
3006: /* swap back pointers so that variables follow the user-provided numbering */
3007: std::swap(C, data->aux);
3008: std::swap(uis, data->is);
3009: PetscCall(MatDestroy(&C));
3010: PetscCall(ISDestroy(&uis));
3011: }
3012: if (algebraic) PetscCall(MatDestroy(&data->aux));
3013: if (unsorted && unsorted != is[0]) {
3014: PetscCall(ISCopy(unsorted, data->is));
3015: PetscCall(ISDestroy(&unsorted));
3016: }
3017: if (PetscDefined(USE_DEBUG)) {
3018: PetscCheck((data->is && dis) || (!data->is && !dis), PETSC_COMM_SELF, PETSC_ERR_PLIB, "An IS pointer is NULL but not the other: input IS pointer (%p) v. output IS pointer (%p)", (void *)dis, (void *)data->is);
3019: if (data->is) {
3020: PetscCall(ISEqualUnsorted(data->is, dis, &flg));
3021: PetscCall(ISDestroy(&dis));
3022: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal");
3023: }
3024: PetscCheck((data->aux && daux) || (!data->aux && !daux), PETSC_COMM_SELF, PETSC_ERR_PLIB, "A Mat pointer is NULL but not the other: input Mat pointer (%p) v. output Mat pointer (%p)", (void *)daux, (void *)data->aux);
3025: if (data->aux) {
3026: PetscCall(MatMultEqual(data->aux, daux, 20, &flg));
3027: PetscCall(MatDestroy(&daux));
3028: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal");
3029: }
3030: }
3031: PetscFunctionReturn(PETSC_SUCCESS);
3032: }
3034: /*@
3035: PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type.
3037: Collective
3039: Input Parameters:
3040: + pc - preconditioner context
3041: - type - coarse correction type, see `PCHPDDMCoarseCorrectionType`
3043: Options Database Key:
3044: . -pc_hpddm_coarse_correction (deflated|additive|balanced|none|deflated_reversed) - type of coarse correction to apply
3046: Level: intermediate
3048: .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
3049: @*/
3050: PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type)
3051: {
3052: PetscFunctionBegin;
3055: PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type));
3056: PetscFunctionReturn(PETSC_SUCCESS);
3057: }
3059: /*@
3060: PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type.
3062: Input Parameter:
3063: . pc - preconditioner context
3065: Output Parameter:
3066: . type - coarse correction type, see `PCHPDDMCoarseCorrectionType`
3068: Level: intermediate
3070: .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
3071: @*/
3072: PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type)
3073: {
3074: PetscFunctionBegin;
3076: if (type) {
3077: PetscAssertPointer(type, 2);
3078: PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type));
3079: }
3080: PetscFunctionReturn(PETSC_SUCCESS);
3081: }
3083: static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type)
3084: {
3085: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3087: PetscFunctionBegin;
3088: data->correction = type;
3089: PetscFunctionReturn(PETSC_SUCCESS);
3090: }
3092: static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type)
3093: {
3094: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3096: PetscFunctionBegin;
3097: *type = data->correction;
3098: PetscFunctionReturn(PETSC_SUCCESS);
3099: }
3101: /*@
3102: PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared.
3104: Input Parameters:
3105: + pc - preconditioner context
3106: - share - whether the `KSP` should be shared or not
3108: Note:
3109: This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped
3110: when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`.
3112: Level: advanced
3114: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()`
3115: @*/
3116: PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share)
3117: {
3118: PetscFunctionBegin;
3120: PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share));
3121: PetscFunctionReturn(PETSC_SUCCESS);
3122: }
3124: /*@
3125: PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared.
3127: Input Parameter:
3128: . pc - preconditioner context
3130: Output Parameter:
3131: . share - whether the `KSP` is shared or not
3133: Note:
3134: 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
3135: when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`.
3137: Level: advanced
3139: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()`
3140: @*/
3141: PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share)
3142: {
3143: PetscFunctionBegin;
3145: if (share) {
3146: PetscAssertPointer(share, 2);
3147: PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share));
3148: }
3149: PetscFunctionReturn(PETSC_SUCCESS);
3150: }
3152: static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share)
3153: {
3154: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3156: PetscFunctionBegin;
3157: data->share = share;
3158: PetscFunctionReturn(PETSC_SUCCESS);
3159: }
3161: static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share)
3162: {
3163: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3165: PetscFunctionBegin;
3166: *share = data->share;
3167: PetscFunctionReturn(PETSC_SUCCESS);
3168: }
3170: /*@
3171: PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator.
3173: Input Parameters:
3174: + pc - preconditioner context
3175: . is - index set of the local deflation matrix
3176: - U - deflation sequential matrix stored as a `MATSEQDENSE`
3178: Level: advanced
3180: .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()`
3181: @*/
3182: PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U)
3183: {
3184: PetscFunctionBegin;
3188: PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U));
3189: PetscFunctionReturn(PETSC_SUCCESS);
3190: }
3192: static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U)
3193: {
3194: PetscFunctionBegin;
3195: PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE));
3196: PetscFunctionReturn(PETSC_SUCCESS);
3197: }
3199: PetscErrorCode HPDDMLoadDL_Private(PetscBool *found)
3200: {
3201: PetscBool flg;
3202: char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN];
3204: PetscFunctionBegin;
3205: PetscAssertPointer(found, 1);
3206: PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir)));
3207: PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr));
3208: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir));
3209: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3210: #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */
3211: if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */
3212: PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir)));
3213: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir));
3214: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3215: }
3216: #endif
3217: if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */
3218: PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg));
3219: if (flg) { /* if both PETSc and SLEPc are configured with --download-hpddm but PETSc has been configured without --download-slepc, one must ensure that libslepc is loaded before libhpddm_petsc */
3220: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir));
3221: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3222: PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir);
3223: PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib));
3224: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */
3225: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3226: }
3227: }
3228: PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib);
3229: PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib));
3230: PetscFunctionReturn(PETSC_SUCCESS);
3231: }
3233: /*MC
3234: PCHPDDM - Interface with the HPDDM library.
3236: This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`.
3237: It may be viewed as an alternative to spectral
3238: AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020`
3240: The matrix used for building the preconditioner (Pmat) may be unassembled (`MATIS`), assembled (`MATAIJ`, `MATBAIJ`, or `MATSBAIJ`), hierarchical (`MATHTOOL`), `MATNORMAL`, `MATNORMALHERMITIAN`, or `MATSCHURCOMPLEMENT` (when `PCHPDDM` is used as part of an outer `PCFIELDSPLIT`).
3242: For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using
3243: `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`.
3245: Options Database Keys:
3246: + -pc_hpddm_define_subdomains (true|false) - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()`
3247: (not relevant with an unassembled Pmat)
3248: . -pc_hpddm_has_neumann (true|false) - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()`
3249: - -pc_hpddm_coarse_correction type - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` default is `deflated`
3251: Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes.
3252: .vb
3253: -pc_hpddm_levels_%d_pc_
3254: -pc_hpddm_levels_%d_ksp_
3255: -pc_hpddm_levels_%d_eps_
3256: -pc_hpddm_levels_%d_p
3257: -pc_hpddm_levels_%d_mat_type
3258: -pc_hpddm_coarse_
3259: -pc_hpddm_coarse_p
3260: -pc_hpddm_coarse_mat_type
3261: -pc_hpddm_coarse_mat_filter
3262: .ve
3264: 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
3265: -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij` will use 10 deflation vectors per subdomain on the fine "level 1",
3266: aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2",
3267: and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`).
3269: 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_absolute val`. The default `-pc_hpddm_coarse_p value` is 1, meaning that the coarse operator is aggregated on a single process.
3271: Level: intermediate
3273: Notes:
3274: This preconditioner requires that PETSc is built with SLEPc (`--download-slepc`).
3276: By default, the underlying concurrent eigenproblems
3277: are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf.
3278: {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As
3279: 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_nev 10
3280: -pc_hpddm_levels_1_st_type sinvert`. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in
3281: SLEPc documentation since they are specific to `PCHPDDM`.
3282: .vb
3283: -pc_hpddm_levels_1_st_share_sub_ksp
3284: -pc_hpddm_levels_%d_eps_threshold_absolute
3285: -pc_hpddm_levels_1_eps_use_inertia
3286: .ve
3288: The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is
3289: used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse
3290: correction are associated to eigenvalues whose magnitude are lower or equal than `-pc_hpddm_levels_N_eps_threshold_absolute`. When using an `EPS` which cannot
3291: determine a priori the proper `-pc_hpddm_levels_N_eps_nev` such that all wanted eigenmodes are retrieved, it is possible to get an estimation of the
3292: correct value using the third option from the list, `-pc_hpddm_levels_1_eps_use_inertia`, see `MatGetInertia()`. In that case, there is no need
3293: to supply `-pc_hpddm_levels_1_eps_nev`. This last option also only applies to the fine-level (N = 1) eigensolver.
3295: See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent`
3297: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`,
3298: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`,
3299: `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()`
3300: M*/
3301: PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc)
3302: {
3303: PC_HPDDM *data;
3304: PetscBool found;
3306: PetscFunctionBegin;
3307: if (!loadedSym) {
3308: PetscCall(HPDDMLoadDL_Private(&found));
3309: if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym));
3310: }
3311: PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc");
3312: PetscCall(PetscNew(&data));
3313: pc->data = data;
3314: data->Neumann = PETSC_BOOL3_UNKNOWN;
3315: pc->ops->reset = PCReset_HPDDM;
3316: pc->ops->destroy = PCDestroy_HPDDM;
3317: pc->ops->setfromoptions = PCSetFromOptions_HPDDM;
3318: pc->ops->setup = PCSetUp_HPDDM;
3319: pc->ops->apply = PCApply_HPDDM<false>;
3320: pc->ops->matapply = PCMatApply_HPDDM<false>;
3321: pc->ops->applytranspose = PCApply_HPDDM<true>;
3322: pc->ops->matapplytranspose = PCMatApply_HPDDM<true>;
3323: pc->ops->view = PCView_HPDDM;
3324: pc->ops->presolve = PCPreSolve_HPDDM;
3326: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM));
3327: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM));
3328: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM));
3329: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM));
3330: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM));
3331: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM));
3332: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM));
3333: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM));
3334: PetscFunctionReturn(PETSC_SUCCESS);
3335: }
3337: /*@
3338: PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`.
3340: Level: developer
3342: .seealso: [](ch_ksp), `PetscInitialize()`
3343: @*/
3344: PetscErrorCode PCHPDDMInitializePackage(void)
3345: {
3346: char ename[32];
3348: PetscFunctionBegin;
3349: if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
3350: PCHPDDMPackageInitialized = PETSC_TRUE;
3351: PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage));
3352: /* general events registered once during package initialization */
3353: /* some of these events are not triggered in libpetsc, */
3354: /* but rather directly in libhpddm_petsc, */
3355: /* which is in charge of performing the following operations */
3357: /* domain decomposition structure from Pmat sparsity pattern */
3358: PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc));
3359: /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */
3360: PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP));
3361: /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */
3362: PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP));
3363: /* next level construction using PtAP and PtBP (not triggered in libpetsc) */
3364: PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next));
3365: static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high");
3366: for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) {
3367: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i));
3368: /* events during a PCSetUp() at level #i _except_ the assembly */
3369: /* of the Galerkin operator of the coarser level #(i + 1) */
3370: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1]));
3371: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i));
3372: /* events during a PCApply() at level #i _except_ */
3373: /* the KSPSolve() of the coarser level #(i + 1) */
3374: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1]));
3375: }
3376: PetscFunctionReturn(PETSC_SUCCESS);
3377: }
3379: /*@
3380: PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`.
3382: Level: developer
3384: .seealso: [](ch_ksp), `PetscFinalize()`
3385: @*/
3386: PetscErrorCode PCHPDDMFinalizePackage(void)
3387: {
3388: PetscFunctionBegin;
3389: PCHPDDMPackageInitialized = PETSC_FALSE;
3390: PetscFunctionReturn(PETSC_SUCCESS);
3391: }
3393: static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y)
3394: {
3395: Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */
3396: /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */
3397: Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */
3398: /* [ A_10 A_11 ] R_1^T A_12 x [ ] */
3399: PetscFunctionBegin;
3400: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3401: PetscCall(VecSet(h->v, 0.0));
3402: PetscCall(VecGetSubVector(h->v, h->is[0], &sub));
3403: PetscCall(MatMult(h->A[0], x, sub));
3404: PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub));
3405: PetscCall(KSPSolve(h->ksp, h->v, h->v));
3406: PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y));
3407: PetscFunctionReturn(PETSC_SUCCESS);
3408: }
3410: static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x)
3411: {
3412: Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */
3413: Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */
3415: PetscFunctionBegin;
3416: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3417: PetscCall(VecSet(h->v, 0.0));
3418: PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y));
3419: PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v));
3420: PetscCall(VecGetSubVector(h->v, h->is[0], &sub));
3421: PetscCall(MatMultTranspose(h->A[0], sub, x));
3422: PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub));
3423: PetscFunctionReturn(PETSC_SUCCESS);
3424: }
3426: static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *)
3427: {
3428: Harmonic h;
3429: Mat A, B;
3430: Vec a, b;
3432: PetscFunctionBegin;
3433: PetscCall(MatShellGetContext(S, static_cast<void *>(&h)));
3434: PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A));
3435: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B));
3436: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3437: PetscCall(MatDenseGetColumnVecRead(A, i, &a));
3438: PetscCall(MatDenseGetColumnVecWrite(B, i, &b));
3439: PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a));
3440: PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b));
3441: PetscCall(MatDenseRestoreColumnVecRead(A, i, &a));
3442: }
3443: PetscCall(MatDestroy(&A));
3444: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A));
3445: PetscCall(KSPMatSolve(h->ksp, B, A));
3446: PetscCall(MatDestroy(&B));
3447: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3448: PetscCall(MatDenseGetColumnVecRead(A, i, &a));
3449: PetscCall(MatDenseGetColumnVecWrite(Y, i, &b));
3450: PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b));
3451: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b));
3452: PetscCall(MatDenseRestoreColumnVecRead(A, i, &a));
3453: }
3454: PetscCall(MatDestroy(&A));
3455: PetscFunctionReturn(PETSC_SUCCESS);
3456: }
3458: static PetscErrorCode MatProduct_AtB_Harmonic(Mat S, Mat Y, Mat X, void *)
3459: {
3460: Harmonic h;
3461: Mat A, B;
3462: Vec a, b;
3464: PetscFunctionBegin;
3465: PetscCall(MatShellGetContext(S, static_cast<void *>(&h)));
3466: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, Y->cmap->n, nullptr, &A));
3467: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3468: PetscCall(MatDenseGetColumnVecRead(Y, i, &b));
3469: PetscCall(MatDenseGetColumnVecWrite(A, i, &a));
3470: PetscCall(VecISCopy(a, h->is[1], SCATTER_FORWARD, b));
3471: PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a));
3472: PetscCall(MatDenseRestoreColumnVecRead(Y, i, &b));
3473: }
3474: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B));
3475: PetscCall(KSPMatSolveTranspose(h->ksp, A, B));
3476: PetscCall(MatDestroy(&A));
3477: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->A[0]->rmap->n, B->cmap->n, nullptr, &A));
3478: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3479: PetscCall(MatDenseGetColumnVecRead(B, i, &b));
3480: PetscCall(MatDenseGetColumnVecWrite(A, i, &a));
3481: PetscCall(VecISCopy(b, h->is[0], SCATTER_REVERSE, a));
3482: PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a));
3483: PetscCall(MatDenseRestoreColumnVecRead(B, i, &b));
3484: }
3485: PetscCall(MatDestroy(&B));
3486: PetscCall(MatTransposeMatMult(h->A[0], A, MAT_REUSE_MATRIX, PETSC_CURRENT, &X));
3487: PetscCall(MatDestroy(&A));
3488: PetscFunctionReturn(PETSC_SUCCESS);
3489: }
3491: static PetscErrorCode MatDestroy_Harmonic(Mat A)
3492: {
3493: Harmonic h;
3495: PetscFunctionBegin;
3496: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3497: for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i));
3498: PetscCall(PetscFree(h->is));
3499: PetscCall(VecDestroy(&h->v));
3500: for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i));
3501: PetscCall(PetscFree(h->A));
3502: PetscCall(KSPDestroy(&h->ksp));
3503: PetscCall(PetscFree(h));
3504: PetscFunctionReturn(PETSC_SUCCESS);
3505: }