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