Actual source code: ex254.c
1: static char help[] = "Test MatSetValuesCOO() for MPIAIJ and its subclasses \n\n";
3: #include <petscmat.h>
5: /*
6: Checks that MatSetValuesCOO() can still be called after MatShift() on a matrix with m local rows and n local columns, whose COO
7: pattern already contains every diagonal entry so the shift creates no new nonzeros
8: */
9: static PetscErrorCode TestShiftCOO(PetscInt m, PetscInt n)
10: {
11: Mat D;
12: PetscInt rstart, rend, N, ncoo;
13: PetscInt *coo_i, *coo_j;
14: PetscScalar readback, expected;
15: PetscScalar *coo_v;
17: PetscFunctionBeginUser;
18: PetscCallMPI(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
19: PetscCall(MatCreate(PETSC_COMM_WORLD, &D));
20: PetscCall(MatSetSizes(D, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
21: PetscCall(MatSetFromOptions(D));
22: PetscCall(MatGetOwnershipRange(D, &rstart, &rend));
23: ncoo = 2 * (rend - rstart);
24: PetscCall(PetscMalloc3(ncoo, &coo_i, ncoo, &coo_j, ncoo, &coo_v));
25: for (PetscInt r = rstart, k = 0; r < rend; r++) {
26: coo_i[k] = r;
27: coo_j[k++] = r % N;
28: coo_i[k] = r;
29: coo_j[k++] = (r + 1) % N;
30: }
31: for (PetscInt k = 0; k < ncoo; k++) coo_v[k] = 1.0;
32: PetscCall(MatSetPreallocationCOO(D, ncoo, coo_i, coo_j));
33: PetscCall(MatSetValuesCOO(D, coo_v, INSERT_VALUES));
34: PetscCall(MatShift(D, 2.0));
35: PetscCall(MatShift(D, 2.0));
36: PetscCall(MatSetValuesCOO(D, coo_v, ADD_VALUES));
37: for (PetscInt r = rstart; r < rend; r++) {
38: PetscCall(MatGetValue(D, r, r % N, &readback));
39: expected = r < N ? 6.0 : 2.0;
40: if (readback != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSetValuesCOO() after MatShift() gave %g in row %" PetscInt_FMT ", expected %g\n", (double)PetscRealPart(readback), r, (double)PetscRealPart(expected)));
41: }
42: PetscCall(PetscFree3(coo_i, coo_j, coo_v));
43: PetscCall(MatDestroy(&D));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: int main(int argc, char **args)
48: {
49: Mat A, B, C;
50: const PetscInt M = 18, N = 18;
51: PetscBool equal, isHypre;
52: PetscScalar *vals;
53: PetscBool flg = PETSC_FALSE, freecoo = PETSC_FALSE, missing_diagonal = PETSC_FALSE;
54: PetscInt ncoos = 1;
55: PetscMPIInt size;
57: // clang-format off
58: /* Construct 18 x 18 matrices, which are big enough to have complex communication patterns but still small enough for debugging */
59: // i0/j0[] has a dense diagonal
60: PetscInt i0[] = {7, 7, 8, 8, 9, 16, 17, 9, 10, 1, 1, -2, 2, 3, 3, 14, 4, 5, 10, 13, 9, 9, 10, 1, 0, 0, 5, 5, 6, 6, 13, 13, 14, -14, 4, 4, 5, 11, 11, 12, 15, 15, 16};
61: PetscInt j0[] = {7, 6, 8, 4, 9, 16, 17, 16, 10, 2, 1, 3, 2, 4, 3, 14, 4, 5, 15, 13, 10, 16, 11, 2, 0, 1, 5, -11, 0, 6, 15, 17, 11, 13, 4, 8, 2, 11, 17, 12, 3, 15, 9};
63: // i0/j0[] miss some diagonals
64: PetscInt i1[] = {8, 5, 15, 16, 6, 13, 4, 17, 8, 9, 9, 10, -6, 12, 7, 3, -4, 1, 1, 2, 5, 5, 6, 14, 17, 8, 9, 9, 10, 4, 5, 10, 11, 1, 2};
65: PetscInt j1[] = {2, 3, 16, 9, 5, 17, 1, 13, 4, 10, 16, 11, -5, 12, 1, 7, -1, 2, 7, 3, 6, 11, 0, 11, 13, 4, 10, 16, 11, 8, -2, 15, 12, 7, 3};
67: PetscInt i2[] = {3, 4, 1, 10, 0, 1, 1, 2, 1, 1, 2, 2, 3, 3, 4, 4, 1, 2, 5, 5, 6, 4, 17, 0, 1, 1, 8, 5, 5, 6, 4, 7, 8, 5};
68: PetscInt j2[] = {7, 1, 2, 11, 5, 2, 7, 3, 2, 7, 3, 8, 4, 9, 3, 5, 7, 3, 6, 11, 0, 1, 13, 5, 2, 7, 4, 6, 11, 0, 1, 3, 4, 2};
69: // clang-format on
71: typedef struct {
72: PetscInt *i, *j, n;
73: } coo_data;
75: coo_data coos[3] = {
76: {i0, j0, PETSC_STATIC_ARRAY_LENGTH(i0)},
77: {i1, j1, PETSC_STATIC_ARRAY_LENGTH(i1)},
78: {i2, j2, PETSC_STATIC_ARRAY_LENGTH(i2)}
79: };
80: coo_data mycoo;
82: PetscFunctionBeginUser;
83: PetscCall(PetscInitialize(&argc, &args, NULL, help));
84: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
85: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ignore_remote", &flg, NULL));
86: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ncoos", &ncoos, NULL));
87: PetscCall(PetscOptionsGetBool(NULL, NULL, "-missing_diagonal", &missing_diagonal, NULL));
89: mycoo.n = 0;
90: if (ncoos > 1) {
91: PetscLayout map;
93: freecoo = PETSC_TRUE;
94: PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD, &map));
95: PetscCall(PetscLayoutSetSize(map, ncoos));
96: PetscCall(PetscLayoutSetUp(map));
97: PetscCall(PetscLayoutGetLocalSize(map, &ncoos));
98: for (PetscInt i = 0; i < ncoos; i++) mycoo.n += coos[i % 3].n;
99: PetscCall(PetscMalloc2(mycoo.n, &mycoo.i, mycoo.n, &mycoo.j));
100: mycoo.n = 0;
101: for (PetscInt i = 0; i < ncoos; i++) {
102: PetscCall(PetscArraycpy(mycoo.i + mycoo.n, coos[i % 3].i, coos[i % 3].n));
103: PetscCall(PetscArraycpy(mycoo.j + mycoo.n, coos[i % 3].j, coos[i % 3].n));
104: mycoo.n += coos[i % 3].n;
105: }
106: PetscCall(PetscLayoutDestroy(&map));
107: } else if (ncoos == 1 && PetscGlobalRank < 3) mycoo = coos[PetscGlobalRank];
109: if (missing_diagonal && PetscGlobalRank == 0) mycoo = coos[1];
111: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
112: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
113: PetscCall(MatSetType(A, MATAIJ));
114: // Do not preallocate A to also test MatHash with MAT_IGNORE_OFF_PROC_ENTRIES
115: // PetscCall(MatSeqAIJSetPreallocation(A, 2, NULL));
116: // PetscCall(MatMPIAIJSetPreallocation(A, 2, NULL, 2, NULL));
117: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
118: PetscCall(MatSetOption(A, MAT_IGNORE_OFF_PROC_ENTRIES, flg));
120: PetscCall(PetscMalloc1(mycoo.n, &vals));
121: for (PetscInt k = 0; k < mycoo.n; k++) {
122: vals[k] = mycoo.j[k];
123: PetscCall(MatSetValue(A, mycoo.i[k], mycoo.j[k], vals[k], ADD_VALUES));
124: }
125: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(MatViewFromOptions(A, NULL, "-a_view"));
129: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
130: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, N));
131: PetscCall(MatSetFromOptions(B));
132: PetscCall(MatSetOption(B, MAT_IGNORE_OFF_PROC_ENTRIES, flg));
133: PetscCall(MatSetPreallocationCOO(B, mycoo.n, mycoo.i, mycoo.j));
135: /* Test with ADD_VALUES on a zeroed matrix */
136: PetscCall(MatSetValuesCOO(B, vals, ADD_VALUES));
137: PetscCall(MatMultEqual(A, B, 10, &equal));
138: if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSetValuesCOO() failed\n"));
139: PetscCall(MatViewFromOptions(B, NULL, "-b_view"));
141: /* Test with MatDuplicate on a zeroed matrix */
142: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &C));
143: PetscCall(MatSetValuesCOO(C, vals, ADD_VALUES));
144: PetscCall(MatMultEqual(A, C, 10, &equal));
145: if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSetValuesCOO() on duplicated matrix failed\n"));
146: PetscCall(MatViewFromOptions(C, NULL, "-c_view"));
148: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &isHypre));
150: /* Test MatSetValues() on a matrix preallocated with MatSetPreallocationCOO() with no intervening
151: MatSetValuesCOO() that would perform matrix assembly.
153: MATHYPRE is left out because this sequence trips a separate defect in hypre's device IJ
154: assemble, which reads the off-diagonal column map without mirroring it to the device
155: (hypre-space/hypre#1625). */
156: if (!isHypre) {
157: Mat D;
158: PetscInt rstart, rend, remote_col, ncoo;
159: PetscInt *coo_i, *coo_j;
160: PetscScalar readback;
161: PetscScalar *coo_v;
162: PetscBool isHostAIJ;
164: PetscCall(MatCreate(PETSC_COMM_WORLD, &D));
165: PetscCall(MatSetSizes(D, PETSC_DECIDE, PETSC_DECIDE, M, N));
166: PetscCall(MatSetFromOptions(D));
167: PetscCall(MatGetOwnershipRange(D, &rstart, &rend));
169: /* A diagonal entry and one in the first column past this rank's rows, which another rank owns
170: whenever there is more than one */
171: remote_col = rend % N;
172: ncoo = 2 * (rend - rstart);
173: PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j));
174: for (PetscInt r = rstart, k = 0; r < rend; r++) {
175: coo_i[k] = r;
176: coo_j[k] = r;
177: k++;
178: coo_i[k] = r;
179: coo_j[k] = remote_col;
180: k++;
181: }
182: PetscCall(MatSetPreallocationCOO(D, ncoo, coo_i, coo_j));
183: PetscCall(PetscFree2(coo_i, coo_j));
185: for (PetscInt r = rstart; r < rend; r++) {
186: PetscCall(MatSetValue(D, r, r, 2.0, INSERT_VALUES));
187: PetscCall(MatSetValue(D, r, remote_col, 5.0, INSERT_VALUES));
188: }
189: PetscCall(MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY));
190: PetscCall(MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY));
192: PetscCall(MatGetValue(D, rstart, remote_col, &readback));
193: if (readback != 5.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSetValues() after MatSetPreallocationCOO() gave %g, expected 5\n", (double)PetscRealPart(readback)));
195: /* MatSetValues() above only wrote into the preallocated nonzero pattern, so MatSetValuesCOO() remains usable */
196: PetscCall(PetscMalloc1(ncoo, &coo_v));
197: for (PetscInt k = 0; k < ncoo; k++) coo_v[k] = (k % 2) ? 7.0 : 0.0;
198: PetscCall(MatSetValuesCOO(D, coo_v, INSERT_VALUES));
199: PetscCall(MatGetValue(D, rstart, remote_col, &readback));
200: if (readback != 7.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSetValuesCOO() after MatSetValues() gave %g, expected 7\n", (double)PetscRealPart(readback)));
202: /* Zeroing no rows leaves the nonzero pattern, and so MatSetValuesCOO(), intact */
203: PetscCall(MatZeroRows(D, 0, NULL, 1.0, NULL, NULL));
204: PetscCall(MatSetValuesCOO(D, coo_v, ADD_VALUES));
205: PetscCall(MatGetValue(D, rstart, remote_col, &readback));
206: if (readback != 14.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSetValuesCOO() after an empty MatZeroRows() gave %g, expected 14\n", (double)PetscRealPart(readback)));
208: /* MatAXPY() with DIFFERENT_NONZERO_PATTERN of a matrix whose nonzeros are a subset of those of D leaves the nonzero pattern,
209: and so MatSetValuesCOO(), intact. Device implementations of MatAXPY(), such as that of MATAIJKOKKOS, replace the matrix,
210: so only the host types MATSEQAIJ and MATMPIAIJ are tested. */
211: PetscCall(PetscObjectTypeCompareAny((PetscObject)D, &isHostAIJ, MATSEQAIJ, MATMPIAIJ, ""));
212: if (isHostAIJ) {
213: Mat Diag;
215: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, rend - rstart, rend - rstart, M, N, 1, NULL, 0, NULL, &Diag));
216: for (PetscInt r = rstart; r < rend; r++) PetscCall(MatSetValue(Diag, r, r, 1.0, INSERT_VALUES));
217: PetscCall(MatAssemblyBegin(Diag, MAT_FINAL_ASSEMBLY));
218: PetscCall(MatAssemblyEnd(Diag, MAT_FINAL_ASSEMBLY));
219: PetscCall(MatAXPY(D, 1.0, Diag, DIFFERENT_NONZERO_PATTERN));
220: PetscCall(MatDestroy(&Diag));
221: PetscCall(MatSetValuesCOO(D, coo_v, INSERT_VALUES));
222: PetscCall(MatGetValue(D, rstart, remote_col, &readback));
223: if (readback != 7.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSetValuesCOO() after MatAXPY() gave %g, expected 7\n", (double)PetscRealPart(readback)));
224: }
225: PetscCall(PetscFree(coo_v));
226: PetscCall(MatDestroy(&D));
227: }
229: /* Test that aij->diag is set up correctly on COO matrices */
230: if (!isHypre) { // TODO: MATHYPRE currently does not support MatSetValues
231: PetscCall(MatShift(A, 2.0));
232: PetscCall(MatShift(B, 2.0));
233: PetscCall(MatMultEqual(A, B, 10, &equal));
234: if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift() on a duplicated COO matrix failed\n"));
236: /* A square matrix whose last process owns no rows, which on one process is a 0 x 0 matrix */
237: PetscCall(TestShiftCOO(PetscGlobalRank == size - 1 ? 0 : 3, PetscGlobalRank == size - 1 ? 0 : 3));
238: /* A matrix with more rows than columns */
239: PetscCall(TestShiftCOO(3, 2));
240: }
242: PetscCall(PetscFree(vals));
243: if (freecoo) PetscCall(PetscFree2(mycoo.i, mycoo.j));
244: PetscCall(MatDestroy(&A));
245: PetscCall(MatDestroy(&B));
246: PetscCall(MatDestroy(&C));
248: PetscCall(PetscFinalize());
249: return 0;
250: }
252: /*TEST
254: testset:
255: output_file: output/empty.out
256: nsize: {{1 2 3}}
257: args: -ignore_remote {{0 1}} -missing_diagonal {{0 1}}
258: filter: grep -v type | grep -v "Mat Object"
260: test:
261: suffix: kokkos
262: requires: kokkos_kernels
263: args: -mat_type aijkokkos
265: test:
266: suffix: cuda
267: requires: cuda
268: args: -mat_type aijcusparse
270: test:
271: suffix: hip
272: requires: hip
273: args: -mat_type aijhipsparse
275: test:
276: suffix: aij
277: args: -mat_type aij
279: test:
280: suffix: hypre
281: requires: hypre
282: args: -mat_type hypre
284: testset:
285: output_file: output/empty.out
286: nsize: 1
287: args: -ncoos 3
288: filter: grep -v type | grep -v "Mat Object"
290: test:
291: suffix: 2_kokkos
292: requires: kokkos_kernels
293: args: -mat_type aijkokkos
295: test:
296: suffix: 2_cuda
297: requires: cuda
298: args: -mat_type aijcusparse
300: test:
301: suffix: 2_hip
302: requires: hip
303: args: -mat_type aijhipsparse
305: test:
306: suffix: 2_aij
307: args: -mat_type aij
309: test:
310: suffix: 2_hypre
311: requires: hypre
312: args: -mat_type hypre
314: TEST*/