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, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = 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, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = 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, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr>
1113: static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type x, Type y)
1114: {
1115: PetscFunctionBegin;
1116: PetscCall(VecISCopy(std::get<2>(*p)[0], std::get<1>(*p), SCATTER_FORWARD, x));
1117: if (!T) PetscCall(PCApply(factor, std::get<2>(*p)[0], std::get<2>(*p)[1]));
1118: else PetscCall(PCApplyTranspose(factor, std::get<2>(*p)[0], std::get<2>(*p)[1]));
1119: PetscCall(VecISCopy(std::get<2>(*p)[1], std::get<1>(*p), SCATTER_REVERSE, y));
1120: PetscFunctionReturn(PETSC_SUCCESS);
1121: }
1123: template <class Type, bool = false, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr>
1124: static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type X, Type Y)
1125: {
1126: Mat B[2];
1127: Vec x, y;
1129: PetscFunctionBegin;
1130: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B));
1131: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B + 1));
1132: for (PetscInt i = 0; i < X->cmap->n; ++i) {
1133: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
1134: PetscCall(MatDenseGetColumnVecWrite(B[0], i, &y));
1135: PetscCall(VecISCopy(y, std::get<1>(*p), SCATTER_FORWARD, x));
1136: PetscCall(MatDenseRestoreColumnVecWrite(B[0], i, &y));
1137: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
1138: }
1139: PetscCall(PCMatApply(factor, B[0], B[1]));
1140: PetscCall(MatDestroy(B));
1141: for (PetscInt i = 0; i < X->cmap->n; ++i) {
1142: PetscCall(MatDenseGetColumnVecRead(B[1], i, &x));
1143: PetscCall(MatDenseGetColumnVecWrite(Y, i, &y));
1144: PetscCall(VecISCopy(x, std::get<1>(*p), SCATTER_REVERSE, y));
1145: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y));
1146: PetscCall(MatDenseRestoreColumnVecRead(B[1], i, &x));
1147: }
1148: PetscCall(MatDestroy(B + 1));
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: template <class Type = Vec, bool T = false>
1153: static PetscErrorCode PCApply_Schur(PC pc, Type x, Type y)
1154: {
1155: PC factor;
1156: Mat A;
1157: MatSolverType type;
1158: PetscBool flg;
1159: std::tuple<KSP, IS, Vec[2]> *p;
1161: PetscFunctionBegin;
1162: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1163: PetscCall(KSPGetPC(std::get<0>(*p), &factor));
1164: PetscCall(PCFactorGetMatSolverType(factor, &type));
1165: PetscCall(PCFactorGetMatrix(factor, &A));
1166: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
1167: if (flg) {
1168: PetscCheck(PetscDefined(HAVE_MUMPS), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType");
1169: PetscCall(MatMumpsSetIcntl(A, 26, 0));
1170: } else {
1171: PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &flg));
1172: PetscCheck(flg && PetscDefined(HAVE_MKL_PARDISO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType");
1173: flg = PETSC_FALSE;
1174: #if PetscDefined(HAVE_MKL_PARDISO)
1175: PetscCall(MatMkl_PardisoSetCntl(A, 70, 1));
1176: #endif
1177: }
1178: PetscCall(PCApply_Schur_Private<Type, T>(p, factor, x, y));
1179: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1));
1180: else {
1181: #if PetscDefined(HAVE_MKL_PARDISO)
1182: PetscCall(MatMkl_PardisoSetCntl(A, 70, 0));
1183: #endif
1184: }
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: static PetscErrorCode PCDestroy_Schur(PC pc)
1189: {
1190: std::tuple<KSP, IS, Vec[2]> *p;
1192: PetscFunctionBegin;
1193: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1194: PetscCall(ISDestroy(&std::get<1>(*p)));
1195: PetscCall(VecDestroy(std::get<2>(*p)));
1196: PetscCall(VecDestroy(std::get<2>(*p) + 1));
1197: PetscCall(PetscFree(p));
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: template <char trans>
1202: static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu)
1203: {
1204: PC pc;
1205: Mat B, X;
1206: PetscInt n, N, j = 0;
1208: PetscFunctionBegin;
1209: PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr));
1210: PetscCall(MatGetLocalSize(B, &n, nullptr));
1211: PetscCall(MatGetSize(B, &N, nullptr));
1212: PetscCall(PetscIntCast(std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx)), &j));
1213: PetscCheck(j > 0 && j < ctx->parent->N, PetscObjectComm((PetscObject)ctx->ksp), PETSC_ERR_PLIB, "No finer level found for PCHPDDM solve");
1214: PetscCall(KSPGetPC(ctx->parent->levels[j - 1]->ksp, &pc));
1215: if (ctx->parent->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr));
1216: if (mu == 1) {
1217: if (!ctx->ksp->vec_rhs) {
1218: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs));
1219: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol));
1220: }
1221: PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs));
1222: if (trans == 'N') PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr));
1223: else {
1224: if (trans == 'H') PetscCall(VecConjugate(ctx->ksp->vec_rhs));
1225: PetscCall(KSPSolveTranspose(ctx->ksp, nullptr, nullptr)); /* TODO: missing KSPSolveHermitianTranspose() */
1226: if (trans == 'H') PetscCall(VecConjugate(ctx->ksp->vec_sol));
1227: }
1228: PetscCall(KSPCheckSolve(ctx->ksp, pc, ctx->ksp->vec_sol));
1229: PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs));
1230: PetscCall(VecResetArray(ctx->ksp->vec_rhs));
1231: } else {
1232: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B));
1233: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X));
1234: if (trans == 'N') PetscCall(KSPMatSolve(ctx->ksp, B, X));
1235: else {
1236: if (trans == 'H') PetscCall(MatConjugate(B));
1237: PetscCall(KSPMatSolveTranspose(ctx->ksp, B, X)); /* TODO: missing KSPMatSolveHermitianTranspose() */
1238: if (trans == 'H') PetscCall(MatConjugate(X));
1239: }
1240: PetscCall(KSPCheckMatSolve(ctx->ksp, pc, X));
1241: PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN));
1242: PetscCall(MatDestroy(&X));
1243: PetscCall(MatDestroy(&B));
1244: }
1245: if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr));
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc)
1250: {
1251: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1253: PetscFunctionBegin;
1254: if (data->setup) {
1255: Mat P;
1256: Vec x, xt = nullptr;
1257: PetscReal t = 0.0, s = 0.0;
1259: PetscCall(PCGetOperators(pc, nullptr, &P));
1260: PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x));
1261: PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx));
1262: }
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[])
1267: {
1268: Mat A;
1269: PetscBool flg;
1271: PetscFunctionBegin;
1272: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n);
1273: /* previously composed Mat */
1274: PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A));
1275: PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat");
1276: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSCHURCOMPLEMENT, &flg)); /* MATSCHURCOMPLEMENT has neither a MatDuplicate() nor a MatCopy() implementation */
1277: if (scall == MAT_INITIAL_MATRIX) {
1278: PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */
1279: if (!flg) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat));
1280: } else if (!flg) PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN));
1281: if (flg) {
1282: PetscCall(MatDestroy(*submat)); /* previously created Mat has to be destroyed */
1283: (*submat)[0] = A;
1284: PetscCall(PetscObjectReference((PetscObject)A));
1285: }
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted)
1290: {
1291: PetscErrorCodeFn *op;
1293: PetscFunctionBegin;
1294: /* previously-composed Mat */
1295: PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C));
1296: PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op));
1297: /* see https://mailman.cels.anl.gov/archives/list/petsc-dev@lists.mcs.anl.gov/message/22HXNMER6OU7N7CXWA2LJAY6RZPGQYXT/ */
1298: PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (PetscErrorCodeFn *)PCHPDDMCreateSubMatrices_Private));
1299: if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */
1300: PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */
1301: PetscCall(PCSetUp(pc));
1302: /* reset MatCreateSubMatrices() */
1303: PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op));
1304: PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p)
1309: {
1310: IS perm;
1311: const PetscInt *ptr;
1312: PetscInt *compressed, size, bs;
1313: std::map<PetscInt, PetscInt> order;
1314: PetscBool flg;
1316: PetscFunctionBegin;
1319: PetscCall(ISGetLocalSize(is, &size));
1320: PetscCall(ISGetBlockSize(is, &bs));
1321: PetscCall(ISSorted(is, &flg));
1322: if (!flg) {
1323: PetscCall(ISGetIndices(is, &ptr));
1324: /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */
1325: for (PetscInt n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs));
1326: PetscCall(ISRestoreIndices(is, &ptr));
1327: size /= bs;
1328: if (out_C) {
1329: PetscCall(PetscMalloc1(size, &compressed));
1330: for (const std::pair<const PetscInt, PetscInt> &i : order) *compressed++ = i.second;
1331: compressed -= size;
1332: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, compressed, PETSC_OWN_POINTER, &perm));
1333: PetscCall(ISSetPermutation(perm));
1334: /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */
1335: PetscCall(MatPermute(in_C, perm, perm, out_C));
1336: if (p) *p = perm;
1337: else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */
1338: }
1339: if (out_is) {
1340: PetscCall(PetscMalloc1(size, &compressed));
1341: for (const std::pair<const PetscInt, PetscInt> &i : order) *compressed++ = i.first;
1342: compressed -= size;
1343: /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */
1344: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, compressed, PETSC_OWN_POINTER, out_is));
1345: }
1346: } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */
1347: if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C));
1348: if (out_is) {
1349: PetscCall(PetscObjectTypeCompare((PetscObject)in_is, ISBLOCK, &flg));
1350: if (flg) PetscCall(ISDuplicate(in_is, out_is));
1351: else {
1352: PetscCall(ISGetIndices(is, &ptr));
1353: if (bs > 1) {
1354: size /= bs;
1355: PetscCall(PetscMalloc1(size, &compressed));
1356: for (PetscInt n = 0; n < size; ++n) compressed[n] = ptr[n * bs] / bs;
1357: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, compressed, PETSC_OWN_POINTER, out_is));
1358: } else PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), 1, size, ptr, PETSC_COPY_VALUES, out_is));
1359: PetscCall(ISRestoreIndices(is, &ptr));
1360: PetscCall(ISSetInfo(*out_is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1361: }
1362: }
1363: }
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10, Mat *B01 = nullptr)
1368: {
1369: Mat T, U = nullptr, B = nullptr;
1370: IS z;
1371: PetscBool flg, conjugate = PETSC_FALSE;
1373: PetscFunctionBegin;
1374: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
1375: if (B01) *B01 = nullptr;
1376: if (flg) {
1377: 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));
1378: PetscCall(MatTransposeGetMat(A10, &U));
1379: } else {
1380: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1381: if (flg) {
1382: 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));
1383: PetscCall(MatHermitianTransposeGetMat(A10, &U));
1384: conjugate = PETSC_TRUE;
1385: }
1386: }
1387: if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T));
1388: else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T));
1389: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg));
1390: if (flg) {
1391: 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));
1392: PetscCall(MatTransposeGetMat(A01, &A01));
1393: PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B));
1394: A01 = B;
1395: } else {
1396: PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1397: if (flg) {
1398: 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));
1399: PetscCall(MatHermitianTransposeGetMat(A01, &A01));
1400: PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B));
1401: A01 = B;
1402: }
1403: }
1404: PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, &flg));
1405: if (flg) {
1406: PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, &flg));
1407: if (flg) {
1408: PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */
1409: if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */
1410: if (B01) PetscCall(MatDuplicate(T, MAT_COPY_VALUES, B01));
1411: PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE));
1412: PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */
1413: }
1414: PetscCall(MatMultEqual(A01, T, 20, &flg));
1415: if (!B01) PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T");
1416: else {
1417: PetscCall(PetscInfo(pc, "A01 and A10^T are equal? %s\n", PetscBools[flg]));
1418: if (!flg) {
1419: if (z) PetscCall(MatDestroy(&T));
1420: else *B01 = T;
1421: flg = PETSC_TRUE;
1422: } else PetscCall(MatDestroy(B01));
1423: }
1424: PetscCall(ISDestroy(&z));
1425: }
1426: }
1427: if (!flg) PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent layouts, cannot test for equality\n"));
1428: if (!B01 || !*B01) PetscCall(MatDestroy(&T));
1429: else if (conjugate) PetscCall(MatConjugate(T));
1430: PetscCall(MatDestroy(&B));
1431: PetscFunctionReturn(PETSC_SUCCESS);
1432: }
1434: static PetscErrorCode PCHPDDMCheckInclusion_Private(PC pc, IS is, IS is_local, PetscBool check)
1435: {
1436: IS intersect;
1437: const char *str = "IS of the auxiliary Mat does not include all local rows of A";
1438: PetscBool equal;
1440: PetscFunctionBegin;
1441: PetscCall(ISIntersect(is, is_local, &intersect));
1442: PetscCall(ISEqualUnsorted(is_local, intersect, &equal));
1443: PetscCall(ISDestroy(&intersect));
1444: if (check) PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", str);
1445: else if (!equal) PetscCall(PetscInfo(pc, "%s\n", str));
1446: PetscFunctionReturn(PETSC_SUCCESS);
1447: }
1449: static PetscErrorCode PCHPDDMCheckMatStructure_Private(PC pc, Mat A, Mat B)
1450: {
1451: Mat X, Y;
1452: const PetscInt *i[2], *j[2];
1453: PetscBool flg = PETSC_TRUE;
1455: PetscFunctionBegin;
1456: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &X)); /* no common way to compare sparsity pattern, so just convert to MATSEQAIJ */
1457: 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) */
1458: PetscCall(MatSeqAIJGetCSRAndMemType(X, &i[0], &j[0], nullptr, nullptr));
1459: PetscCall(MatSeqAIJGetCSRAndMemType(Y, &i[1], &j[1], nullptr, nullptr));
1460: for (PetscInt row = 0; (row < X->rmap->n) && flg; ++row) {
1461: const PetscInt n = i[0][row + 1] - i[0][row];
1463: for (PetscInt k = i[1][row], location; k < i[1][row + 1]; ++k) {
1464: PetscCall(PetscFindInt(j[1][k], n, j[0] + i[0][row], &location));
1465: if (location < 0) {
1466: flg = PETSC_FALSE;
1467: break;
1468: }
1469: }
1470: }
1471: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)pc)));
1472: 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");
1473: PetscCall(MatDestroy(&Y));
1474: PetscCall(MatDestroy(&X));
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub)
1479: {
1480: IS is;
1482: PetscFunctionBegin;
1483: if (!flg) {
1484: if (algebraic) {
1485: PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is));
1486: PetscCall(ISDestroy(&is));
1487: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr));
1488: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr));
1489: }
1490: PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub));
1491: }
1492: PetscFunctionReturn(PETSC_SUCCESS);
1493: }
1495: static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat Q, IS *is, Mat *sub[], PetscBool block)
1496: {
1497: IS icol[3], irow[2];
1498: Mat *M;
1499: Mat P = Q;
1500: PetscReal *ptr;
1501: PetscInt *idx, p = 0, bs = P->cmap->bs;
1502: PetscBool flg;
1504: PetscFunctionBegin;
1505: /* MatCreateSubMatrices_MPISBAIJ() may return a rectangular MATSEQSBAIJ containing only the explicitly stored upper-triangular entries */
1506: /* of the selected rows. But the missing lower-triangular entries are needed by MatGetColumnNorms(), MatGetRowSum(), and MatMatMult() */
1507: PetscCall(PetscObjectTypeCompare((PetscObject)Q, MATMPISBAIJ, &flg));
1508: if (flg) PetscCall(MatConvert(Q, MATBAIJ, MAT_INITIAL_MATRIX, &P));
1509: PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2));
1510: PetscCall(ISSetBlockSize(icol[2], bs));
1511: PetscCall(ISSetIdentity(icol[2]));
1512: PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M));
1513: if (flg) {
1514: PetscCall(MatDestroy(&P));
1515: P = Q; /* continue using the caller-owned (MATMPISBAIJ) matrix */
1516: }
1517: PetscCall(ISDestroy(icol + 2));
1518: PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow));
1519: PetscCall(ISSetBlockSize(irow[0], bs));
1520: PetscCall(ISSetIdentity(irow[0]));
1521: if (!block) {
1522: PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx));
1523: PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr));
1524: /* check for nonzero columns so that M[0] may be expressed in compact form */
1525: for (PetscInt n = 0; n < P->cmap->N; n += bs) {
1526: if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs;
1527: }
1528: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1));
1529: PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
1530: PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2));
1531: irow[1] = irow[0];
1532: /* 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 */
1533: icol[0] = is[0];
1534: PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub));
1535: PetscCall(ISDestroy(icol + 1));
1536: PetscCall(PetscFree2(ptr, idx));
1537: PetscCall(MatPropagateSymmetryOptions(P, (*sub)[0]));
1538: if (flg) PetscCall(MatConvert((*sub)[0], MATSBAIJ, MAT_INPLACE_MATRIX, sub[0]));
1539: /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */
1540: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2]));
1541: /* Mat used in eq. (3.1) of [2022b] */
1542: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1]));
1543: } else {
1544: Mat aux;
1546: PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
1547: /* diagonal block of the overlapping rows */
1548: PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub));
1549: PetscCall(MatPropagateSymmetryOptions(P, (*sub)[0]));
1550: if (flg && bs == 1) PetscCall(MatConvert((*sub)[0], MATSBAIJ, MAT_INPLACE_MATRIX, sub[0]));
1551: PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux));
1552: aux->spd = PETSC_BOOL3_UNKNOWN; /* the auxiliary Mat need not be SPD */
1553: PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1554: if (bs == 1) { /* scalar case */
1555: Vec sum[2];
1557: PetscCall(MatCreateVecs(aux, sum, sum + 1));
1558: PetscCall(MatGetRowSum(M[0], sum[0]));
1559: PetscCall(MatGetRowSum(aux, sum[1]));
1560: /* off-diagonal block row sum (full rows - diagonal block rows) */
1561: PetscCall(VecAXPY(sum[0], -1.0, sum[1]));
1562: /* subdomain matrix plus off-diagonal block row sum */
1563: PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES));
1564: PetscCall(VecDestroy(sum));
1565: PetscCall(VecDestroy(sum + 1));
1566: } else { /* vectorial case */
1567: /* TODO: missing MatGetValuesBlocked(), so the code below is */
1568: /* an extension of the scalar case for when bs > 1, but it could */
1569: /* be more efficient by avoiding all these MatMatMult() */
1570: Mat sum[2], ones;
1571: PetscScalar *ptr;
1573: aux->symmetry_eternal = PETSC_FALSE;
1574: aux->symmetric = PETSC_BOOL3_UNKNOWN;
1575: aux->hermitian = PETSC_BOOL3_UNKNOWN;
1576: PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr));
1577: PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones));
1578: for (PetscInt n = 0; n < M[0]->cmap->n; n += bs) {
1579: for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0;
1580: }
1581: PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum));
1582: PetscCall(MatDestroy(&ones));
1583: PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones));
1584: PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n));
1585: PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum + 1));
1586: PetscCall(MatDestroy(&ones));
1587: PetscCall(PetscFree(ptr));
1588: /* off-diagonal block row sum (full rows - diagonal block rows) */
1589: PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN));
1590: PetscCall(MatDestroy(sum + 1));
1591: /* re-order values to be consistent with MatSetValuesBlocked() */
1592: /* equivalent to MatTranspose() which does not truly handle */
1593: /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */
1594: PetscCall(MatDenseGetArrayWrite(sum[0], &ptr));
1595: HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs);
1596: /* subdomain matrix plus off-diagonal block row sum */
1597: for (PetscInt n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES));
1598: PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY));
1599: PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY));
1600: PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr));
1601: PetscCall(MatDestroy(sum));
1602: }
1603: PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
1604: /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */
1605: PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux));
1606: }
1607: PetscCall(ISDestroy(irow));
1608: PetscCall(MatDestroySubMatrices(1, &M));
1609: PetscFunctionReturn(PETSC_SUCCESS);
1610: }
1612: static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y)
1613: {
1614: Mat A;
1615: MatSolverType type;
1616: IS is[2];
1617: PetscBool flg;
1618: std::pair<PC, Vec[2]> *p;
1620: PetscFunctionBegin;
1621: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1622: if (p->second[0]) { /* in case of a centralized Schur complement, some processes may have no local operator */
1623: PetscCall(PCGetOperators(p->first, &A, nullptr));
1624: PetscCall(MatNestGetISs(A, is, nullptr));
1625: PetscCall(PetscObjectTypeCompareAny((PetscObject)p->first, &flg, PCLU, PCCHOLESKY, ""));
1626: if (flg) { /* partial solve currently only makes sense with exact factorizations */
1627: PetscCall(PCFactorGetMatSolverType(p->first, &type));
1628: PetscCall(PCFactorGetMatrix(p->first, &A));
1629: if (A->schur) {
1630: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
1631: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */
1632: } else flg = PETSC_FALSE;
1633: }
1634: PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */
1635: PetscCall(PCApply(p->first, p->second[0], p->second[1]));
1636: PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */
1637: if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */
1638: }
1639: PetscFunctionReturn(PETSC_SUCCESS);
1640: }
1642: static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer)
1643: {
1644: std::pair<PC, Vec[2]> *p;
1646: PetscFunctionBegin;
1647: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1648: PetscCall(PCView(p->first, viewer));
1649: PetscFunctionReturn(PETSC_SUCCESS);
1650: }
1652: static PetscErrorCode PCDestroy_Nest(PC pc)
1653: {
1654: std::pair<PC, Vec[2]> *p;
1656: PetscFunctionBegin;
1657: PetscCall(PCShellGetContext(pc, static_cast<void *>(&p)));
1658: PetscCall(VecDestroy(p->second));
1659: PetscCall(VecDestroy(p->second + 1));
1660: PetscCall(PCDestroy(&p->first));
1661: PetscCall(PetscFree(p));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: template <bool T = false>
1666: static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y)
1667: {
1668: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
1670: PetscFunctionBegin;
1671: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1672: PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */
1673: PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD));
1674: if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */
1675: else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1]));
1676: PetscCall(VecSet(y, 0.0));
1677: 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 */
1678: PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE));
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: static PetscErrorCode MatDestroy_Schur(Mat A)
1683: {
1684: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
1686: PetscFunctionBegin;
1687: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1688: PetscCall(VecDestroy(std::get<2>(*ctx)));
1689: PetscCall(VecDestroy(std::get<2>(*ctx) + 1));
1690: PetscCall(PetscFree(ctx));
1691: PetscFunctionReturn(PETSC_SUCCESS);
1692: }
1694: static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y)
1695: {
1696: PC pc;
1697: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1699: PetscFunctionBegin;
1700: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1701: pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc;
1702: 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 */
1703: PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */
1704: PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */
1705: PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */
1706: PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */
1707: } else {
1708: PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */
1709: PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */
1710: PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */
1711: PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */
1712: }
1713: PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */
1714: PetscFunctionReturn(PETSC_SUCCESS);
1715: }
1717: static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer)
1718: {
1719: PetscBool ascii;
1720: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1722: PetscFunctionBegin;
1723: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
1724: if (ascii) {
1725: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1726: 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)"));
1727: PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */
1728: }
1729: PetscFunctionReturn(PETSC_SUCCESS);
1730: }
1732: static PetscErrorCode MatDestroy_SchurCorrection(Mat A)
1733: {
1734: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx;
1736: PetscFunctionBegin;
1737: PetscCall(MatShellGetContext(A, static_cast<void *>(&ctx)));
1738: PetscCall(VecDestroy(std::get<3>(*ctx)));
1739: PetscCall(VecDestroy(std::get<3>(*ctx) + 1));
1740: PetscCall(VecDestroy(std::get<3>(*ctx) + 2));
1741: PetscCall(PCDestroy(std::get<0>(*ctx) + 1));
1742: PetscCall(PetscFree(ctx));
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: static PetscErrorCode PCPostSolve_SchurPreLeastSquares(PC, KSP, Vec, Vec x)
1747: {
1748: PetscFunctionBegin;
1749: PetscCall(VecScale(x, -1.0));
1750: PetscFunctionReturn(PETSC_SUCCESS);
1751: }
1753: static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context)
1754: {
1755: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context);
1757: PetscFunctionBegin;
1758: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) {
1759: PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2]));
1760: 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 */
1761: }
1762: PetscFunctionReturn(PETSC_SUCCESS);
1763: }
1765: static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context)
1766: {
1767: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context);
1769: PetscFunctionBegin;
1770: 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 */
1771: else {
1772: PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2]));
1773: PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */
1774: }
1775: PetscFunctionReturn(PETSC_SUCCESS);
1776: }
1778: static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec);
1779: static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec);
1780: static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *);
1781: static PetscErrorCode MatProduct_AtB_Harmonic(Mat, Mat, Mat, void *);
1782: static PetscErrorCode MatDestroy_Harmonic(Mat);
1784: static PetscErrorCode PCHPDDMComputeScaling_Private(PC_HPDDM *data, IS is, Vec scaling)
1785: {
1786: IS perm;
1787: PetscInt n, m;
1789: PetscFunctionBegin;
1790: PetscCall(ISGetLocalSize(is, &n));
1791: PetscCall(VecGetSize(data->scaling, &m));
1792: PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Overlapping subdomain size differs from the supplied scaling vector size");
1793: PetscCall(ISEmbed(is, data->is, PETSC_TRUE, &perm));
1794: PetscCall(ISGetLocalSize(perm, &m));
1795: 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);
1796: PetscCall(VecISCopy(data->scaling, perm, SCATTER_REVERSE, scaling));
1797: PetscCall(ISDestroy(&perm));
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: }
1801: static PetscErrorCode PCSetUp_HPDDM(PC pc)
1802: {
1803: PC_HPDDM *data = (PC_HPDDM *)pc->data;
1804: PC inner;
1805: KSP *ksp = nullptr;
1806: Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S, ev, daux = nullptr;
1807: Vec xin, v;
1808: std::vector<Vec> initial;
1809: IS is[1], loc, uis = data->is, unsorted = nullptr, dis = nullptr;
1810: ISLocalToGlobalMapping l2g;
1811: char prefix[256];
1812: const char *pcpre;
1813: PetscInt n, requested, reused = 0, overlap = -1;
1814: MatStructure structure = UNKNOWN_NONZERO_PATTERN;
1815: PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE;
1816: DM dm;
1817: std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr;
1819: PetscFunctionBegin;
1820: if (!data->levels) PetscCall(PetscInfo(pc, "No level allocated, defaulting to a single level, PCSetFromOptions() should be called before PCSetUp() to avoid this\n"));
1821: PetscCall(PCHPDDMInitializeLevels_Private(data));
1822: requested = data->N;
1823: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
1824: PetscCall(PCGetOperators(pc, &A, &P));
1825: if (data->scaling) {
1826: PetscCheck(data->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "User-supplied scaling requires an IS, passed to PCHPDDMSetScaling() or PCHPDDMSetDeflationMat()");
1827: if (data->deflation) {
1828: PetscCheck(data->N == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "User-supplied scaling with a deflation space requires exactly two levels");
1829: PetscCall(MatGetSize(data->aux, nullptr, &n));
1830: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &n, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1831: PetscCheck(n > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "User-supplied scaling with a deflation space requires at least one deflation vector globally");
1832: } 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");
1833: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)P, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, ""));
1834: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "User-supplied scaling requires an assembled MATAIJ, MATBAIJ, or MATSBAIJ PMat");
1835: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr));
1836: PetscCheck(subdomains, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "User-supplied scaling requires -%spc_hpddm_define_subdomains true", pcpre ? pcpre : "");
1837: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse"));
1838: n = 0;
1839: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, prefix, "-pc_asm_overlap", &n, nullptr));
1840: 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);
1841: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr));
1842: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_block_splitting", &block, nullptr));
1843: 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 : "");
1844: }
1845: if (!data->levels[0]->ksp) {
1846: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp));
1847: PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel));
1848: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse"));
1849: PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix));
1850: PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY));
1851: } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled && data->levels[0]->ksp->pc->reusepreconditioner) {
1852: /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */
1853: /* then just propagate the appropriate flag to the coarser levels */
1854: for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1855: /* the following KSP and PC may be NULL for some processes, hence the check */
1856: if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE));
1857: if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE));
1858: }
1859: /* early bail out because there is nothing to do */
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: } else {
1862: /* reset coarser levels */
1863: for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
1864: 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) {
1865: reused = data->N - n;
1866: break;
1867: }
1868: PetscCall(KSPDestroy(&data->levels[n]->ksp));
1869: PetscCall(PCDestroy(&data->levels[n]->pc));
1870: }
1871: /* check if some coarser levels are being reused */
1872: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
1873: const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &i_0;
1875: if (!data->deflation && *addr != 0 && reused != data->N - 1) {
1876: /* reuse previously computed eigenvectors */
1877: ev = data->levels[0]->P->getMat();
1878: if (ev) {
1879: initial.reserve(*addr);
1880: for (n = 0; n < *addr; ++n) {
1881: PetscCall(MatDenseGetColumnVecRead(ev, n, &xin));
1882: PetscCall(VecDuplicate(xin, &v));
1883: PetscCall(VecCopy(xin, v));
1884: initial.emplace_back(v);
1885: PetscCall(MatDenseRestoreColumnVecRead(ev, n, &xin));
1886: }
1887: }
1888: }
1889: }
1890: data->N -= reused;
1891: PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P));
1893: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis));
1894: if (!data->is && !ismatis) {
1895: PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr;
1896: PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr;
1897: void *uctx = nullptr;
1899: /* first see if we can get the data from the DM */
1900: PetscCall(MatGetDM(P, &dm));
1901: if (!dm) PetscCall(MatGetDM(A, &dm));
1902: if (!dm) PetscCall(PCGetDM(pc, &dm));
1903: if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */
1904: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create));
1905: if (create) {
1906: PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx));
1907: if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */
1908: }
1909: }
1910: if (!create) {
1911: if (!uis) {
1912: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis));
1913: PetscCall(PetscObjectReference((PetscObject)uis));
1914: }
1915: if (!uaux) {
1916: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux));
1917: PetscCall(PetscObjectReference((PetscObject)uaux));
1918: }
1919: /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */
1920: if (!uis) {
1921: PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis));
1922: PetscCall(PetscObjectReference((PetscObject)uis));
1923: }
1924: if (!uaux) {
1925: PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux));
1926: PetscCall(PetscObjectReference((PetscObject)uaux));
1927: }
1928: }
1929: PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx));
1930: PetscCall(MatDestroy(&uaux));
1931: PetscCall(ISDestroy(&uis));
1932: }
1934: if (!ismatis) {
1935: PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc));
1936: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_block_splitting", &block, nullptr));
1937: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr));
1938: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg));
1939: if (data->is || flg) {
1940: if (block || overlap != -1) {
1941: PetscCall(ISDestroy(&data->is));
1942: PetscCall(MatDestroy(&data->aux));
1943: } else if (flg) {
1944: PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO;
1946: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg));
1947: if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) {
1948: PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */
1949: PetscCall(MatDestroy(&data->aux));
1950: } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) {
1951: PetscContainer container = nullptr;
1953: PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container));
1954: if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */
1955: PC_HPDDM *data_00;
1956: KSP ksp, inner_ksp;
1957: PC pc_00;
1958: Mat A11;
1959: Vec d = nullptr;
1960: PetscReal norm;
1961: const PetscInt *ranges;
1962: PetscMPIInt size;
1963: char *prefix;
1965: PetscCall(MatSchurComplementGetKSP(P, &ksp));
1966: PetscCall(KSPGetPC(ksp, &pc_00));
1967: PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg));
1968: 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 : "",
1969: ((PetscObject)pc_00)->type_name, PCHPDDM);
1970: data_00 = (PC_HPDDM *)pc_00->data;
1971: 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],
1972: data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix);
1973: 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() : "");
1974: PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg));
1975: 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,
1976: ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM);
1977: 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) */
1978: PetscCall(MatSchurComplementGetSubMatrices(P, nullptr, nullptr, nullptr, nullptr, &A11));
1979: PetscCall(MatGetOwnershipRanges(P, &ranges));
1980: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)P), &size));
1981: flg = PetscBool(std::find_if(ranges, ranges + size + 1, [&](PetscInt v) { return v != ranges[0] && v != ranges[size]; }) == ranges + size + 1); /* are all local matrices but one of dimension 0 (centralized Schur complement)? */
1982: if (!flg) {
1983: if (PetscDefined(USE_DEBUG) || !data->is) {
1984: Mat A01, A10, B = nullptr, C = nullptr, *sub;
1986: PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr));
1987: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
1988: if (flg) {
1989: PetscCall(MatTransposeGetMat(A10, &C));
1990: PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B));
1991: } else {
1992: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
1993: if (flg) {
1994: PetscCall(MatHermitianTransposeGetMat(A10, &C));
1995: PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B));
1996: }
1997: }
1998: if (flg)
1999: 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));
2000: if (!B) {
2001: B = A10;
2002: PetscCall(PetscObjectReference((PetscObject)B));
2003: } else if (!data->is) {
2004: PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
2005: if (!flg) C = A01;
2006: else
2007: 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));
2008: }
2009: PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis));
2010: PetscCall(ISSetIdentity(uis));
2011: if (!data->is) {
2012: if (!C) PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C));
2013: else PetscCall(PetscObjectReference((PetscObject)C));
2014: PetscCall(ISDuplicate(data_00->is, is));
2015: PetscCall(MatIncreaseOverlap(A, 1, is, 1));
2016: PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2017: PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub));
2018: PetscCall(MatDestroy(&C));
2019: PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C));
2020: PetscCall(MatDestroySubMatrices(1, &sub));
2021: PetscCall(MatFindNonzeroRows(C, &data->is));
2022: PetscCheck(data->is, PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "No empty row, which likely means that some rows of A_10 are dense");
2023: PetscCall(MatDestroy(&C));
2024: PetscCall(ISDestroy(is));
2025: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc));
2026: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_FALSE));
2027: PetscCall(ISExpand(data->is, loc, is));
2028: PetscCall(ISDestroy(&loc));
2029: PetscCall(ISDestroy(&data->is));
2030: data->is = is[0];
2031: is[0] = nullptr;
2032: }
2033: if (PetscDefined(USE_DEBUG)) {
2034: PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10));
2035: 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 */
2036: PetscCall(ISDestroy(&uis));
2037: PetscCall(ISDuplicate(data->is, &uis));
2038: PetscCall(ISSort(uis));
2039: PetscCall(ISComplement(uis, 0, B->rmap->N, is));
2040: PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C));
2041: PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr));
2042: PetscCall(ISDestroy(is));
2043: PetscCall(MatMultEqual(sub[0], C, 20, &flg));
2044: 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 */
2045: PetscCall(MatDestroy(&C));
2046: PetscCall(MatDestroySubMatrices(1, &sub));
2047: }
2048: PetscCall(ISDestroy(&uis));
2049: PetscCall(MatDestroy(&B));
2050: }
2051: flg = PETSC_FALSE;
2052: if (!data->aux && A11) {
2053: Mat D;
2055: PetscCall(MatCreateVecs(A11, &d, nullptr));
2056: PetscCall(MatGetDiagonal(A11, d));
2057: PetscCall(PetscObjectTypeCompareAny((PetscObject)A11, &flg, MATDIAGONAL, MATCONSTANTDIAGONAL, ""));
2058: if (!flg) {
2059: PetscCall(MatCreateDiagonal(d, &D));
2060: PetscCall(MatMultEqual(A11, D, 20, &flg));
2061: PetscCall(MatDestroy(&D));
2062: }
2063: 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"));
2064: }
2065: if (PetscDefined(USE_DEBUG) || (data->Neumann != PETSC_BOOL3_TRUE && !flg)) {
2066: if (A11) PetscCall(MatNorm(A11, NORM_INFINITY, &norm));
2067: else norm = 0.0;
2068: if (data->Neumann != PETSC_BOOL3_TRUE && !flg) {
2069: 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 : "");
2070: 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"));
2071: PetscCall(MatDestroy(&data->aux));
2072: flg = PETSC_TRUE;
2073: }
2074: }
2075: 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 */
2076: PetscSF scatter;
2077: const PetscScalar *read;
2078: PetscScalar *write, *diagonal = nullptr;
2080: PetscCall(MatDestroy(&data->aux));
2081: PetscCall(ISGetLocalSize(data->is, &n));
2082: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &xin));
2083: PetscCall(VecDuplicate(xin, &v));
2084: PetscCall(VecScatterCreate(xin, data->is, v, nullptr, &scatter));
2085: PetscCall(VecSet(v, 1.0));
2086: PetscCall(VecSet(xin, 1.0));
2087: PetscCall(VecScatterBegin(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE));
2088: PetscCall(VecScatterEnd(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */
2089: PetscCall(PetscSFDestroy(&scatter));
2090: if (d) {
2091: PetscCall(VecScatterCreate(d, data->is, v, nullptr, &scatter));
2092: PetscCall(VecScatterBegin(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD));
2093: PetscCall(VecScatterEnd(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD));
2094: PetscCall(PetscSFDestroy(&scatter));
2095: PetscCall(VecDestroy(&d));
2096: PetscCall(PetscMalloc1(n, &diagonal));
2097: PetscCall(VecGetArrayRead(v, &read));
2098: PetscCallCXX(std::copy_n(read, n, diagonal));
2099: PetscCall(VecRestoreArrayRead(v, &read));
2100: }
2101: PetscCall(VecDestroy(&v));
2102: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v));
2103: PetscCall(VecGetArrayRead(xin, &read));
2104: PetscCall(VecGetArrayWrite(v, &write));
2105: 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];
2106: PetscCall(PetscFree(diagonal));
2107: PetscCall(VecRestoreArrayRead(xin, &read));
2108: PetscCall(VecRestoreArrayWrite(v, &write));
2109: PetscCall(VecDestroy(&xin));
2110: PetscCall(MatCreateDiagonal(v, &data->aux));
2111: PetscCall(VecDestroy(&v));
2112: }
2113: uis = data->is;
2114: uaux = data->aux;
2115: PetscCall(PetscObjectReference((PetscObject)uis));
2116: PetscCall(PetscObjectReference((PetscObject)uaux));
2117: PetscCall(PetscStrallocpy(pcpre, &prefix));
2118: PetscCall(PCSetOptionsPrefix(pc, nullptr));
2119: PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */
2120: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */
2121: PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n));
2122: PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2));
2123: PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat));
2124: PetscCall(KSPSetOptionsPrefix(inner_ksp, prefix));
2125: PetscCall(KSPAppendOptionsPrefix(inner_ksp, "pc_hpddm_"));
2126: PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE));
2127: PetscCall(KSPSetFromOptions(inner_ksp));
2128: PetscCall(KSPGetPC(inner_ksp, &inner));
2129: PetscCall(PCSetOptionsPrefix(inner, nullptr));
2130: PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */
2131: PetscCall(PCKSPSetKSP(pc, inner_ksp));
2132: std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */
2133: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1]));
2134: PetscCall(PCSetOptionsPrefix(pc, prefix)); /* both PC share the same prefix so that the outer PC can be reset with PCSetFromOptions() */
2135: PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix));
2136: PetscCall(PetscFree(prefix));
2137: PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat));
2138: PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM));
2139: 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 */
2140: if (flg) static_cast<PC_HPDDM *>(std::get<0>(*ctx)[1]->data)->Neumann = PETSC_BOOL3_TRUE;
2141: else if (PetscDefined(USE_DEBUG) && norm > PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0)) {
2142: /* no check when A11 is near zero */
2143: PetscCall(MatCreateSubMatrices(A11, 1, &uis, &uis, MAT_INITIAL_MATRIX, &sub));
2144: PetscCall(PCHPDDMCheckMatStructure_Private(pc, sub[0], uaux));
2145: PetscCall(MatDestroySubMatrices(1, &sub));
2146: }
2147: PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1]));
2148: PetscCall(PetscObjectDereference((PetscObject)uis));
2149: PetscCall(PetscObjectDereference((PetscObject)uaux));
2150: 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 */
2151: PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_SchurCorrection));
2152: PetscCall(MatShellSetOperation(S, MATOP_VIEW, (PetscErrorCodeFn *)MatView_SchurCorrection));
2153: PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_SchurCorrection));
2154: PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx))));
2155: if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) {
2156: PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx));
2157: } else { /* no support for PC_SYMMETRIC */
2158: 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]);
2159: }
2160: PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx));
2161: PetscCall(PetscObjectContainerCompose((PetscObject)std::get<0>(*ctx)[1], "_PCHPDDM_Schur", ctx, nullptr));
2162: PetscCall(PCSetUp(std::get<0>(*ctx)[1]));
2163: PetscCall(KSPSetOperators(inner_ksp, S, S));
2164: PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1));
2165: PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2));
2166: PetscCall(PetscObjectDereference((PetscObject)inner_ksp));
2167: PetscCall(PetscObjectDereference((PetscObject)S));
2168: } else {
2169: std::get<0>(*ctx)[0] = pc_00;
2170: PetscCall(PetscObjectContainerCompose((PetscObject)pc, "_PCHPDDM_Schur", ctx, nullptr));
2171: 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 */
2172: if (A11) {
2173: PetscCall(MatGetDiagonalBlock(A11, &data->aux));
2174: PetscCall(PetscObjectReference((PetscObject)data->aux));
2175: } else {
2176: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, P->rmap->n, P->rmap->n, 0, nullptr, &data->aux));
2177: PetscCall(MatAssemblyBegin(data->aux, MAT_FINAL_ASSEMBLY));
2178: PetscCall(MatAssemblyEnd(data->aux, MAT_FINAL_ASSEMBLY));
2179: }
2180: PetscCall(PCSetUp(pc));
2181: }
2182: for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it));
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */
2185: PetscCall(PetscContainerGetPointer(container, static_cast<void *>(&ctx)));
2186: }
2187: }
2188: }
2189: }
2190: if (!data->is && data->N > 1) {
2191: char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */
2193: PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, ""));
2194: if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) {
2195: Mat B;
2197: PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre));
2198: if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED;
2199: PetscCall(MatDestroy(&B));
2200: } else {
2201: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg));
2202: if (flg) {
2203: Mat A00, P00, A01, A10, A11, B, N;
2204: PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES;
2206: PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11));
2207: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg));
2208: if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) {
2209: Mat B01;
2210: Vec diagonal = nullptr;
2211: const PetscScalar *array;
2212: MatSchurComplementAinvType type;
2214: PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10, &B01));
2215: if (A11) {
2216: PetscCall(MatCreateVecs(A11, &diagonal, nullptr));
2217: PetscCall(MatGetDiagonal(A11, diagonal));
2218: }
2219: PetscCall(MatCreateVecs(P00, &v, nullptr));
2220: PetscCall(MatSchurComplementGetAinvType(P, &type));
2221: 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",
2222: ((PetscObject)P)->prefix ? ((PetscObject)P)->prefix : "", MatSchurComplementAinvTypes[type]);
2223: if (type != MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
2224: if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
2225: PetscCall(MatGetRowSum(P00, v));
2226: if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00));
2227: PetscCall(MatDestroy(&P00));
2228: PetscCall(VecGetArrayRead(v, &array));
2229: PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00));
2230: PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2231: for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES));
2232: PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY));
2233: PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY));
2234: PetscCall(VecRestoreArrayRead(v, &array));
2235: PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */
2236: PetscCall(MatDestroy(&P00));
2237: } else PetscCall(MatGetDiagonal(P00, v));
2238: PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */
2239: PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */
2240: PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B));
2241: PetscCall(MatDiagonalScale(B, v, nullptr));
2242: if (B01) PetscCall(MatDiagonalScale(B01, v, nullptr));
2243: } else {
2244: Mat D00;
2245: MatType type;
2247: PetscCall(MatCreate(PetscObjectComm((PetscObject)A00), &D00));
2248: PetscCall(MatSetType(D00, MATAIJ));
2249: PetscCall(MatSetOptionsPrefix(D00, ((PetscObject)A00)->prefix));
2250: PetscCall(MatAppendOptionsPrefix(D00, "block_diagonal_"));
2251: PetscCall(MatSetFromOptions(D00)); /* for setting -mat_block_size dynamically */
2252: 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 */
2253: PetscCall(MatSetBlockSizesFromMats(B, D00, D00));
2254: PetscCall(MatInvertBlockDiagonalMat(B, D00));
2255: PetscCall(MatDestroy(&B));
2256: PetscCall(MatGetType(A01, &type)); /* cache MatType */
2257: PetscCall(MatConvert(A01, MATAIJ, MAT_INPLACE_MATRIX, &A01)); /* MatProduct is not versatile enough to fallback gracefully if no implementation found, so MatConvert() */
2258: PetscCall(MatMatMult(D00, A01, MAT_INITIAL_MATRIX, PETSC_CURRENT, &B));
2259: PetscCall(MatDestroy(&D00));
2260: PetscCall(MatConvert(A01, type, MAT_INPLACE_MATRIX, &A01)); /* reset to previous MatType */
2261: PetscCall(MatConvert(B, type, MAT_INPLACE_MATRIX, &B));
2262: if (!B01) { /* symmetric case */
2263: B01 = A01;
2264: PetscCall(PetscObjectReference((PetscObject)B01));
2265: }
2266: }
2267: if (B01 && B01 != A01) PetscCall(MatSetBlockSizesFromMats(B01, A01, A01)); /* TODO: remove this line once Firedrake is fixed */
2268: PetscCall(VecDestroy(&v));
2269: PetscCall(MatCreateNormalHermitian(B, &N));
2270: PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal, B01));
2271: PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg));
2272: if (!flg) {
2273: PetscCall(MatDestroy(&P));
2274: P = N;
2275: PetscCall(PetscObjectReference((PetscObject)P));
2276: }
2277: if (diagonal) {
2278: PetscCall(MatSetOption(P, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); /* may have missing diagonal entries */
2279: PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES));
2280: PetscCall(PCSetOperators(pc, P, P)); /* replace P by A01^T inv(diag(P00)) A01 - diag(P11) */
2281: PetscCall(VecDestroy(&diagonal));
2282: } else PetscCall(PCSetOperators(pc, B01 ? P : N, P)); /* replace P by A01^T inv(diag(P00)) A01 */
2283: pc->ops->postsolve = PCPostSolve_SchurPreLeastSquares; /* PCFIELDSPLIT expect a KSP for (P11 - A10 inv(diag(P00)) A01) */
2284: PetscCall(MatDestroy(&N)); /* but a PC for (A10 inv(diag(P00)) A10 - P11) is setup instead */
2285: PetscCall(MatDestroy(&P)); /* so the sign of the solution must be flipped */
2286: PetscCall(MatDestroy(&B));
2287: } else
2288: 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 : "");
2289: for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it));
2290: PetscFunctionReturn(PETSC_SUCCESS);
2291: } else {
2292: PetscCall(PetscOptionsGetString(((PetscObject)pc)->options, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr));
2293: PetscCall(PetscStrcmp(type, PCMAT, &algebraic));
2294: 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 : "");
2295: if (overlap != -1) {
2296: 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 : "");
2297: PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-%spc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", pcpre ? pcpre : "", overlap);
2298: }
2299: if (block || overlap != -1) algebraic = PETSC_TRUE;
2300: PetscCheck(!data->scaling || !algebraic, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "User-supplied scaling cannot be combined with an internally assembled auxiliary problem");
2301: if (algebraic) {
2302: PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is));
2303: PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1));
2304: PetscCall(ISSort(data->is));
2305: } else
2306: 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 : ""));
2307: }
2308: }
2309: }
2310: }
2311: if (PetscDefined(USE_DEBUG)) {
2312: if (data->is) PetscCall(ISDuplicate(data->is, &dis));
2313: if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux));
2314: }
2315: if (data->is || (ismatis && data->N > 1)) {
2316: if (ismatis) {
2317: PetscCall(MatISGetLocalMat(P, &N));
2318: PetscCall(PetscObjectTypeCompareAny((PetscObject)N, &flg, MATSEQBAIJ, MATSEQSBAIJ, ""));
2319: PetscCall(MatISRestoreLocalMat(P, &N));
2320: PetscCall(MatConvert(P, flg ? MATMPIBAIJ : MATMPIAIJ, MAT_INITIAL_MATRIX, &C));
2321: PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr));
2322: PetscCall(PetscObjectReference((PetscObject)P));
2323: PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C));
2324: std::swap(C, P);
2325: PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n));
2326: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc));
2327: PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0]));
2328: PetscCall(ISDestroy(&loc));
2329: /* the auxiliary Mat is _not_ the local Neumann matrix */
2330: /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */
2331: data->Neumann = PETSC_BOOL3_FALSE;
2332: structure = SAME_NONZERO_PATTERN;
2333: } else {
2334: is[0] = data->is;
2335: if (algebraic || ctx) subdomains = PETSC_TRUE;
2336: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr));
2337: if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre);
2338: if (PetscBool3ToBool(data->Neumann)) {
2339: PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_block_splitting and -%spc_hpddm_has_neumann", pcpre ? pcpre : "", pcpre ? pcpre : "");
2340: 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 : "");
2341: 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 : "");
2342: }
2343: if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN;
2344: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc));
2345: }
2346: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2347: PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */
2348: if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */
2349: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : ""));
2350: PetscCall(PetscOptionsSetValue(((PetscObject)pc)->options, prefix, MatStructures[structure]));
2351: }
2352: flg = PETSC_FALSE;
2353: if (data->share) {
2354: data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */
2355: if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : ""));
2356: else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n"));
2357: else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n"));
2358: else if (!algebraic && structure != SAME_NONZERO_PATTERN)
2359: 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]));
2360: else {
2361: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHTOOL, &flg));
2362: if (flg) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since Pmat is of type MATHTOOL\n"));
2363: else data->share = PETSC_TRUE;
2364: }
2365: if (!data->share) {
2366: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_share_sub_ksp", pcpre ? pcpre : ""));
2367: PetscCall(PetscOptionsClearValue(((PetscObject)pc)->options, prefix));
2368: }
2369: }
2370: if (!ismatis) {
2371: if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted));
2372: else unsorted = is[0];
2373: }
2374: if ((ctx || data->N > 1) && (data->aux || ismatis || algebraic)) {
2375: PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level");
2376: PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2377: if (ismatis) {
2378: /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */
2379: PetscCall(MatIncreaseOverlap(P, 1, is, 1));
2380: PetscCall(ISDestroy(&data->is));
2381: data->is = is[0];
2382: } else {
2383: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_TRUE));
2384: if (!ctx && overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private));
2385: }
2386: if (algebraic && overlap == -1) {
2387: PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block));
2388: if (block) {
2389: PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux));
2390: PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr));
2391: }
2392: } else if (!ctx) {
2393: if (PetscBool3ToBool(data->Neumann)) sub = &data->aux;
2394: else {
2395: PetscBool flg;
2397: if (overlap != -1) {
2398: Harmonic h;
2399: Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */
2400: IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */
2401: const PetscInt *i[2], bs = P->cmap->bs; /* with a GEVP: [ A_00 A_01 ] */
2402: PetscInt n[2], location; /* [ A_10 A_11 A_12 ] */
2403: std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */
2405: do {
2406: PetscCall(ISDuplicate(data->is, ov));
2407: if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1));
2408: PetscCall(ISDuplicate(ov[0], ov + 1));
2409: PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1));
2410: PetscCall(ISGetLocalSize(ov[0], n));
2411: PetscCall(ISGetLocalSize(ov[1], n + 1));
2412: flg = PetscBool(n[0] == n[1] && n[0] != P->rmap->n);
2413: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
2414: if (flg) {
2415: PetscCall(ISDestroy(ov));
2416: PetscCall(ISDestroy(ov + 1));
2417: PetscCheck(--overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No oversampling possible");
2418: PetscCall(PetscInfo(pc, "Supplied -%spc_hpddm_harmonic_overlap parameter is too large, it has been decreased to %" PetscInt_FMT "\n", pcpre ? pcpre : "", overlap));
2419: } else break;
2420: } while (1);
2421: PetscCall(PetscNew(&h));
2422: h->ksp = nullptr;
2423: PetscCall(PetscCalloc1(2, &h->A));
2424: PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-eps_nev", &flg));
2425: if (!flg) {
2426: PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-svd_nsv", &flg));
2427: if (!flg) PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, prefix, "-svd_threshold_relative", &flg));
2428: } else flg = PETSC_FALSE;
2429: PetscCall(ISSort(ov[0]));
2430: if (!flg) PetscCall(ISSort(ov[1]));
2431: PetscCall(PetscCalloc1(5, &h->is));
2432: PetscCall(MatCreateSubMatrices(P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */
2433: for (PetscInt j = 0; j < 2; ++j) PetscCall(ISGetIndices(ov[j], i + j));
2434: v[1].reserve((n[1] - n[0]) / bs);
2435: for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */
2436: PetscCall(ISLocate(ov[0], i[1][j], &location));
2437: if (location < 0) v[1].emplace_back(j / bs);
2438: }
2439: if (!flg) {
2440: h->A[1] = a[0];
2441: PetscCall(PetscObjectReference((PetscObject)h->A[1]));
2442: v[0].reserve((n[0] - P->rmap->n) / bs);
2443: for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */
2444: PetscCall(ISLocate(loc, i[1][j], &location));
2445: if (location < 0) {
2446: PetscCall(ISLocate(ov[0], i[1][j], &location));
2447: if (location >= 0) v[0].emplace_back(j / bs);
2448: }
2449: }
2450: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows));
2451: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4));
2452: PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */
2453: PetscCall(ISDestroy(&rows));
2454: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows));
2455: PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */
2456: PetscCall(ISDestroy(&rows));
2457: v[0].clear();
2458: PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3));
2459: PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2));
2460: }
2461: v[0].reserve((n[0] - P->rmap->n) / bs);
2462: for (PetscInt j = 0; j < n[0]; j += bs) {
2463: PetscCall(ISLocate(loc, i[0][j], &location));
2464: if (location < 0) v[0].emplace_back(j / bs);
2465: }
2466: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows));
2467: for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j));
2468: if (flg) {
2469: PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &stride));
2470: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols));
2471: PetscCall(MatCreateSubMatrix(a[0], stride, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */
2472: PetscCall(ISDestroy(&cols));
2473: PetscCall(ISDestroy(&stride));
2474: PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg));
2475: if (flg) { /* initial Pmat was MATSBAIJ, convert back to the same format since this submatrix is square */
2476: PetscCall(MatSetOption(A0, MAT_SYMMETRIC, PETSC_TRUE));
2477: PetscCall(MatConvert(A0, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &A0));
2478: }
2479: flg = PETSC_TRUE;
2480: PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2));
2481: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols));
2482: PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */
2483: PetscCall(ISDestroy(&cols));
2484: }
2485: PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride));
2486: PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is));
2487: PetscCall(ISDestroy(&stride));
2488: PetscCall(ISDestroy(&rows));
2489: PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1));
2490: if (subdomains) {
2491: if (!data->levels[0]->pc) {
2492: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc));
2493: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2494: PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix));
2495: PetscCall(PCSetOperators(data->levels[0]->pc, A, P));
2496: }
2497: PetscCall(PCSetType(data->levels[0]->pc, PCASM));
2498: if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc));
2499: PetscCall(PCSetModifySubMatrices(data->levels[0]->pc, pc->modifysubmatrices, pc->modifysubmatricesP));
2500: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE));
2501: if (!flg) ++overlap;
2502: if (data->share) {
2503: PetscInt n = -1;
2505: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp));
2506: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2507: if (flg) {
2508: h->ksp = ksp[0];
2509: PetscCall(PetscObjectReference((PetscObject)h->ksp));
2510: }
2511: }
2512: }
2513: if (!h->ksp) {
2514: PetscBool share = data->share;
2516: PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp));
2517: PetscCall(KSPSetType(h->ksp, KSPPREONLY));
2518: PetscCall(KSPSetOperators(h->ksp, A0, A0));
2519: do {
2520: if (!data->share) {
2521: share = PETSC_FALSE;
2522: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_"));
2523: PetscCall(KSPSetOptionsPrefix(h->ksp, prefix));
2524: PetscCall(KSPSetFromOptions(h->ksp));
2525: } else {
2526: MatSolverType type;
2528: PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp[0]->pc, &data->share, PCLU, PCCHOLESKY, ""));
2529: if (data->share) {
2530: PetscCall(PCFactorGetMatSolverType(ksp[0]->pc, &type));
2531: if (!type) {
2532: if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(ksp[0]->pc, MATSOLVERMUMPS));
2533: else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(ksp[0]->pc, MATSOLVERMKL_PARDISO));
2534: else data->share = PETSC_FALSE;
2535: if (data->share) PetscCall(PCSetFromOptions(ksp[0]->pc));
2536: } else {
2537: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share));
2538: if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share));
2539: }
2540: if (data->share) {
2541: std::tuple<KSP, IS, Vec[2]> *p;
2543: PetscCall(PCFactorGetMatrix(ksp[0]->pc, &A));
2544: PetscCall(MatFactorSetSchurIS(A, h->is[4]));
2545: PetscCall(KSPSetUp(ksp[0]));
2546: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : ""));
2547: PetscCall(KSPSetOptionsPrefix(h->ksp, prefix));
2548: PetscCall(KSPSetFromOptions(h->ksp));
2549: PetscCall(PCSetType(h->ksp->pc, PCSHELL));
2550: PetscCall(PetscNew(&p));
2551: std::get<0>(*p) = ksp[0];
2552: PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p)));
2553: PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1));
2554: PetscCall(PCShellSetContext(h->ksp->pc, p));
2555: PetscCall(PCShellSetApply(h->ksp->pc, PCApply_Schur));
2556: PetscCall(PCShellSetApplyTranspose(h->ksp->pc, PCApply_Schur<Vec, true>));
2557: PetscCall(PCShellSetMatApply(h->ksp->pc, PCApply_Schur<Mat>));
2558: PetscCall(PCShellSetDestroy(h->ksp->pc, PCDestroy_Schur));
2559: }
2560: }
2561: if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n"));
2562: }
2563: } 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 */
2564: }
2565: PetscCall(ISDestroy(ov));
2566: PetscCall(ISDestroy(ov + 1));
2567: if (overlap == 1 && subdomains && flg) {
2568: *subA = A0;
2569: sub = subA;
2570: } else PetscCall(MatDestroy(&A0));
2571: PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux));
2572: PetscCall(MatSetVecType(data->aux, h->A[0]->defaultvectype));
2573: PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */
2574: PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr));
2575: PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Harmonic));
2576: PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Harmonic));
2577: PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE));
2578: PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AtB, nullptr, MatProduct_AtB_Harmonic, nullptr, MATDENSE, MATDENSE));
2579: PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Harmonic));
2580: PetscCall(MatDestroySubMatrices(1, &a));
2581: }
2582: if (overlap != 1 || !subdomains) {
2583: PetscCall(MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub));
2584: if (ismatis) {
2585: PetscCall(MatISGetLocalMat(C, &N));
2586: PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg));
2587: if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub));
2588: PetscCall(MatISRestoreLocalMat(C, &N));
2589: }
2590: }
2591: }
2592: }
2593: if (data->N > 1) {
2594: /* Vec holding the partition of unity */
2595: if (!data->levels[0]->D) {
2596: PetscCall(ISGetLocalSize(data->is, &n));
2597: PetscCall(VecCreate(PETSC_COMM_SELF, &data->levels[0]->D));
2598: PetscCall(VecSetSizes(data->levels[0]->D, n, n));
2599: PetscCall(VecSetType(data->levels[0]->D, A->defaultvectype));
2600: }
2601: if (data->share && overlap == -1) {
2602: Mat D;
2603: IS perm = nullptr;
2604: PetscInt size = -1;
2606: if (!data->levels[0]->pc) {
2607: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : ""));
2608: PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc));
2609: PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix));
2610: PetscCall(PCSetOperators(data->levels[0]->pc, A, P));
2611: }
2612: PetscCall(PCSetType(data->levels[0]->pc, PCASM));
2613: if (!ctx) {
2614: if (!data->levels[0]->pc->setupcalled) {
2615: IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */
2617: PetscCall(ISDuplicate(is[0], &sorted));
2618: PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc));
2619: PetscCall(PetscObjectDereference((PetscObject)sorted));
2620: }
2621: PetscCall(PCSetFromOptions(data->levels[0]->pc));
2622: PetscCall(PCSetModifySubMatrices(data->levels[0]->pc, pc->modifysubmatrices, pc->modifysubmatricesP));
2623: if (block) {
2624: PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm));
2625: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic));
2626: } else PetscCall(PCSetUp(data->levels[0]->pc));
2627: PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp));
2628: if (size != 1) {
2629: data->share = PETSC_FALSE;
2630: PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size);
2631: PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n"));
2632: PetscCall(ISDestroy(&unsorted));
2633: unsorted = is[0];
2634: } else {
2635: const char *matpre;
2636: PetscBool cmp[4];
2638: if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm));
2639: if (perm) { /* unsorted input IS */
2640: if (!PetscBool3ToBool(data->Neumann) && !block) {
2641: PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */
2642: PetscCall(MatHeaderReplace(sub[0], &D));
2643: }
2644: if (data->B) { /* see PCHPDDMSetRHSMat() */
2645: PetscCall(MatPermute(data->B, perm, perm, &D));
2646: PetscCall(MatHeaderReplace(data->B, &D));
2647: }
2648: PetscCall(ISDestroy(&perm));
2649: }
2650: PetscCall(KSPGetOperators(ksp[0], subA, subA + 1));
2651: PetscCall(PetscObjectReference((PetscObject)subA[0]));
2652: PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D));
2653: PetscCall(MatGetOptionsPrefix(subA[1], &matpre));
2654: PetscCall(MatSetOptionsPrefix(D, matpre));
2655: PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp));
2656: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1));
2657: if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2));
2658: else cmp[2] = PETSC_FALSE;
2659: if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3));
2660: else cmp[3] = PETSC_FALSE;
2661: 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);
2662: if (!cmp[0] && !cmp[2]) {
2663: if (!block) {
2664: if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckMatStructure_Private(pc, D, C));
2665: PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN));
2666: } else {
2667: structure = DIFFERENT_NONZERO_PATTERN;
2668: PetscCall(MatAXPY(D, 1.0, data->aux, structure));
2669: }
2670: } else {
2671: Mat mat[2];
2673: if (cmp[0]) {
2674: PetscCall(MatNormalGetMat(D, mat));
2675: PetscCall(MatNormalGetMat(C, mat + 1));
2676: } else {
2677: PetscCall(MatNormalHermitianGetMat(D, mat));
2678: PetscCall(MatNormalHermitianGetMat(C, mat + 1));
2679: }
2680: PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN));
2681: }
2682: PetscCall(MatPropagateSymmetryOptions(C, D));
2683: PetscCall(MatDestroy(&C));
2684: C = D;
2685: /* swap pointers so that variables stay consistent throughout PCSetUp() */
2686: std::swap(C, data->aux);
2687: std::swap(uis, data->is);
2688: swap = PETSC_TRUE;
2689: }
2690: }
2691: }
2692: }
2693: if (ctx) {
2694: PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data;
2695: PC s;
2696: Mat A00, P00, A01 = nullptr, A10, N, b[4];
2697: IS sorted, is[2], *is_00;
2698: MatSolverType type;
2699: std::pair<PC, Vec[2]> *p;
2701: n = -1;
2702: PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp));
2703: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2704: PetscCall(KSPGetOperators(ksp[0], subA, subA + 1));
2705: PetscCall(ISGetLocalSize(data_00->is, &n));
2706: if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) {
2707: PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr));
2708: PetscCall(ISGetLocalSize(*is_00, &n));
2709: 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);
2710: } else is_00 = &data_00->is;
2711: PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */
2712: std::swap(C, data->aux);
2713: std::swap(uis, data->is);
2714: swap = PETSC_TRUE;
2715: PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, nullptr));
2716: std::get<1>(*ctx)[1] = A10;
2717: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg));
2718: if (flg) PetscCall(MatTransposeGetMat(A10, &A01));
2719: else {
2720: PetscBool flg;
2722: PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg));
2723: if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01));
2724: }
2725: 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 */
2726: PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */
2727: if (!A01) {
2728: PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2729: PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub));
2730: b[2] = sub[0];
2731: PetscCall(PetscObjectReference((PetscObject)sub[0]));
2732: PetscCall(MatDestroySubMatrices(1, &sub));
2733: PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg));
2734: A10 = nullptr;
2735: if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10));
2736: else {
2737: PetscBool flg;
2739: PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg));
2740: if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10));
2741: }
2742: if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub));
2743: else {
2744: if (flg) PetscCall(MatCreateTranspose(b[2], b + 1));
2745: else PetscCall(MatCreateHermitianTranspose(b[2], b + 1));
2746: }
2747: } else {
2748: PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
2749: PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub));
2750: if (flg) PetscCall(MatCreateTranspose(*sub, b + 2));
2751: else PetscCall(MatCreateHermitianTranspose(*sub, b + 2));
2752: }
2753: if (A01 || !A10) {
2754: b[1] = sub[0];
2755: PetscCall(PetscObjectReference((PetscObject)sub[0]));
2756: }
2757: PetscCall(MatDestroySubMatrices(1, &sub));
2758: PetscCall(ISDestroy(&sorted));
2759: b[3] = data->aux;
2760: PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S));
2761: PetscCall(MatSchurComplementSetKSP(S, ksp[0]));
2762: if (data->N != 1) {
2763: PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */
2764: PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc));
2765: PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */
2766: s = data->levels[0]->pc;
2767: } else {
2768: is[0] = data->is;
2769: PetscCall(PetscObjectReference((PetscObject)is[0]));
2770: PetscCall(PetscObjectReference((PetscObject)b[3]));
2771: PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */
2772: data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */
2773: PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */
2774: PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc));
2775: PetscCall(ISDestroy(is));
2776: PetscCall(ISDestroy(&loc));
2777: s = pc;
2778: }
2779: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */
2780: PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp));
2781: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n);
2782: PetscCall(KSPGetPC(ksp[0], &inner));
2783: PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */
2784: b[0] = subA[0];
2785: 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 */
2786: if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3]));
2787: PetscCall(PetscObjectDereference((PetscObject)b[1]));
2788: PetscCall(PetscObjectDereference((PetscObject)b[2]));
2789: PetscCall(PCCreate(PETSC_COMM_SELF, &s));
2790: PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix));
2791: PetscCall(PCSetOptionsPrefix(inner, nullptr));
2792: PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE));
2793: PetscCall(PCSetType(s, PCLU));
2794: 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 */
2795: PetscCall(PCSetFromOptions(s));
2796: PetscCall(PCFactorGetMatSolverType(s, &type));
2797: PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg));
2798: PetscCall(MatGetLocalSize(P, &n, nullptr));
2799: if (flg || n == 0) {
2800: PetscCall(PCSetOperators(s, N, N));
2801: if (n) {
2802: PetscCall(PCFactorGetMatrix(s, b));
2803: PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix));
2804: n = -1;
2805: PetscCall(PetscOptionsGetInt(((PetscObject)pc)->options, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr));
2806: if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */
2807: PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */
2808: PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */
2809: }
2810: } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */
2811: } else {
2812: MatInfo info;
2814: PetscCall(MatGetInfo(b[3], MAT_LOCAL, &info));
2815: if (info.nz_used == 0.0) {
2816: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, b[3]->rmap->n, b[3]->cmap->n, 1, nullptr, b + 3));
2817: for (PetscInt i = 0; i < b[3]->rmap->n; ++i) PetscCall(MatSetValue(b[3], i, i, 0.0, INSERT_VALUES));
2818: PetscCall(MatAssemblyBegin(b[3], MAT_FINAL_ASSEMBLY));
2819: PetscCall(MatAssemblyEnd(b[3], MAT_FINAL_ASSEMBLY));
2820: PetscCall(MatNestSetSubMat(N, 1, 1, b[3]));
2821: PetscCall(PetscObjectDereference((PetscObject)b[3]));
2822: }
2823: PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b));
2824: PetscCall(PCSetOperators(s, N, *b));
2825: PetscCall(PetscObjectDereference((PetscObject)*b));
2826: PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, ""));
2827: if (flg) {
2828: 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 */
2829: if (info.nz_used == 0.0) {
2830: PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, ""));
2831: if (flg) {
2832: PetscCall(MatFactorGetSolverType(*b, &type));
2833: PetscCall(PetscStrcmp(type, MATSOLVERPETSC, &flg));
2834: if (flg) {
2835: PetscCall(PetscOptionsHasName(((PetscObject)s)->options, ((PetscObject)s)->prefix, "-pc_factor_mat_ordering_type", &flg));
2836: if (!flg) PetscCall(PCFactorSetMatOrderingType(s, MATORDERINGNATURAL));
2837: }
2838: }
2839: flg = PETSC_TRUE;
2840: }
2841: }
2842: }
2843: PetscCall(PetscNew(&p));
2844: p->first = s;
2845: if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1));
2846: else p->second[0] = p->second[1] = nullptr;
2847: PetscCall(PCShellSetContext(inner, p));
2848: PetscCall(PCShellSetApply(inner, PCApply_Nest));
2849: PetscCall(PCShellSetView(inner, PCView_Nest));
2850: PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest));
2851: PetscCall(PetscObjectDereference((PetscObject)N));
2852: if (!data) {
2853: PetscCall(MatDestroy(&S));
2854: PetscCall(ISDestroy(&unsorted));
2855: PetscCall(MatDestroy(&C));
2856: PetscCall(ISDestroy(&uis));
2857: PetscCall(PetscFree(ctx));
2858: if (PetscDefined(USE_DEBUG)) {
2859: PetscCall(ISDestroy(&dis));
2860: PetscCall(MatDestroy(&daux));
2861: }
2862: PetscFunctionReturn(PETSC_SUCCESS);
2863: }
2864: }
2865: if (!data->levels[0]->scatter) {
2866: PetscCall(MatCreateVecs(P, &xin, nullptr));
2867: if (ismatis) PetscCall(MatDestroy(&P));
2868: PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter));
2869: PetscCall(VecDestroy(&xin));
2870: }
2871: if (data->levels[0]->P) {
2872: /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */
2873: PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], !pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE));
2874: }
2875: if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>();
2876: if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr));
2877: else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr));
2878: /* HPDDM internal data structure */
2879: PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels));
2880: if (data->scaling && data->deflation) {
2881: const PetscScalar *array;
2882: PetscReal *d;
2884: PetscCall(PCHPDDMComputeScaling_Private(data, data->is, data->levels[0]->D));
2885: PetscCall(VecGetArrayRead(data->levels[0]->D, &array));
2886: /* HPDDM owns mutable real storage; in complex builds it is separate from the PETSc Vec */
2887: d = const_cast<PetscReal *>(data->levels[0]->P->getScaling());
2888: for (PetscInt i = 0; i < data->levels[0]->P->getDof(); ++i) d[i] = PetscRealPart(array[i]);
2889: PetscCall(VecRestoreArrayRead(data->levels[0]->D, &array));
2890: }
2891: if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr));
2892: /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */
2893: if (!ctx) {
2894: if (data->deflation || overlap != -1) weighted = data->aux;
2895: else if (!data->B) {
2896: PetscBool cmp;
2898: PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted));
2899: PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, ""));
2900: if (cmp) flg = PETSC_FALSE;
2901: PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D));
2902: /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */
2903: /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */
2904: PetscCall(MatPropagateSymmetryOptions(sub[0], weighted));
2905: if (PetscDefined(USE_DEBUG) && PetscBool3ToBool(data->Neumann)) {
2906: Mat *sub, A[2];
2907: PetscReal norm[2];
2909: PetscCall(MatCreateSubMatrices(P, 1, &data->is, &data->is, MAT_INITIAL_MATRIX, &sub));
2910: PetscCall(MatDiagonalScale(sub[0], data->levels[0]->D, data->levels[0]->D));
2911: 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 */
2912: PetscCall(MatConvert(weighted, MATSEQAIJ, MAT_INITIAL_MATRIX, A + 1));
2913: PetscCall(MatAXPY(A[0], -1.0, A[1], UNKNOWN_NONZERO_PATTERN));
2914: PetscCall(MatNorm(A[0], NORM_FROBENIUS, norm));
2915: if (norm[0]) {
2916: PetscCall(MatNorm(A[1], NORM_FROBENIUS, norm + 1));
2917: 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 : "");
2918: }
2919: PetscCall(MatDestroySubMatrices(1, &sub));
2920: for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(A + i));
2921: }
2922: } else weighted = data->B;
2923: } else weighted = nullptr;
2924: /* SLEPc is used inside the loaded symbol */
2925: 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));
2926: if (!ctx && data->share && overlap == -1) {
2927: Mat st[2];
2929: PetscCheck(ksp && ksp[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing shared subdomain KSP");
2930: PetscCall(KSPGetOperators(ksp[0], st, st + 1));
2931: PetscCall(MatCopy(subA[0], st[0], structure));
2932: if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN));
2933: PetscCall(PetscObjectDereference((PetscObject)subA[0]));
2934: }
2935: if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr));
2936: if (ismatis) PetscCall(MatISGetLocalMat(C, &N));
2937: else N = data->aux;
2938: if (!ctx) P = sub[0];
2939: else P = S;
2940: /* going through the grid hierarchy */
2941: for (n = 1; n < data->N; ++n) {
2942: if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr));
2943: /* method composed in the loaded symbol since there, SLEPc is used as well */
2944: PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels));
2945: if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr));
2946: }
2947: /* reset to NULL to avoid any faulty use */
2948: PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr));
2949: if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr));
2950: else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */
2951: for (n = 0; n < data->N - 1; ++n) {
2952: /* HPDDM internal work buffers */
2953: if (data->levels[n]->P) PetscCallCXX(data->levels[n]->P->start());
2954: }
2955: if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub));
2956: if (ismatis) data->is = nullptr;
2957: for (n = 0; n < data->N - 1 + (reused > 0); ++n) {
2958: if (data->levels[n]->P) {
2959: PC spc;
2961: /* force the PC to be PCSHELL to do the coarse grid corrections */
2962: PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE));
2963: PetscCall(KSPGetPC(data->levels[n]->ksp, &spc));
2964: PetscCall(PCSetType(spc, PCSHELL));
2965: PetscCall(PCShellSetContext(spc, data->levels[n]));
2966: PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell));
2967: PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell));
2968: PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell));
2969: PetscCall(PCShellSetApplyTranspose(spc, PCApplyTranspose_HPDDMShell));
2970: PetscCall(PCShellSetMatApplyTranspose(spc, PCMatApplyTranspose_HPDDMShell));
2971: if (ctx && n == 0) {
2972: Mat Amat, Pmat;
2973: PetscInt m, M;
2974: std::tuple<Mat, PetscSF, Vec[2]> *ctx;
2976: PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat));
2977: PetscCall(MatGetLocalSize(Pmat, &m, nullptr));
2978: PetscCall(MatGetSize(Pmat, &M, nullptr));
2979: PetscCall(PetscNew(&ctx));
2980: std::get<0>(*ctx) = S;
2981: std::get<1>(*ctx) = data->levels[n]->scatter;
2982: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat));
2983: PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Schur<false>));
2984: PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMult_Schur<true>));
2985: PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Schur));
2986: PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1));
2987: PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat));
2988: PetscCall(PetscObjectDereference((PetscObject)Amat));
2989: }
2990: PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell));
2991: if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc));
2992: if (n < reused) {
2993: PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE));
2994: PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE));
2995: }
2996: PetscCall(PCSetUp(spc));
2997: }
2998: }
2999: if (ctx) PetscCall(MatDestroy(&S));
3000: if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr));
3001: } else flg = reused ? PETSC_FALSE : PETSC_TRUE;
3002: if (!ismatis && subdomains) {
3003: if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner));
3004: else inner = data->levels[0]->pc;
3005: if (inner) {
3006: if (!inner->setupcalled) {
3007: PetscCall(PCSetType(inner, PCASM));
3008: if (data->scaling) PetscCall(PCASMSetType(inner, PC_ASM_WEIGHTED));
3009: }
3010: PetscCall(PCSetFromOptions(inner));
3011: PetscCall(PCSetModifySubMatrices(inner, pc->modifysubmatrices, pc->modifysubmatricesP));
3012: PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg));
3013: if (flg) {
3014: if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */
3015: IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */
3017: PetscCall(ISDuplicate(is[0], &sorted));
3018: PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc));
3019: PetscCall(PetscObjectDereference((PetscObject)sorted));
3020: }
3021: if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */
3022: PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr));
3023: PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic));
3024: PetscCall(PetscObjectDereference((PetscObject)P));
3025: }
3026: }
3027: }
3028: if (data->N > 1) {
3029: if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub));
3030: if (overlap == 1) PetscCall(MatDestroy(subA));
3031: }
3032: }
3033: PetscCall(ISDestroy(&loc));
3034: } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */
3035: if (requested != data->N + reused) {
3036: 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,
3037: data->N, pcpre ? pcpre : ""));
3038: 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 : "",
3039: data->N, pcpre ? pcpre : "", data->N));
3040: /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */
3041: for (n = data->N - 1; n < requested - 1; ++n) {
3042: if (data->levels[n]->P) {
3043: PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE));
3044: PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0]));
3045: PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1]));
3046: PetscCall(MatDestroy(data->levels[n]->V));
3047: PetscCall(MatDestroy(data->levels[n]->V + 1));
3048: PetscCall(MatDestroy(data->levels[n]->V + 2));
3049: PetscCall(VecDestroy(&data->levels[n]->D));
3050: PetscCall(PetscSFDestroy(&data->levels[n]->scatter));
3051: }
3052: }
3053: if (reused) {
3054: for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) {
3055: PetscCall(KSPDestroy(&data->levels[n]->ksp));
3056: PetscCall(PCDestroy(&data->levels[n]->pc));
3057: }
3058: }
3059: 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,
3060: data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N, pcpre ? pcpre : "", data->N);
3061: }
3062: /* these solvers are created after PCSetFromOptions() is called */
3063: if (pc->setfromoptionscalled) {
3064: for (n = 0; n < data->N; ++n) {
3065: if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp));
3066: if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc));
3067: }
3068: pc->setfromoptionscalled = 0;
3069: }
3070: data->N += reused;
3071: if (data->share && swap) {
3072: /* swap back pointers so that variables follow the user-provided numbering */
3073: std::swap(C, data->aux);
3074: std::swap(uis, data->is);
3075: PetscCall(MatDestroy(&C));
3076: PetscCall(ISDestroy(&uis));
3077: }
3078: if (data->scaling) {
3079: if (!data->levels[0]->pc) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); /* single level: inspect directly the PC attached to the KSP */
3080: else inner = data->levels[0]->pc;
3081: PetscCall(PetscObjectTypeCompare((PetscObject)inner, PCASM, &flg));
3082: if (flg) {
3083: PCASMType type;
3085: PetscCall(PCASMGetType(inner, &type));
3086: flg = (type == PC_ASM_WEIGHTED ? PETSC_TRUE : PETSC_FALSE);
3087: }
3088: if (flg) { /* forward the scaling to PC_ASM_WEIGHTED, otherwise it only scales the coarse space */
3089: IS *is;
3091: PetscCall(PCSetUp(inner));
3092: PetscCall(PCASMGetLocalSubdomains(inner, &n, &is, nullptr));
3093: PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "User-supplied scaling requires one overlapping subdomain per process");
3094: PetscCall(PCASMGetLocalSubmatrices(inner, nullptr, &sub));
3095: PetscCall(MatCreateVecs(sub[0], &v, nullptr));
3096: PetscCall(PCHPDDMComputeScaling_Private(data, is[0], v));
3097: PetscCall(PCASMWeightedSetScaling(inner, 1, &v));
3098: PetscCall(VecDestroy(&v));
3099: }
3100: }
3101: if (algebraic) PetscCall(MatDestroy(&data->aux));
3102: if (unsorted && unsorted != is[0]) {
3103: PetscCall(ISCopy(unsorted, data->is));
3104: PetscCall(ISDestroy(&unsorted));
3105: }
3106: if (PetscDefined(USE_DEBUG)) {
3107: 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);
3108: if (data->is) {
3109: PetscCall(ISEqualUnsorted(data->is, dis, &flg));
3110: PetscCall(ISDestroy(&dis));
3111: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal");
3112: }
3113: 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);
3114: if (data->aux) {
3115: PetscCall(MatMultEqual(data->aux, daux, 20, &flg));
3116: PetscCall(MatDestroy(&daux));
3117: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal");
3118: }
3119: }
3120: PetscFunctionReturn(PETSC_SUCCESS);
3121: }
3123: /*@
3124: PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type.
3126: Collective
3128: Input Parameters:
3129: + pc - preconditioner context
3130: - type - coarse correction type, see `PCHPDDMCoarseCorrectionType`
3132: Options Database Key:
3133: . -pc_hpddm_coarse_correction (deflated|additive|balanced|none|deflated_reversed) - type of coarse correction to apply
3135: Level: intermediate
3137: .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
3138: @*/
3139: PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type)
3140: {
3141: PetscFunctionBegin;
3144: PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type));
3145: PetscFunctionReturn(PETSC_SUCCESS);
3146: }
3148: /*@
3149: PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type.
3151: Input Parameter:
3152: . pc - preconditioner context
3154: Output Parameter:
3155: . type - coarse correction type, see `PCHPDDMCoarseCorrectionType`
3157: Level: intermediate
3159: .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType`
3160: @*/
3161: PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type)
3162: {
3163: PetscFunctionBegin;
3165: if (type) {
3166: PetscAssertPointer(type, 2);
3167: PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type));
3168: }
3169: PetscFunctionReturn(PETSC_SUCCESS);
3170: }
3172: static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type)
3173: {
3174: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3176: PetscFunctionBegin;
3177: data->correction = type;
3178: PetscFunctionReturn(PETSC_SUCCESS);
3179: }
3181: static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type)
3182: {
3183: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3185: PetscFunctionBegin;
3186: *type = data->correction;
3187: PetscFunctionReturn(PETSC_SUCCESS);
3188: }
3190: /*@
3191: PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared.
3193: Input Parameters:
3194: + pc - preconditioner context
3195: - share - whether the `KSP` should be shared or not
3197: Note:
3198: This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped
3199: when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`.
3201: Level: advanced
3203: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()`
3204: @*/
3205: PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share)
3206: {
3207: PetscFunctionBegin;
3209: PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share));
3210: PetscFunctionReturn(PETSC_SUCCESS);
3211: }
3213: /*@
3214: PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared.
3216: Input Parameter:
3217: . pc - preconditioner context
3219: Output Parameter:
3220: . share - whether the `KSP` is shared or not
3222: Note:
3223: 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
3224: when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`.
3226: Level: advanced
3228: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()`
3229: @*/
3230: PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share)
3231: {
3232: PetscFunctionBegin;
3234: if (share) {
3235: PetscAssertPointer(share, 2);
3236: PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share));
3237: }
3238: PetscFunctionReturn(PETSC_SUCCESS);
3239: }
3241: static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share)
3242: {
3243: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3245: PetscFunctionBegin;
3246: data->share = share;
3247: PetscFunctionReturn(PETSC_SUCCESS);
3248: }
3250: static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share)
3251: {
3252: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3254: PetscFunctionBegin;
3255: *share = data->share;
3256: PetscFunctionReturn(PETSC_SUCCESS);
3257: }
3259: /*@
3260: PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator.
3262: Logically Collective; `U` must be `NULL` on all processes or on none
3264: Input Parameters:
3265: + pc - preconditioner context
3266: . is - index set of the local deflation matrix, or `NULL` to use the stored one
3267: - U - deflation sequential matrix stored as a `MATSEQDENSE`, or `NULL` to remove the stored one
3269: Level: advanced
3271: Notes:
3272: The `is` argument is shared with `PCHPDDMSetScaling()`. Both routines can be called in any order, or
3273: only one of them. If `is` is not `NULL`, it replaces the stored one, and a stored scaling is kept and
3274: must have the same size. If `is` is `NULL`, the stored one, which must exist, is used. When `U` is
3275: `NULL`, `is` must be `NULL` too, and the stored one is discarded unless a scaling still uses it; remove
3276: the deflation space before `PCSetUp()`, or call `PCReset()`, which discards the deflation space, the
3277: scaling, and their index set.
3279: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetScaling()`, `PCDeflationSetSpace()`, `PCMGSetRestriction()`
3280: @*/
3281: PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U)
3282: {
3283: PetscFunctionBegin;
3287: PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U));
3288: PetscFunctionReturn(PETSC_SUCCESS);
3289: }
3291: static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U)
3292: {
3293: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3294: Vec scaling = data->scaling;
3296: PetscFunctionBegin;
3297: if (U) {
3298: if (!is) is = data->is;
3299: PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No stored IS, supply one to PCHPDDMSetDeflationMat() or PCHPDDMSetScaling()");
3300: if (scaling) {
3301: PetscInt n, m;
3303: PetscCall(ISGetLocalSize(is, &n));
3304: PetscCall(VecGetSize(scaling, &m));
3305: 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);
3306: PetscCall(PetscObjectReference((PetscObject)scaling));
3307: }
3308: /* the scaling shares the IS, but PCHPDDMSetAuxiliaryMat_Private() destroys it if it resets the PC */
3309: PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE));
3310: PetscCall(VecDestroy(&data->scaling));
3311: data->scaling = scaling;
3312: } else {
3313: PetscCheck(!is, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Pass a NULL IS to remove the deflation space");
3314: PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Remove the deflation space before PCSetUp() or after PCReset()");
3315: if (data->deflation) {
3316: PetscCall(MatDestroy(&data->aux));
3317: data->deflation = PETSC_FALSE;
3318: }
3319: if (!data->aux && !scaling) PetscCall(ISDestroy(&data->is));
3320: }
3321: PetscFunctionReturn(PETSC_SUCCESS);
3322: }
3324: /*@
3325: PCHPDDMSetScaling - Sets a diagonal scaling for supplied deflation vectors and, by default, the fine-level `PCASM`
3327: Logically Collective; `scaling` must be `NULL` on all processes or on none
3329: Input Parameters:
3330: + pc - preconditioner context
3331: . is - index set of the local scaling, or `NULL` to use the stored one
3332: - scaling - local diagonal scaling, or `NULL` to remove the stored one and restore the default Boolean scaling before setup
3334: Level: advanced
3336: Notes:
3337: The `is` argument is shared with `PCHPDDMSetDeflationMat()`. Both routines can be called in any
3338: order before `PCSetUp()`, and the second call can provide `NULL` for `is`.
3340: Setting a scaling sets the default `PCASMType` of the fine-level `PCASM` to `PC_ASM_WEIGHTED`. It can be changed with
3341: `-pc_hpddm_levels_1_pc_asm_type`, or `-pc_hpddm_coarse_pc_asm_type` with a single level; the scaling then only applies
3342: to the coarse correction, if any.
3344: Currently, this cannot be used when `PCHPDDM` generates the second level with `PCHPDDMSetAuxiliaryMat()`.
3346: .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetDeflationMat()`, `PCASMWeightedSetScaling()`, `PCHPDDMSetCoarseCorrectionType()`
3347: @*/
3348: PetscErrorCode PCHPDDMSetScaling(PC pc, IS is, Vec scaling)
3349: {
3350: PetscFunctionBegin;
3354: PetscTryMethod(pc, "PCHPDDMSetScaling_C", (PC, IS, Vec), (pc, is, scaling));
3355: PetscFunctionReturn(PETSC_SUCCESS);
3356: }
3358: static PetscErrorCode PCHPDDMSetScaling_HPDDM(PC pc, IS is, Vec scaling)
3359: {
3360: PC_HPDDM *data = (PC_HPDDM *)pc->data;
3362: PetscFunctionBegin;
3363: PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Set the scaling before PCSetUp() or after PCReset()");
3364: if (scaling) {
3365: const PetscScalar *array;
3366: PetscInt n, m;
3367: PetscMPIInt size;
3368: PetscBool valid = PETSC_TRUE;
3370: PetscCheck(!data->aux || data->deflation, PETSC_COMM_SELF, PETSC_ERR_SUP, "User-supplied scaling cannot be combined with an auxiliary matrix supplied with PCHPDDMSetAuxiliaryMat()");
3371: if (!is) is = data->is;
3372: PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No stored IS, supply one to PCHPDDMSetScaling() or PCHPDDMSetDeflationMat()");
3373: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)scaling), &size));
3374: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Scaling must have a single-process communicator");
3375: PetscCall(VecGetSize(scaling, &n));
3376: PetscCall(ISGetLocalSize(is, &m));
3377: PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Scaling has size %" PetscInt_FMT ", expected %" PetscInt_FMT, n, m);
3378: if (data->aux && is != data->is) {
3379: PetscCall(MatGetLocalSize(data->aux, &n, nullptr));
3380: 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);
3381: }
3382: PetscCall(VecGetArrayRead(scaling, &array));
3383: for (PetscInt i = 0; (i < m) && valid; ++i) valid = PetscBool(PetscImaginaryPart(array[i]) == 0 && PetscRealPart(array[i]) >= 0 && !PetscIsInfOrNanScalar(array[i]));
3384: PetscCall(VecRestoreArrayRead(scaling, &array));
3385: PetscCheck(valid, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scaling must contain finite, real, nonnegative entries");
3386: PetscCall(PetscObjectReference((PetscObject)scaling));
3387: PetscCall(PetscObjectReference((PetscObject)is));
3388: PetscCall(ISDestroy(&data->is));
3389: data->is = is;
3390: } else {
3391: PetscCheck(!is, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Pass a NULL IS to remove previously stored scaling");
3392: if (!data->aux) PetscCall(ISDestroy(&data->is));
3393: }
3394: PetscCall(VecDestroy(&data->scaling));
3395: data->scaling = scaling;
3396: PetscFunctionReturn(PETSC_SUCCESS);
3397: }
3399: PetscErrorCode HPDDMLoadDL_Private(PetscBool *found)
3400: {
3401: PetscBool flg;
3402: char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN];
3404: PetscFunctionBegin;
3405: PetscAssertPointer(found, 1);
3406: PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir)));
3407: PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr));
3408: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir));
3409: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3410: #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */
3411: if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */
3412: PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir)));
3413: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir));
3414: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3415: }
3416: #endif
3417: if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */
3418: PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg));
3419: 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 */
3420: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir));
3421: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3422: PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir);
3423: PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib));
3424: PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */
3425: PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found));
3426: }
3427: }
3428: PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib);
3429: PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib));
3430: PetscFunctionReturn(PETSC_SUCCESS);
3431: }
3433: /*MC
3434: PCHPDDM - Interface with the HPDDM library.
3436: This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`.
3437: It may be viewed as an alternative to spectral
3438: AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020`
3440: 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`).
3442: For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using
3443: `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`.
3445: Options Database Keys:
3446: + -pc_hpddm_define_subdomains (true|false) - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()`
3447: (not relevant with an unassembled Pmat)
3448: . -pc_hpddm_has_neumann (true|false) - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()`
3449: - -pc_hpddm_coarse_correction type - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` default is `deflated`
3451: Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes.
3452: .vb
3453: -pc_hpddm_levels_%d_pc_
3454: -pc_hpddm_levels_%d_ksp_
3455: -pc_hpddm_levels_%d_eps_
3456: -pc_hpddm_levels_%d_p
3457: -pc_hpddm_levels_%d_mat_type
3458: -pc_hpddm_coarse_
3459: -pc_hpddm_coarse_p
3460: -pc_hpddm_coarse_mat_type
3461: -pc_hpddm_coarse_mat_filter
3462: .ve
3464: 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
3465: -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij` will use 10 deflation vectors per subdomain on the fine "level 1",
3466: aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2",
3467: and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`).
3469: 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.
3471: Level: intermediate
3473: Notes:
3474: This preconditioner requires that PETSc is built with SLEPc (`--download-slepc`).
3476: By default, the underlying concurrent eigenproblems
3477: are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf.
3478: {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As
3479: 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
3480: -pc_hpddm_levels_1_st_type sinvert`. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in
3481: SLEPc documentation since they are specific to `PCHPDDM`.
3482: .vb
3483: -pc_hpddm_levels_1_st_share_sub_ksp
3484: -pc_hpddm_levels_%d_eps_threshold_absolute
3485: -pc_hpddm_levels_1_eps_use_inertia
3486: .ve
3488: The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is
3489: used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse
3490: 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
3491: 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
3492: 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
3493: to supply `-pc_hpddm_levels_1_eps_nev`. This last option also only applies to the fine-level (N = 1) eigensolver.
3495: See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent`
3497: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`,
3498: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetScaling()`, `PCHPDDMSetSTShareSubKSP()`,
3499: `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()`
3500: M*/
3501: PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc)
3502: {
3503: PC_HPDDM *data;
3504: PetscBool found;
3506: PetscFunctionBegin;
3507: if (!loadedSym) {
3508: PetscCall(HPDDMLoadDL_Private(&found));
3509: if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym));
3510: }
3511: PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc");
3512: PetscCall(PetscNew(&data));
3513: pc->data = data;
3514: data->Neumann = PETSC_BOOL3_UNKNOWN;
3515: pc->ops->reset = PCReset_HPDDM;
3516: pc->ops->destroy = PCDestroy_HPDDM;
3517: pc->ops->setfromoptions = PCSetFromOptions_HPDDM;
3518: pc->ops->setup = PCSetUp_HPDDM;
3519: pc->ops->apply = PCApply_HPDDM<false>;
3520: pc->ops->matapply = PCMatApply_HPDDM<false>;
3521: pc->ops->applytranspose = PCApply_HPDDM<true>;
3522: pc->ops->matapplytranspose = PCMatApply_HPDDM<true>;
3523: pc->ops->view = PCView_HPDDM;
3524: pc->ops->presolve = PCPreSolve_HPDDM;
3526: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM));
3527: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM));
3528: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM));
3529: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM));
3530: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM));
3531: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM));
3532: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM));
3533: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM));
3534: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetScaling_C", PCHPDDMSetScaling_HPDDM));
3535: PetscFunctionReturn(PETSC_SUCCESS);
3536: }
3538: /*@
3539: PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`.
3541: Level: developer
3543: .seealso: [](ch_ksp), `PetscInitialize()`
3544: @*/
3545: PetscErrorCode PCHPDDMInitializePackage(void)
3546: {
3547: char ename[32];
3549: PetscFunctionBegin;
3550: if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
3551: PCHPDDMPackageInitialized = PETSC_TRUE;
3552: PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage));
3553: /* general events registered once during package initialization */
3554: /* some of these events are not triggered in libpetsc, */
3555: /* but rather directly in libhpddm_petsc, */
3556: /* which is in charge of performing the following operations */
3558: /* domain decomposition structure from Pmat sparsity pattern */
3559: PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc));
3560: /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */
3561: PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP));
3562: /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */
3563: PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP));
3564: /* next level construction using PtAP and PtBP (not triggered in libpetsc) */
3565: PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next));
3566: static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high");
3567: for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) {
3568: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i));
3569: /* events during a PCSetUp() at level #i _except_ the assembly */
3570: /* of the Galerkin operator of the coarser level #(i + 1) */
3571: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1]));
3572: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i));
3573: /* events during a PCApply() at level #i _except_ */
3574: /* the KSPSolve() of the coarser level #(i + 1) */
3575: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1]));
3576: }
3577: PetscFunctionReturn(PETSC_SUCCESS);
3578: }
3580: /*@
3581: PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`.
3583: Level: developer
3585: .seealso: [](ch_ksp), `PetscFinalize()`
3586: @*/
3587: PetscErrorCode PCHPDDMFinalizePackage(void)
3588: {
3589: PetscFunctionBegin;
3590: PCHPDDMPackageInitialized = PETSC_FALSE;
3591: PetscFunctionReturn(PETSC_SUCCESS);
3592: }
3594: static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y)
3595: {
3596: Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */
3597: /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */
3598: Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */
3599: /* [ A_10 A_11 ] R_1^T A_12 x [ ] */
3600: PetscFunctionBegin;
3601: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3602: PetscCall(VecSet(h->v, 0.0));
3603: PetscCall(VecGetSubVector(h->v, h->is[0], &sub));
3604: PetscCall(MatMult(h->A[0], x, sub));
3605: PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub));
3606: PetscCall(KSPSolve(h->ksp, h->v, h->v));
3607: PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y));
3608: PetscFunctionReturn(PETSC_SUCCESS);
3609: }
3611: static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x)
3612: {
3613: Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */
3614: Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */
3616: PetscFunctionBegin;
3617: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3618: PetscCall(VecSet(h->v, 0.0));
3619: PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y));
3620: PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v));
3621: PetscCall(VecGetSubVector(h->v, h->is[0], &sub));
3622: PetscCall(MatMultTranspose(h->A[0], sub, x));
3623: PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub));
3624: PetscFunctionReturn(PETSC_SUCCESS);
3625: }
3627: static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *)
3628: {
3629: Harmonic h;
3630: Mat A, B;
3631: Vec a, b;
3633: PetscFunctionBegin;
3634: PetscCall(MatShellGetContext(S, static_cast<void *>(&h)));
3635: PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A));
3636: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B));
3637: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3638: PetscCall(MatDenseGetColumnVecRead(A, i, &a));
3639: PetscCall(MatDenseGetColumnVecWrite(B, i, &b));
3640: PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a));
3641: PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b));
3642: PetscCall(MatDenseRestoreColumnVecRead(A, i, &a));
3643: }
3644: PetscCall(MatDestroy(&A));
3645: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A));
3646: PetscCall(KSPMatSolve(h->ksp, B, A));
3647: PetscCall(MatDestroy(&B));
3648: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3649: PetscCall(MatDenseGetColumnVecRead(A, i, &a));
3650: PetscCall(MatDenseGetColumnVecWrite(Y, i, &b));
3651: PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b));
3652: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b));
3653: PetscCall(MatDenseRestoreColumnVecRead(A, i, &a));
3654: }
3655: PetscCall(MatDestroy(&A));
3656: PetscFunctionReturn(PETSC_SUCCESS);
3657: }
3659: static PetscErrorCode MatProduct_AtB_Harmonic(Mat S, Mat Y, Mat X, void *)
3660: {
3661: Harmonic h;
3662: Mat A, B;
3663: Vec a, b;
3665: PetscFunctionBegin;
3666: PetscCall(MatShellGetContext(S, static_cast<void *>(&h)));
3667: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, Y->cmap->n, nullptr, &A));
3668: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3669: PetscCall(MatDenseGetColumnVecRead(Y, i, &b));
3670: PetscCall(MatDenseGetColumnVecWrite(A, i, &a));
3671: PetscCall(VecISCopy(a, h->is[1], SCATTER_FORWARD, b));
3672: PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a));
3673: PetscCall(MatDenseRestoreColumnVecRead(Y, i, &b));
3674: }
3675: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B));
3676: PetscCall(KSPMatSolveTranspose(h->ksp, A, B));
3677: PetscCall(MatDestroy(&A));
3678: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->A[0]->rmap->n, B->cmap->n, nullptr, &A));
3679: for (PetscInt i = 0; i < A->cmap->n; ++i) {
3680: PetscCall(MatDenseGetColumnVecRead(B, i, &b));
3681: PetscCall(MatDenseGetColumnVecWrite(A, i, &a));
3682: PetscCall(VecISCopy(b, h->is[0], SCATTER_REVERSE, a));
3683: PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a));
3684: PetscCall(MatDenseRestoreColumnVecRead(B, i, &b));
3685: }
3686: PetscCall(MatDestroy(&B));
3687: PetscCall(MatTransposeMatMult(h->A[0], A, MAT_REUSE_MATRIX, PETSC_CURRENT, &X));
3688: PetscCall(MatDestroy(&A));
3689: PetscFunctionReturn(PETSC_SUCCESS);
3690: }
3692: static PetscErrorCode MatDestroy_Harmonic(Mat A)
3693: {
3694: Harmonic h;
3696: PetscFunctionBegin;
3697: PetscCall(MatShellGetContext(A, static_cast<void *>(&h)));
3698: for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i));
3699: PetscCall(PetscFree(h->is));
3700: PetscCall(VecDestroy(&h->v));
3701: for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i));
3702: PetscCall(PetscFree(h->A));
3703: PetscCall(KSPDestroy(&h->ksp));
3704: PetscCall(PetscFree(h));
3705: PetscFunctionReturn(PETSC_SUCCESS);
3706: }