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