Actual source code: ex14.c
1: static char help[] = "Tests that PCApplyTranspose() is the transpose of PCApply() for PCMG.\n\n";
3: /*
4: Checks that PCApplyTranspose() is the transpose of PCApply() for PCMG, for every cycle type,
5: with distinct up and down smoothers, and with a restriction that is not the transpose of the
6: interpolation (-distinct_restriction). The operator is a nonsymmetric 1D upwind
7: advection-diffusion stencil, the hierarchy is built by hand, and the smoothers come from the
8: options database. Two checks are made: the bilinear identity y . (M x) == x . (M^T y) on
9: random vectors, and M and M^T assembled column by column and compared entry for entry.
11: -n must be of the form m * 2^(levels-1) + 1.
12: */
14: #include <petscksp.h>
16: /*
17: Writes the 1D upwind advection-diffusion stencil into A. The advection makes the sub and
18: super diagonals differ, i.e., A != A^T, and the shift keeps every row strictly diagonally
19: dominant.
20: */
21: static PetscErrorCode SetOperatorValues1D(Mat A, PetscInt n, PetscReal advection, PetscReal shift)
22: {
23: PetscInt row_start, row_end;
25: PetscFunctionBeginUser;
26: PetscCall(MatGetOwnershipRange(A, &row_start, &row_end));
27: for (PetscInt i = row_start; i < row_end; i++) {
28: PetscInt cols[3], n_cols = 0;
29: PetscScalar vals[3];
31: if (i > 0) {
32: cols[n_cols] = i - 1;
33: vals[n_cols] = -1.0 - advection;
34: n_cols++;
35: }
36: cols[n_cols] = i;
37: vals[n_cols] = 2.0 + advection + shift;
38: n_cols++;
39: if (i < n - 1) {
40: cols[n_cols] = i + 1;
41: vals[n_cols] = -1.0;
42: n_cols++;
43: }
44: PetscCall(MatSetValues(A, 1, &i, n_cols, cols, vals, INSERT_VALUES));
45: }
46: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
47: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode BuildOperator(PetscInt n, PetscReal advection, PetscReal shift, Mat *A_out)
52: {
53: Mat A;
55: PetscFunctionBeginUser;
56: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
57: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
58: PetscCall(MatSetFromOptions(A));
59: PetscCall(MatSeqAIJSetPreallocation(A, 3, NULL));
60: PetscCall(MatMPIAIJSetPreallocation(A, 3, NULL, 2, NULL));
61: PetscCall(MatSetUp(A));
62: PetscCall(SetOperatorValues1D(A, n, advection, shift));
63: *A_out = A;
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*
68: Standard 1D linear interpolation from n_coarse points to n_fine = 2 * n_coarse - 1 points -
69: fine point 2j takes coarse point j, fine point 2j+1 the average of coarse points j and j+1.
70: */
71: static PetscErrorCode BuildInterpolation(PetscInt n_fine, PetscInt n_coarse, Mat *P_out)
72: {
73: Mat P;
74: PetscInt row_start, row_end;
76: PetscFunctionBeginUser;
77: PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
78: PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n_fine, n_coarse));
79: PetscCall(MatSetType(P, MATAIJ));
80: PetscCall(MatSeqAIJSetPreallocation(P, 2, NULL));
81: PetscCall(MatMPIAIJSetPreallocation(P, 2, NULL, 2, NULL));
82: PetscCall(MatGetOwnershipRange(P, &row_start, &row_end));
84: for (PetscInt i = row_start; i < row_end; i++) {
85: PetscInt cols[2], n_cols;
86: PetscScalar vals[2];
88: if (i % 2 == 0) {
89: n_cols = 1;
90: cols[0] = i / 2;
91: vals[0] = 1.0;
92: } else {
93: n_cols = 2;
94: cols[0] = (i - 1) / 2;
95: cols[1] = (i + 1) / 2;
96: vals[0] = 0.5;
97: vals[1] = 0.5;
98: }
99: PetscCall(MatSetValues(P, 1, &i, n_cols, cols, vals, INSERT_VALUES));
100: }
101: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
102: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
104: *P_out = P;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: /*
109: A restriction that is deliberately not the transpose of the interpolation, R = D P^T with D
110: a non-constant diagonal. A transposed cycle that restricted with R^T where it should have
111: used P, or interpolated with P where it should have used R^T, is invisible when R = P^T.
112: */
113: static PetscErrorCode BuildRestriction(Mat P, Mat *R_out)
114: {
115: Mat R;
116: Vec d;
117: PetscInt row_start, row_end, n_coarse;
119: PetscFunctionBeginUser;
120: PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));
121: PetscCall(MatGetSize(R, &n_coarse, NULL));
122: /* The left Vec of R, so the diagonal matches the row layout MatDiagonalScale() scales */
123: PetscCall(MatCreateVecs(R, NULL, &d));
124: PetscCall(VecGetOwnershipRange(d, &row_start, &row_end));
125: for (PetscInt j = row_start; j < row_end; j++) {
126: PetscScalar value = 1.0 + 0.5 * ((PetscReal)j / (PetscReal)n_coarse);
128: PetscCall(VecSetValue(d, j, value, INSERT_VALUES));
129: }
130: PetscCall(VecAssemblyBegin(d));
131: PetscCall(VecAssemblyEnd(d));
132: PetscCall(MatDiagonalScale(R, d, NULL));
133: PetscCall(VecDestroy(&d));
135: *R_out = R;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*
140: Builds the multigrid hierarchy by hand, with no DM, so that the interpolation and the
141: restriction are exactly what this test wants them to be. The coarse operators are formed by
142: PCMG itself as R A P.
143: */
144: static PetscErrorCode BuildHierarchy(PC pc, PetscInt levels, PetscInt n, PetscBool distinct_restriction)
145: {
146: PetscInt *level_size;
148: PetscFunctionBeginUser;
149: PetscCall(PetscMalloc1(levels, &level_size));
150: level_size[levels - 1] = n;
151: for (PetscInt l = levels - 2; l >= 0; l--) level_size[l] = (level_size[l + 1] + 1) / 2;
153: PetscCall(PCMGSetLevels(pc, levels, NULL));
154: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
156: for (PetscInt l = 1; l < levels; l++) {
157: Mat P;
159: PetscCall(BuildInterpolation(level_size[l], level_size[l - 1], &P));
160: PetscCall(PCMGSetInterpolation(pc, l, P));
161: /* With no restriction set PCMG uses P^T */
162: if (distinct_restriction) {
163: Mat R;
165: PetscCall(BuildRestriction(P, &R));
166: PetscCall(PCMGSetRestriction(pc, l, R));
167: PetscCall(MatDestroy(&R));
168: }
169: PetscCall(MatDestroy(&P));
170: }
172: PetscCall(PetscFree(level_size));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*
177: The transpose identity itself. Uses VecTDot() rather than VecDot() so that this is the
178: bilinear form in both real and complex builds - PCApplyTranspose() is the true transpose,
179: not the Hermitian transpose.
180: */
181: static PetscErrorCode CheckTransposeIdentity(PC pc, Mat A, PetscRandom rand, PetscInt n_pairs, PetscReal tol)
182: {
183: Vec x, y, mx, mty;
185: PetscFunctionBeginUser;
186: PetscCall(MatCreateVecs(A, &x, &mx));
187: PetscCall(MatCreateVecs(A, &y, &mty));
189: for (PetscInt i = 0; i < n_pairs; i++) {
190: PetscScalar lhs, rhs;
191: PetscReal diff, denom;
193: /* Random vectors - constant ones would hide a row/column mixup */
194: PetscCall(VecSetRandom(x, rand));
195: PetscCall(VecSetRandom(y, rand));
197: PetscCall(PCApply(pc, x, mx));
198: PetscCall(PCApplyTranspose(pc, y, mty));
200: PetscCall(VecTDot(mx, y, &lhs));
201: PetscCall(VecTDot(x, mty, &rhs));
203: diff = PetscAbsScalar(lhs - rhs);
204: denom = PetscMax(PetscAbsScalar(lhs), PetscAbsScalar(rhs));
205: if (denom < 1.0) denom = 1.0;
207: PetscCheck(diff / denom <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "bilinear identity check: y.(Mx) != x.(M^T y) on pair %" PetscInt_FMT ", relative difference %g > %g", i, (double)(diff / denom), (double)tol);
208: }
210: PetscCall(VecDestroy(&x));
211: PetscCall(VecDestroy(&y));
212: PetscCall(VecDestroy(&mx));
213: PetscCall(VecDestroy(&mty));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /*
218: Builds the preconditioner out explicitly, one column at a time, either as M or as M^T.
219: */
220: static PetscErrorCode BuildExplicit(PC pc, Mat A, PetscInt n, PetscBool transpose, Mat *out)
221: {
222: Mat dense;
223: Vec e;
224: PetscInt local_rows;
226: PetscFunctionBeginUser;
227: PetscCall(MatGetLocalSize(A, &local_rows, NULL));
228: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), local_rows, PETSC_DECIDE, n, n, NULL, &dense));
229: PetscCall(MatCreateVecs(A, NULL, &e));
231: for (PetscInt j = 0; j < n; j++) {
232: Vec col;
234: PetscCall(VecZeroEntries(e));
235: PetscCall(VecSetValue(e, j, 1.0, INSERT_VALUES));
236: PetscCall(VecAssemblyBegin(e));
237: PetscCall(VecAssemblyEnd(e));
239: PetscCall(MatDenseGetColumnVecWrite(dense, j, &col));
240: if (transpose) PetscCall(PCApplyTranspose(pc, e, col));
241: else PetscCall(PCApply(pc, e, col));
242: PetscCall(MatDenseRestoreColumnVecWrite(dense, j, &col));
243: }
244: PetscCall(VecDestroy(&e));
246: *out = dense;
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*
251: The sharp version of the check - build M and M^T out column by column and compare every
252: entry, rather than the single number the bilinear identity gives.
253: */
254: static PetscErrorCode CheckTransposeExplicitly(PC pc, Mat A, PetscInt n, PetscReal tol)
255: {
256: Mat m, mt, mt_transposed;
257: PetscReal diff, scale;
259: PetscFunctionBeginUser;
260: PetscCall(BuildExplicit(pc, A, n, PETSC_FALSE, &m));
261: PetscCall(BuildExplicit(pc, A, n, PETSC_TRUE, &mt));
263: /* (M^T)^T has to be M, entry for entry */
264: PetscCall(MatTranspose(mt, MAT_INITIAL_MATRIX, &mt_transposed));
265: PetscCall(MatNorm(m, NORM_FROBENIUS, &scale));
266: PetscCall(MatAXPY(mt_transposed, -1.0, m, SAME_NONZERO_PATTERN));
267: PetscCall(MatNorm(mt_transposed, NORM_FROBENIUS, &diff));
269: if (scale < 1.0) scale = 1.0;
270: PetscCheck(diff / scale <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "explicit transpose check: the explicit PCApplyTranspose() matrix is not the transpose of the explicit PCApply() matrix, relative Frobenius difference %g > %g", (double)(diff / scale), (double)tol);
272: PetscCall(MatDestroy(&m));
273: PetscCall(MatDestroy(&mt));
274: PetscCall(MatDestroy(&mt_transposed));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: int main(int argc, char **args)
279: {
280: Mat A;
281: PC pc;
282: PetscRandom rand;
283: PetscInt n = 33, levels = 3, n_pairs = 5, coarsening = 1;
284: PetscReal advection = 1.0, shift = 0.1, tol = PETSC_SMALL;
285: PetscBool distinct_restriction = PETSC_FALSE;
287: PetscFunctionBeginUser;
288: PetscCall(PetscInitialize(&argc, &args, NULL, help));
289: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
290: PetscCall(PetscOptionsGetInt(NULL, NULL, "-levels", &levels, NULL));
291: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_pairs", &n_pairs, NULL));
292: PetscCall(PetscOptionsGetReal(NULL, NULL, "-advection", &advection, NULL));
293: PetscCall(PetscOptionsGetReal(NULL, NULL, "-shift", &shift, NULL));
294: PetscCall(PetscOptionsGetReal(NULL, NULL, "-check_tol", &tol, NULL));
295: PetscCall(PetscOptionsGetBool(NULL, NULL, "-distinct_restriction", &distinct_restriction, NULL));
297: PetscCheck(levels > 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "-levels must be at least 2, not %" PetscInt_FMT, levels);
298: for (PetscInt l = 0; l < levels - 1; l++) coarsening *= 2;
299: PetscCheck(n > 2 * coarsening && (n - 1) % coarsening == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "-n %" PetscInt_FMT " must be larger than %" PetscInt_FMT " and of the form m * 2^(levels-1) + 1 so every level is a valid coarsening of the one above it", n, 2 * coarsening);
301: PetscCall(BuildOperator(n, advection, shift, &A));
303: /* Seeded so a failure is reproducible */
304: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
305: PetscCall(PetscRandomSetFromOptions(rand));
306: PetscCall(PetscRandomSetSeed(rand, 314159));
307: PetscCall(PetscRandomSeed(rand));
309: /* The cycle type and the smoothers all come from the options database */
310: PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
311: PetscCall(PCSetType(pc, PCMG));
312: PetscCall(PCSetOperators(pc, A, A));
313: PetscCall(BuildHierarchy(pc, levels, n, distinct_restriction));
314: PetscCall(PCSetFromOptions(pc));
315: PetscCall(PCSetUp(pc));
316: /* The PC is not driven by a KSP, so -pc_view has to be asked for explicitly */
317: PetscCall(PCViewFromOptions(pc, NULL, "-pc_view"));
319: PetscCall(CheckTransposeIdentity(pc, A, rand, n_pairs, tol));
320: PetscCall(CheckTransposeExplicitly(pc, A, n, tol));
321: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "PCApplyTranspose() is the transpose of PCApply()\n"));
323: PetscCall(PCDestroy(&pc));
324: PetscCall(PetscRandomDestroy(&rand));
325: PetscCall(MatDestroy(&A));
326: PetscCall(PetscFinalize());
327: return 0;
328: }
330: /*TEST
332: testset:
333: args: -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -mg_levels_ksp_max_it 2 -mg_coarse_ksp_type richardson -mg_coarse_pc_type jacobi -mg_coarse_ksp_max_it 4 -mg_coarse_ksp_norm_type none
334: output_file: output/ex14_1.out
335: nsize: {{1 2}}
337: test:
338: suffix: cycles
339: args: -pc_mg_type {{additive multiplicative full kaskade}shared output} -distinct_restriction {{0 1}shared output}
341: test:
342: suffix: cycles_distinct_smoothup
343: args: -pc_mg_type {{additive multiplicative full kaskade}shared output} -distinct_restriction {{0 1}shared output} -pc_mg_distinct_smoothup -mg_levels_up_ksp_type richardson -mg_levels_up_pc_type jacobi -mg_levels_up_ksp_max_it 1 -mg_levels_up_ksp_richardson_scale 0.7
345: test:
346: suffix: w_cycle
347: args: -pc_mg_type multiplicative -pc_mg_cycle_type w -distinct_restriction {{0 1}shared output}
349: test:
350: suffix: w_cycle_distinct_smoothup
351: args: -pc_mg_type multiplicative -pc_mg_cycle_type w -distinct_restriction {{0 1}shared output} -pc_mg_distinct_smoothup -mg_levels_up_ksp_type richardson -mg_levels_up_pc_type jacobi -mg_levels_up_ksp_max_it 1 -mg_levels_up_ksp_richardson_scale 0.7
353: test:
354: suffix: multiplicative_cycles
355: args: -pc_mg_type multiplicative -pc_mg_multiplicative_cycles 2 -distinct_restriction {{0 1}shared output}
357: test:
358: suffix: multiplicative_cycles_distinct_smoothup
359: args: -pc_mg_type multiplicative -pc_mg_multiplicative_cycles 2 -distinct_restriction {{0 1}shared output} -pc_mg_distinct_smoothup -mg_levels_up_ksp_type richardson -mg_levels_up_pc_type jacobi -mg_levels_up_ksp_max_it 1 -mg_levels_up_ksp_richardson_scale 0.7
361: test:
362: suffix: lu_coarse
363: nsize: 1
364: output_file: output/ex14_1.out
365: args: -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -mg_levels_ksp_max_it 2 -mg_coarse_ksp_type preonly -mg_coarse_pc_type lu -pc_mg_type {{additive multiplicative full kaskade}shared output} -distinct_restriction {{0 1}shared output}
367: TEST*/