Actual source code: ex18.c
1: static char help[] = "Tests the use of MatZeroRowsColumns() for parallel matrices.\n\
2: Contributed-by: Stephan Kramer <s.kramer@imperial.ac.uk>\n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A;
9: Vec x, rhs, y;
10: PetscInt i, j, k, b, m = 3, n, nlocal = 2, bs = 1, Ii, J;
11: PetscInt *boundary_nodes, nboundary_nodes, *boundary_indices;
12: PetscMPIInt rank, size;
13: PetscScalar v, v0, v1, v2, a0 = 0.1, a, rhsval, *boundary_values, diag = 1.0;
14: PetscReal norm;
15: char convname[64];
16: PetscBool upwind = PETSC_FALSE, nonlocalBC = PETSC_FALSE, zerorhs = PETSC_TRUE, convert = PETSC_FALSE, missingdiag = PETSC_FALSE;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, NULL, help));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
22: n = nlocal * size;
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
25: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonlocal_bc", &nonlocalBC, NULL));
26: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-diag", &diag, NULL));
27: PetscCall(PetscOptionsGetString(NULL, NULL, "-convname", convname, sizeof(convname), &convert));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-zerorhs", &zerorhs, NULL));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-missing_diagonal", &missingdiag, NULL));
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n * bs, m * n * bs));
33: PetscCall(MatSetFromOptions(A));
34: PetscCall(MatSetUp(A));
36: PetscCall(MatCreateVecs(A, NULL, &rhs));
37: PetscCall(VecSetFromOptions(rhs));
38: PetscCall(VecSetUp(rhs));
40: rhsval = 0.0;
41: for (i = 0; i < m; i++) {
42: for (j = nlocal * rank; j < nlocal * (rank + 1); j++) {
43: a = a0;
44: for (b = 0; b < bs; b++) {
45: /* let's start with a 5-point stencil diffusion term */
46: v = -1.0;
47: Ii = (j + n * i) * bs + b;
48: if (i > 0) {
49: J = Ii - n * bs;
50: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
51: }
52: if (i < m - 1) {
53: J = Ii + n * bs;
54: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
55: }
56: if (j > 0) {
57: J = Ii - 1 * bs;
58: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
59: }
60: if (j < n - 1) {
61: J = Ii + 1 * bs;
62: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
63: }
64: v = 4.0;
65: if (!missingdiag || Ii != m * n * bs - 1) PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
66: if (upwind) {
67: /* now add a 2nd order upwind advection term to add a little asymmetry */
68: if (j > 2) {
69: J = Ii - 2 * bs;
70: v2 = 0.5 * a;
71: v1 = -2.0 * a;
72: v0 = 1.5 * a;
73: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v2, ADD_VALUES));
74: } else {
75: /* fall back to 1st order upwind */
76: v1 = -1.0 * a;
77: v0 = 1.0 * a;
78: }
79: if (j > 1) {
80: J = Ii - 1 * bs;
81: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v1, ADD_VALUES));
82: }
83: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v0, ADD_VALUES));
84: a /= 10.; /* use a different velocity for the next component */
85: /* add a coupling to the previous and next components */
86: v = 0.5;
87: if (b > 0) {
88: J = Ii - 1;
89: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
90: }
91: if (b < bs - 1) {
92: J = Ii + 1;
93: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
94: }
95: }
96: /* make up some rhs */
97: PetscCall(VecSetValue(rhs, Ii, rhsval, INSERT_VALUES));
98: rhsval += 1.0;
99: }
100: }
101: }
102: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
103: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
105: if (convert) { /* Test different Mat implementations */
106: Mat B;
108: PetscCall(MatConvert(A, convname, MAT_INITIAL_MATRIX, &B));
109: PetscCall(MatDestroy(&A));
110: A = B;
111: }
113: /* the diagonal shift is only inserted where the diagonal entry is already allocated, so it must not need new nonzeros */
114: if (missingdiag) PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
116: PetscCall(VecAssemblyBegin(rhs));
117: PetscCall(VecAssemblyEnd(rhs));
118: /* set rhs to zero to simplify */
119: if (zerorhs) PetscCall(VecZeroEntries(rhs));
121: if (nonlocalBC) {
122: /*version where boundary conditions are set by processes that don't necessarily own the nodes */
123: if (rank == 0) {
124: nboundary_nodes = size > m ? nlocal : m - size + nlocal;
125: PetscCall(PetscMalloc1(nboundary_nodes, &boundary_nodes));
126: k = 0;
127: for (i = size; i < m; i++, k++) boundary_nodes[k] = n * i;
128: } else if (rank < m) {
129: nboundary_nodes = nlocal + 1;
130: PetscCall(PetscMalloc1(nboundary_nodes, &boundary_nodes));
131: boundary_nodes[0] = rank * n;
132: k = 1;
133: } else {
134: nboundary_nodes = nlocal;
135: PetscCall(PetscMalloc1(nboundary_nodes, &boundary_nodes));
136: k = 0;
137: }
138: for (j = nlocal * rank; j < nlocal * (rank + 1); j++, k++) boundary_nodes[k] = j;
139: } else {
140: /*version where boundary conditions are set by the node owners only */
141: PetscCall(PetscMalloc1(m * n, &boundary_nodes));
142: k = 0;
143: for (j = 0; j < n; j++) {
144: Ii = j;
145: if (Ii >= rank * m * nlocal && Ii < (rank + 1) * m * nlocal) boundary_nodes[k++] = Ii;
146: }
147: for (i = 1; i < m; i++) {
148: Ii = n * i;
149: if (Ii >= rank * m * nlocal && Ii < (rank + 1) * m * nlocal) boundary_nodes[k++] = Ii;
150: }
151: nboundary_nodes = k;
152: }
154: PetscCall(VecDuplicate(rhs, &x));
155: PetscCall(PetscMalloc2(nboundary_nodes * bs, &boundary_indices, nboundary_nodes * bs, &boundary_values));
156: for (k = 0; k < nboundary_nodes; k++) {
157: Ii = boundary_nodes[k] * bs;
158: v = 1.0 * boundary_nodes[k];
159: for (b = 0; b < bs; b++, Ii++) {
160: boundary_indices[k * bs + b] = Ii;
161: boundary_values[k * bs + b] = v;
162: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d %" PetscInt_FMT " %f\n", rank, Ii, (double)PetscRealPart(v)));
163: v += 0.1;
164: }
165: }
166: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
167: PetscCall(VecSetValues(x, nboundary_nodes * bs, boundary_indices, boundary_values, INSERT_VALUES));
168: PetscCall(VecAssemblyBegin(x));
169: PetscCall(VecAssemblyEnd(x));
171: /* We can check the rhs returned by MatZeroColumns by computing y=rhs-A*x and overwriting the boundary entries with boundary values */
172: PetscCall(VecDuplicate(x, &y));
173: PetscCall(MatMult(A, x, y));
174: PetscCall(VecAYPX(y, -1.0, rhs));
175: for (k = 0; k < nboundary_nodes * bs; k++) boundary_values[k] *= diag;
176: PetscCall(VecSetValues(y, nboundary_nodes * bs, boundary_indices, boundary_values, INSERT_VALUES));
177: PetscCall(VecAssemblyBegin(y));
178: PetscCall(VecAssemblyEnd(y));
180: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Matrix A and vector x:\n"));
181: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
182: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
184: PetscCall(MatZeroRowsColumns(A, nboundary_nodes * bs, boundary_indices, diag, x, rhs));
185: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Vector rhs returned by MatZeroRowsColumns\n"));
186: PetscCall(VecView(rhs, PETSC_VIEWER_STDOUT_WORLD));
187: PetscCall(VecAXPY(y, -1.0, rhs));
188: PetscCall(VecNorm(y, NORM_INFINITY, &norm));
189: if (norm > 1.0e-10) {
190: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*** Difference between rhs and y, inf-norm: %f\n", (double)norm));
191: PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
192: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Bug in MatZeroRowsColumns");
193: }
195: PetscCall(PetscFree(boundary_nodes));
196: PetscCall(PetscFree2(boundary_indices, boundary_values));
197: PetscCall(VecDestroy(&x));
198: PetscCall(VecDestroy(&y));
199: PetscCall(VecDestroy(&rhs));
200: PetscCall(MatDestroy(&A));
202: PetscCall(PetscFinalize());
203: return 0;
204: }
206: /*TEST
208: test:
209: diff_args: -j
210: suffix: 0
212: test:
213: diff_args: -j
214: suffix: 1
215: nsize: 2
217: test:
218: diff_args: -j
219: suffix: 10
220: nsize: 2
221: args: -bs 2 -nonlocal_bc
223: test:
224: diff_args: -j
225: suffix: 11
226: nsize: 7
227: args: -bs 2 -nonlocal_bc
229: test:
230: diff_args: -j
231: suffix: 12
232: args: -bs 2 -nonlocal_bc -mat_type baij
234: test:
235: diff_args: -j
236: suffix: 13
237: nsize: 2
238: args: -bs 2 -nonlocal_bc -mat_type baij
240: test:
241: diff_args: -j
242: suffix: 14
243: nsize: 7
244: args: -bs 2 -nonlocal_bc -mat_type baij
246: test:
247: diff_args: -j
248: suffix: 2
249: nsize: 7
251: test:
252: diff_args: -j
253: suffix: 3
254: args: -mat_type baij
256: test:
257: diff_args: -j
258: suffix: 4
259: nsize: 2
260: args: -mat_type baij
262: test:
263: diff_args: -j
264: suffix: 5
265: nsize: 7
266: args: -mat_type baij
268: test:
269: diff_args: -j
270: suffix: 6
271: args: -bs 2 -mat_type baij
273: test:
274: diff_args: -j
275: suffix: 7
276: nsize: 2
277: args: -bs 2 -mat_type baij
279: test:
280: diff_args: -j
281: suffix: 8
282: nsize: 7
283: args: -bs 2 -mat_type baij
285: test:
286: diff_args: -j
287: suffix: 9
288: args: -bs 2 -nonlocal_bc
290: test:
291: diff_args: -j
292: suffix: 15
293: args: -bs 2 -nonlocal_bc -convname shell
295: test:
296: diff_args: -j
297: suffix: 16
298: nsize: 2
299: args: -bs 2 -nonlocal_bc -convname shell
301: test:
302: diff_args: -j
303: suffix: 17
304: args: -bs 2 -nonlocal_bc -convname dense
306: testset:
307: diff_args: -j
308: suffix: full
309: nsize: {{1 3}separate output}
310: args: -diag {{0.12 -0.13}separate output} -convname {{aij shell baij}separate output} -zerorhs 0
312: test:
313: diff_args: -j
314: requires: cuda
315: suffix: cusparse_1
316: nsize: 1
317: args: -diag {{0.12 -0.13}separate output} -convname {{seqaijcusparse mpiaijcusparse}separate output} -zerorhs 0 -mat_type {{seqaijcusparse mpiaijcusparse}separate output}
319: test:
320: diff_args: -j
321: requires: cuda
322: suffix: cusparse_3
323: nsize: 3
324: args: -diag {{0.12 -0.13}separate output} -convname mpiaijcusparse -zerorhs 0 -mat_type mpiaijcusparse
326: test:
327: suffix: 18
328: nsize: 3
329: args: -diag 0.12 -convname {{mpisbaij is}separate output} -zerorhs 0
331: testset:
332: diff_args: -j
333: suffix: missing_diagonal
334: nsize: {{1 3}separate output}
335: args: -missing_diagonal -diag 0.12 -zerorhs 0
337: TEST*/