Actual source code: ex16.c
1: static const char help[] = "Tests PCHPDDM with a user-supplied deflation space and scaling against an explicitly assembled reference operator.\n\n";
3: /*
4: Setup, on two processes:
5: - A is an 8 x 8 tridiagonal matrix, complex in complex builds, with 4 rows per process.
6: - The overlapping subdomain of process p holds the global rows 2p, ..., 2p + 5, see SubdomainRow().
7: - The scaling D_p is given by Scaling() and the single deflation vector U_p by Mode().
8: With -empty_rank0, process 0 supplies no deflation vector.
9: - ModifySubmatrices() shifts the local matrices of the fine-level PCASM, as optimized transmission
10: conditions would, so that the local solvers A_p differ from the restrictions of A.
12: After each setup, PCApply() and PCMatApply() must match the reference operator
13: T = Q + B (I - A Q) with -pc_hpddm_coarse_correction deflated,
14: T = B + Q (I - A B) with -pc_hpddm_coarse_correction deflated_reversed,
15: T = Q + (I - Q A) B (I - A Q) with -pc_hpddm_coarse_correction balanced,
16: where Q = Z (Z^H A Z)^{-1} Z^H is the coarse correction with Z = [R_p^T D_p U_p], and B is the
17: fine-level correction, see FineLevelType.
18: Also checks without deflation (-scaling_only) and then T = B.
19: */
21: #include <petscksp.h>
23: #define N_GLOBAL 8 // number of unknowns
24: #define N_LOCAL 4 // rows owned by each process
25: #define N_SUBDOMAIN 6 // rows of each overlapping subdomain
27: // fine-level corrections that the reference operator can reproduce
28: typedef enum {
29: FINE_ASM_WEIGHTED, // default with a user-supplied scaling: B = sum_p R_p^T D_p A_p^{-1} R_p
30: FINE_ASM_BASIC, // -pc_hpddm_levels_1_pc_asm_type basic: B = sum_p R_p^T A_p^{-1} R_p
31: FINE_JACOBI // -pc_hpddm_levels_1_pc_type jacobi: B = diag(A)^{-1}
32: } FineLevelType;
34: static PetscScalar TestScalar(PetscReal real, PetscReal imag)
35: {
36: #if PetscDefined(USE_COMPLEX)
37: return PetscCMPLX(real, imag);
38: #else
39: (void)imag;
40: return real;
41: #endif
42: }
44: // global row of the i-th row of the overlapping subdomain of process p: rows 0-5 for process 0, rows 2-7 for process 1
45: static PetscInt SubdomainRow(PetscMPIInt p, PetscInt i)
46: {
47: return 2 * p + i;
48: }
50: /*
51: Scaling of process p at a global row. Process 0 scales its rows 0-5 by 1, 1, 1, 0.65, 0.35, 0, and process 1 its
52: rows 2-7 by the complements 0, 0.35, 0.65, 1, 1, 1. The two scalings sum to one on every row and vanish outside each
53: subdomain. They also vanish on the subdomain row that A couples to the outside (row 5 for process 0, row 2 for
54: process 1), as PCHPDDMSetScaling() requires.
55: */
56: static PetscScalar Scaling(PetscMPIInt p, PetscInt row)
57: {
58: const PetscReal process0[N_GLOBAL] = {1.0, 1.0, 1.0, 0.65, 0.35, 0.0, 0.0, 0.0};
60: return p == 0 ? process0[row] : 1.0 - process0[row];
61: }
63: // entry of the deflation vector of process p at a global row; complex in complex builds, which exercises conjugation
64: static PetscScalar Mode(PetscMPIInt p, PetscInt row)
65: {
66: return TestScalar(1.0 + 0.15 * row, 0.3 * (p + 1) - 0.05 * row);
67: }
69: // shift applied to every fine-level local matrix by ModifySubmatrices(), and by the reference operator
70: static PetscScalar LocalShift(void)
71: {
72: return TestScalar(0.35, 0.2);
73: }
75: static PetscErrorCode ModifySubmatrices(PC pc, PetscInt n, const IS rows[], const IS cols[], Mat sub[], PetscCtx ctx)
76: {
77: PetscFunctionBeginUser;
78: for (PetscInt i = 0; i < n; ++i) PetscCall(MatShift(sub[i], LocalShift()));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode CreateOperator(Mat *A)
83: {
84: PetscInt start, end;
86: PetscFunctionBeginUser;
87: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, N_LOCAL, N_LOCAL, N_GLOBAL, N_GLOBAL, 3, NULL, 2, NULL, A));
88: PetscCall(MatGetOwnershipRange(*A, &start, &end));
89: for (PetscInt row = start; row < end; ++row) {
90: PetscCall(MatSetValue(*A, row, row, TestScalar(3.0 + 0.1 * row, 0.25), INSERT_VALUES));
91: if (row > 0) PetscCall(MatSetValue(*A, row, row - 1, TestScalar(-0.9, 0.07), INSERT_VALUES));
92: if (row < N_GLOBAL - 1) PetscCall(MatSetValue(*A, row, row + 1, TestScalar(-0.3, -0.12), INSERT_VALUES));
93: }
94: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
95: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*
100: Inputs of PCHPDDMSetDeflationMat() and PCHPDDMSetScaling(), as a user would build them: the subdomain IS,
101: the local deflation matrix U (one column, or none on process 0 with -empty_rank0), and the local scaling D. All three
102: list the subdomain rows in reverse order, unlike the sorted order used internally by PCASM, so that the test also
103: checks how PCHPDDM reorders the scaling.
104: */
105: static PetscErrorCode CreateDeflation(PetscBool empty_rank0, IS *is, Mat *U, Vec *D)
106: {
107: PetscMPIInt rank;
108: PetscInt rows[N_SUBDOMAIN], nmodes;
109: PetscScalar *u, *d;
111: PetscFunctionBeginUser;
112: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
113: nmodes = (empty_rank0 == PETSC_TRUE && rank == 0) ? 0 : 1; // number of local deflation vectors
114: for (PetscInt i = 0; i < N_SUBDOMAIN; ++i) rows[i] = SubdomainRow(rank, N_SUBDOMAIN - 1 - i);
115: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N_SUBDOMAIN, rows, PETSC_COPY_VALUES, is));
116: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, N_SUBDOMAIN, nmodes, NULL, U));
117: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N_SUBDOMAIN, D));
118: PetscCall(MatDenseGetArray(*U, &u));
119: PetscCall(VecGetArray(*D, &d));
120: for (PetscInt i = 0; i < N_SUBDOMAIN; ++i) {
121: d[i] = Scaling(rank, rows[i]);
122: if (nmodes > 0) u[i] = Mode(rank, rows[i]);
123: }
124: PetscCall(VecRestoreArray(*D, &d));
125: PetscCall(MatDenseRestoreArray(*U, &u));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*
130: Fine-level correction B, see FineLevelType. For PCASM, each process adds its term R_p^T [D_p] A_p^{-1} R_p column by
131: column, with A_p the shifted restriction of A to its subdomain, factored with LU like the PCASM local solver.
132: */
133: static PetscErrorCode CreateFineReference(Mat A, FineLevelType fine, Mat *B)
134: {
135: Mat *sub;
136: KSP ksp;
137: PC lu;
138: Vec x, y;
139: IS is;
140: PetscMPIInt rank;
141: PetscInt rows[N_SUBDOMAIN], start, end;
142: const PetscScalar *values;
143: PetscScalar column[N_SUBDOMAIN];
145: PetscFunctionBeginUser;
146: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
147: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, N_LOCAL, N_LOCAL, N_GLOBAL, N_GLOBAL, N_LOCAL, NULL, N_LOCAL, NULL, B)); // rows may be dense
148: if (fine == FINE_JACOBI) {
149: PetscCall(MatCreateVecs(A, NULL, &x));
150: PetscCall(MatGetDiagonal(A, x));
151: PetscCall(MatGetOwnershipRange(A, &start, &end));
152: PetscCall(VecGetArrayRead(x, &values));
153: for (PetscInt row = start; row < end; ++row) PetscCall(MatSetValue(*B, row, row, (PetscScalar)1.0 / values[row - start], INSERT_VALUES));
154: PetscCall(VecRestoreArrayRead(x, &values));
155: PetscCall(VecDestroy(&x));
156: } else {
157: for (PetscInt i = 0; i < N_SUBDOMAIN; ++i) rows[i] = SubdomainRow(rank, i);
158: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N_SUBDOMAIN, rows, PETSC_COPY_VALUES, &is));
159: PetscCall(MatCreateSubMatrices(A, 1, &is, &is, MAT_INITIAL_MATRIX, &sub));
160: PetscCall(MatShift(sub[0], LocalShift()));
161: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
162: PetscCall(KSPSetType(ksp, KSPPREONLY));
163: PetscCall(KSPSetOperators(ksp, sub[0], sub[0]));
164: PetscCall(KSPGetPC(ksp, &lu));
165: PetscCall(PCSetType(lu, PCLU));
166: PetscCall(MatCreateVecs(sub[0], &x, &y));
167: for (PetscInt j = 0; j < N_SUBDOMAIN; ++j) { // column j of A_p^{-1}, scaled by D_p for PC_ASM_WEIGHTED
168: PetscCall(VecSet(x, 0.0));
169: PetscCall(VecSetValue(x, j, 1.0, INSERT_VALUES));
170: PetscCall(VecAssemblyBegin(x));
171: PetscCall(VecAssemblyEnd(x));
172: PetscCall(KSPSolve(ksp, x, y));
173: PetscCall(VecGetArrayRead(y, &values));
174: for (PetscInt i = 0; i < N_SUBDOMAIN; ++i) column[i] = (fine == FINE_ASM_WEIGHTED ? Scaling(rank, rows[i]) : (PetscScalar)1.0) * values[i];
175: PetscCall(MatSetValues(*B, N_SUBDOMAIN, rows, 1, rows + j, column, ADD_VALUES));
176: PetscCall(VecRestoreArrayRead(y, &values));
177: }
178: PetscCall(VecDestroy(&x));
179: PetscCall(VecDestroy(&y));
180: PetscCall(KSPDestroy(&ksp));
181: PetscCall(MatDestroySubMatrices(1, &sub));
182: PetscCall(ISDestroy(&is));
183: }
184: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
185: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*
190: Coarse correction Q = Z (Z^H A Z)^{-1} Z^H. Column p of Z is the scaled deflation vector R_p^T D_p U_p of process p,
191: which is zero outside subdomain p because Scaling() is. With -empty_rank0, Z only has the column of process 1.
192: */
193: static PetscErrorCode CreateCoarseReference(Mat A, PetscBool empty_rank0, Mat *Q)
194: {
195: Mat Z, Zh, E, EinvZh;
196: KSP ksp;
197: PC redundant;
198: PetscMPIInt rank, first_process = empty_rank0 == PETSC_TRUE ? 1 : 0; // first process with a deflation vector
199: PetscInt start, end;
201: PetscFunctionBeginUser;
202: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
203: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, N_LOCAL, rank < first_process ? 0 : 1, N_GLOBAL, 2 - first_process, 1, NULL, 1, NULL, &Z));
204: PetscCall(MatGetOwnershipRange(Z, &start, &end));
205: for (PetscInt row = start; row < end; ++row)
206: for (PetscMPIInt p = first_process; p < 2; ++p) PetscCall(MatSetValue(Z, row, p - first_process, Scaling(p, row) * Mode(p, row), INSERT_VALUES));
207: PetscCall(MatAssemblyBegin(Z, MAT_FINAL_ASSEMBLY));
208: PetscCall(MatAssemblyEnd(Z, MAT_FINAL_ASSEMBLY));
209: PetscCall(MatHermitianTranspose(Z, MAT_INITIAL_MATRIX, &Zh));
210: PetscCall(MatMatMatMult(Zh, A, Z, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &E));
211: PetscCall(MatConvert(Zh, MATDENSE, MAT_INPLACE_MATRIX, &Zh)); // KSPMatSolve() needs dense right-hand sides
212: PetscCall(MatDuplicate(Zh, MAT_DO_NOT_COPY_VALUES, &EinvZh));
213: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
214: PetscCall(KSPSetType(ksp, KSPPREONLY));
215: PetscCall(KSPSetOperators(ksp, E, E));
216: PetscCall(KSPGetPC(ksp, &redundant));
217: PetscCall(PCSetType(redundant, PCREDUNDANT)); // gathers E on every process and factors it with the sequential PETSc LU
218: PetscCall(KSPMatSolve(ksp, Zh, EinvZh));
219: PetscCall(MatMatMult(Z, EinvZh, MAT_INITIAL_MATRIX, PETSC_DETERMINE, Q));
220: PetscCall(KSPDestroy(&ksp));
221: PetscCall(MatDestroy(&Z));
222: PetscCall(MatDestroy(&Zh));
223: PetscCall(MatDestroy(&E));
224: PetscCall(MatDestroy(&EinvZh));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*
229: Reference two-level operator T = first + second (I - A first), see the description at the top of this file: the
230: deflated correction applies the coarse level first, the reversed one applies the fine level first. The balanced
231: correction also multiplies second (I - A first) on the left by (I - first A).
232: */
233: static PetscErrorCode CreateReference(Mat A, PetscBool empty_rank0, FineLevelType fine, PCHPDDMCoarseCorrectionType correction, Mat *T)
234: {
235: Mat B, Q, first, second, IA;
237: PetscFunctionBeginUser;
238: PetscCheck(correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED || correction == PC_HPDDM_COARSE_CORRECTION_BALANCED, PETSC_COMM_WORLD, PETSC_ERR_SUP, "No reference for PCHPDDMCoarseCorrectionType %s, use %s, %s, or %s", PCHPDDMCoarseCorrectionTypes[correction], PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_DEFLATED], PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED], PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_BALANCED]);
239: PetscCall(CreateFineReference(A, fine, &B));
240: PetscCall(CreateCoarseReference(A, empty_rank0, &Q));
241: if (correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED_REVERSED) {
242: first = B;
243: second = Q;
244: } else {
245: first = Q;
246: second = B;
247: }
248: PetscCall(MatMatMult(A, first, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &IA));
249: PetscCall(MatScale(IA, -1.0));
250: PetscCall(MatShift(IA, 1.0)); // I - A first
251: PetscCall(MatMatMult(second, IA, MAT_INITIAL_MATRIX, PETSC_DETERMINE, T));
252: if (correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) {
253: Mat product;
255: PetscCall(MatDestroy(&IA));
256: PetscCall(MatMatMult(first, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &IA));
257: PetscCall(MatScale(IA, -1.0));
258: PetscCall(MatShift(IA, 1.0)); // I - first A
259: PetscCall(MatMatMult(IA, *T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &product));
260: PetscCall(MatDestroy(T));
261: *T = product;
262: }
263: PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B)); // T is dense, and MatAXPY() needs B in the same format; in place, so first still refers to B
264: PetscCall(MatAXPY(*T, 1.0, first, SAME_NONZERO_PATTERN));
265: PetscCall(MatDestroy(&B));
266: PetscCall(MatDestroy(&Q));
267: PetscCall(MatDestroy(&IA));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*
272: Fine-level correction that PCHPDDM builds from the same options, read with the prefix of its fine level:
273: pc_hpddm_levels_1_, or pc_hpddm_coarse_ with a single level. Only the corrections of FineLevelType have a reference.
274: */
275: static PetscErrorCode GetFineLevelType(const char prefix[], FineLevelType *fine)
276: {
277: char type[256] = PCASM;
278: PCASMType asm_type = PC_ASM_WEIGHTED; // default selected by PCHPDDM with a user-supplied scaling
279: PetscBool asm_pc, jacobi;
281: PetscFunctionBeginUser;
282: PetscCall(PetscOptionsGetString(NULL, prefix, "-pc_type", type, sizeof(type), NULL));
283: PetscCall(PetscOptionsGetEnum(NULL, prefix, "-pc_asm_type", PCASMTypes, (PetscEnum *)&asm_type, NULL));
284: PetscCall(PetscStrcmp(type, PCASM, &asm_pc));
285: PetscCall(PetscStrcmp(type, PCJACOBI, &jacobi));
286: PetscCheck(jacobi == PETSC_TRUE || (asm_pc == PETSC_TRUE && (asm_type == PC_ASM_WEIGHTED || asm_type == PC_ASM_BASIC)), PETSC_COMM_WORLD, PETSC_ERR_SUP, "No reference for this fine-level PC, use PCJACOBI or PCASM of type PC_ASM_WEIGHTED or PC_ASM_BASIC");
287: if (jacobi == PETSC_TRUE) *fine = FINE_JACOBI;
288: else *fine = asm_type == PC_ASM_WEIGHTED ? FINE_ASM_WEIGHTED : FINE_ASM_BASIC;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: static PetscErrorCode CheckVector(Vec actual, Vec expected)
293: {
294: PetscReal error, norm;
296: PetscFunctionBeginUser;
297: PetscCall(VecNorm(expected, NORM_INFINITY, &norm));
298: PetscCall(VecAXPY(expected, -1.0, actual));
299: PetscCall(VecNorm(expected, NORM_INFINITY, &error));
300: PetscCheck(error <= 1000.0 * PETSC_MACHINE_EPSILON * PetscMax(1.0, norm), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "PCHPDDM differs from the reference operator by %g", (double)error);
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /*
305: Compares the PC with the reference operator T: PCApply() on every unit vector, and PCMatApply() on a block of three
306: vectors. Also compares PCApplyTranspose() and PCMatApplyTranspose() with T^T, including with complex scalars.
307: */
308: static PetscErrorCode CheckApplications(PC pc, Mat T)
309: {
310: Vec x, y, expected;
311: PetscInt N, m, start, end;
313: PetscFunctionBeginUser;
314: PetscCall(MatGetSize(T, &N, NULL));
315: PetscCall(MatGetLocalSize(T, &m, NULL));
316: PetscCall(MatCreateVecs(T, &x, &y));
317: PetscCall(VecDuplicate(y, &expected));
318: PetscCall(VecGetOwnershipRange(x, &start, &end));
319: for (PetscInt j = 0; j < N; ++j) {
320: PetscCall(VecSet(x, 0.0));
321: if (start <= j && j < end) PetscCall(VecSetValue(x, j, 1.0, INSERT_VALUES));
322: PetscCall(VecAssemblyBegin(x));
323: PetscCall(VecAssemblyEnd(x));
324: PetscCall(PCApply(pc, x, y));
325: PetscCall(MatMult(T, x, expected));
326: PetscCall(CheckVector(y, expected));
327: PetscCall(PCApplyTranspose(pc, x, y));
328: PetscCall(MatMultTranspose(T, x, expected));
329: PetscCall(CheckVector(y, expected));
330: }
331: {
332: Mat X, Y;
334: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, N, 3, NULL, &X));
335: PetscCall(MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &Y));
336: PetscCall(MatSetRandom(X, NULL));
337: for (PetscInt transpose = 0; transpose < 2; ++transpose) {
338: if (transpose) PetscCall(PCMatApplyTranspose(pc, X, Y));
339: else PetscCall(PCMatApply(pc, X, Y));
340: for (PetscInt j = 0; j < 3; ++j) {
341: Vec xc, yc;
343: PetscCall(MatDenseGetColumnVecRead(X, j, &xc));
344: PetscCall(MatDenseGetColumnVecRead(Y, j, &yc));
345: if (transpose) PetscCall(MatMultTranspose(T, xc, expected));
346: else PetscCall(MatMult(T, xc, expected));
347: PetscCall(CheckVector(yc, expected));
348: PetscCall(MatDenseRestoreColumnVecRead(Y, j, &yc));
349: PetscCall(MatDenseRestoreColumnVecRead(X, j, &xc));
350: }
351: }
352: PetscCall(MatDestroy(&X));
353: PetscCall(MatDestroy(&Y));
354: }
355: PetscCall(VecDestroy(&x));
356: PetscCall(VecDestroy(&y));
357: PetscCall(VecDestroy(&expected));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: // checks that the PC holds the shared IS of PCHPDDMSetDeflationMat() and PCHPDDMSetScaling() exactly when one of them uses it
362: static PetscErrorCode CheckReferences(IS is, PetscBool held)
363: {
364: PetscInt refct;
366: PetscFunctionBeginUser;
367: PetscCall(PetscObjectGetReference((PetscObject)is, &refct));
368: PetscCheck(refct == (held ? 2 : 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "IS has %" PetscInt_FMT " references, expected %d", refct, held ? 2 : 1);
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: int main(int argc, char **argv)
373: {
374: PC pc;
375: Mat A, U, T;
376: IS is;
377: Vec D;
378: PetscMPIInt size;
379: PetscBool empty_rank0 = PETSC_FALSE, scaling_only = PETSC_FALSE;
380: FineLevelType fine = FINE_ASM_WEIGHTED;
381: PCHPDDMCoarseCorrectionType correction;
383: PetscFunctionBeginUser;
384: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
385: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
386: PetscCheck(size == 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This example requires 2 processes");
387: PetscCall(PetscOptionsGetBool(NULL, NULL, "-empty_rank0", &empty_rank0, NULL));
388: PetscCall(PetscOptionsGetBool(NULL, NULL, "-scaling_only", &scaling_only, NULL)); // no deflation space, the scaling only applies to the fine level
389: PetscCall(GetFineLevelType(scaling_only ? "pc_hpddm_coarse_" : "pc_hpddm_levels_1_", &fine));
390: PetscCall(CreateOperator(&A));
391: PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
392: PetscCall(PCSetType(pc, PCHPDDM));
393: /*
394: Three setups of the same PC, each checked against the reference operator:
395: 0. initial setup, after installing and removing the deflation space and scaling in both orders, which must leave
396: no trace, and with the PC holding their shared IS exactly while one of them uses it;
397: 1. new operator values with the same deflation space and scaling: both levels are rebuilt, and the fine-level
398: PCASMType selected at the initial setup must persist;
399: 2. PCReset() discards the deflation space and scaling, which are then installed again, the scaling first.
400: */
401: for (PetscInt pass = 0; pass < 3; ++pass) {
402: if (pass == 1) PetscCall(MatScale(A, 1.1));
403: else {
404: if (pass == 2) PetscCall(PCReset(pc));
405: PetscCall(CreateDeflation(empty_rank0, &is, &U, &D));
406: if (scaling_only) PetscCall(PCHPDDMSetScaling(pc, is, D));
407: else if (pass == 0) {
408: PetscCall(PCHPDDMSetScaling(pc, is, D)); // the scaling alone stores the IS
409: PetscCall(CheckReferences(is, PETSC_TRUE));
410: PetscCall(PCHPDDMSetDeflationMat(pc, NULL, U)); // reuses the stored IS and keeps the scaling
411: PetscCall(PCHPDDMSetScaling(pc, NULL, NULL));
412: PetscCall(CheckReferences(is, PETSC_TRUE)); // still used by the deflation space
413: PetscCall(PCHPDDMSetScaling(pc, NULL, D));
414: PetscCall(PCHPDDMSetDeflationMat(pc, NULL, NULL));
415: PetscCall(CheckReferences(is, PETSC_TRUE)); // still used by the scaling
416: PetscCall(PCHPDDMSetDeflationMat(pc, is, U));
417: PetscCall(PCHPDDMSetScaling(pc, NULL, NULL));
418: PetscCall(PCHPDDMSetDeflationMat(pc, NULL, NULL));
419: PetscCall(CheckReferences(is, PETSC_FALSE));
420: }
421: if (!scaling_only) {
422: if (pass == 2) { // the deflation space reuses the IS stored by the scaling
423: PetscCall(PCHPDDMSetScaling(pc, is, D));
424: PetscCall(PCHPDDMSetDeflationMat(pc, NULL, U));
425: } else {
426: IS copy;
428: PetscCall(PCHPDDMSetDeflationMat(pc, is, U));
429: PetscCall(ISDuplicate(is, ©)); // replace the stored IS while keeping the deflation matrix
430: PetscCall(PCHPDDMSetScaling(pc, copy, D));
431: PetscCall(CheckReferences(is, PETSC_FALSE));
432: PetscCall(CheckReferences(copy, PETSC_TRUE));
433: PetscCall(ISDestroy(©));
434: }
435: }
436: PetscCall(VecDestroy(&D)); // the PC keeps its own reference
437: PetscCall(ISDestroy(&is));
438: PetscCall(MatDestroy(&U));
439: }
440: PetscCall(PCSetOperators(pc, A, A));
441: PetscCall(PCSetModifySubMatrices(pc, ModifySubmatrices, NULL));
442: PetscCall(PCSetFromOptions(pc));
443: PetscCall(PCSetUp(pc));
444: if (scaling_only) PetscCall(CreateFineReference(A, fine, &T));
445: else {
446: PetscCall(PCHPDDMGetCoarseCorrectionType(pc, &correction));
447: PetscCall(CreateReference(A, empty_rank0, fine, correction, &T));
448: }
449: PetscCall(CheckApplications(pc, T));
450: PetscCall(MatDestroy(&T));
451: }
452: PetscCall(PCDestroy(&pc));
453: PetscCall(MatDestroy(&A));
454: PetscCall(PetscFinalize());
455: return 0;
456: }
458: /*TEST
460: testset:
461: # every case runs with and without a deflation vector on process 0
462: nsize: 2
463: requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
464: args: -pc_hpddm_define_subdomains true -pc_hpddm_coarse_mat_type aij -empty_rank0 {{false true}}
465: output_file: output/empty.out
466: test:
467: # default fine level: PCASM of type PC_ASM_WEIGHTED, with the same scaling as the coarse space
468: suffix: weighted
469: args: -pc_hpddm_levels_1_sub_pc_type lu
470: test:
471: # same fine level with the fine-level correction applied first
472: suffix: weighted_reversed
473: args: -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_coarse_correction deflated_reversed
474: test:
475: suffix: balanced
476: args: -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_coarse_correction balanced
477: test:
478: # PCASM of type PC_ASM_BASIC: the scaling only applies to the coarse space
479: # with one subdomain per process, multiplicative composition gives the same operator as additive, but it checks
480: # that the second setup keeps the selected PCASMType instead of reverting to PC_ASM_WEIGHTED
481: suffix: basic
482: args: -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_pc_asm_type basic -pc_hpddm_levels_1_pc_asm_local_type {{additive multiplicative}}
483: test:
484: # non-PCASM fine level: the scaling only applies to the coarse space
485: suffix: jacobi
486: args: -pc_hpddm_levels_1_pc_type jacobi
488: test:
489: # no deflation space: the scaling only applies to the fine-level PCASM of a single-level PC, whose prefix is -pc_hpddm_coarse_
490: suffix: scaling_only
491: nsize: 2
492: requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
493: args: -pc_hpddm_define_subdomains true -scaling_only -pc_hpddm_coarse_sub_pc_type lu
494: output_file: output/empty.out
496: TEST*/