1 static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, At, AAt, T = NULL;
8 Vec x, y, z;
9 ISLocalToGlobalMapping rl2g, cl2g;
10 IS is;
11 PetscLayout rmap, cmap;
12 PetscInt *it, *jt;
13 PetscInt n1 = 11, n2 = 9;
14 PetscInt i1[] = {7, 6, 2, 0, 4, 1, 1, 0, 2, 2, 1, -1, -1};
15 PetscInt j1[] = {1, 4, 3, 5, 3, 3, 4, 5, 0, 3, 1, -1, -1};
16 PetscInt i2[] = {7, 6, 2, 0, 4, 1, 1, 2, 1, -1, -1};
17 PetscInt j2[] = {1, 4, 3, 5, 3, 3, 4, 0, 1, -1, -1};
18 PetscScalar v1[] = {-1., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., PETSC_MAX_REAL, PETSC_MAX_REAL};
19 PetscScalar v2[] = {1., -1., -2., -3., -4., -5., -6., -7., -8., -9., -10., PETSC_MAX_REAL, PETSC_MAX_REAL};
20 PetscInt N = 6, m = 8, M, rstart, cstart, i;
21 PetscMPIInt size;
22 PetscBool loc = PETSC_FALSE;
23 PetscBool locdiag = PETSC_TRUE;
24 PetscBool localapi = PETSC_FALSE;
25 PetscBool neg = PETSC_FALSE;
26 PetscBool ismatis, ismpiaij, ishypre;
27
28 PetscFunctionBeginUser;
29 PetscCall(PetscInitialize(&argc, &args, NULL, help));
30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-neg", &neg, NULL));
31 PetscCall(PetscOptionsGetBool(NULL, NULL, "-loc", &loc, NULL));
32 PetscCall(PetscOptionsGetBool(NULL, NULL, "-locdiag", &locdiag, NULL));
33 PetscCall(PetscOptionsGetBool(NULL, NULL, "-localapi", &localapi, NULL));
34 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
35 if (loc) {
36 if (locdiag) {
37 PetscCall(MatSetSizes(A, m, N, PETSC_DECIDE, PETSC_DECIDE));
38 } else {
39 PetscCall(MatSetSizes(A, m, m + N, PETSC_DECIDE, PETSC_DECIDE));
40 }
41 } else {
42 PetscCall(MatSetSizes(A, m, PETSC_DECIDE, PETSC_DECIDE, N));
43 }
44 PetscCall(MatSetFromOptions(A));
45 PetscCall(MatGetLayouts(A, &rmap, &cmap));
46 PetscCall(PetscLayoutSetUp(rmap));
47 PetscCall(PetscLayoutSetUp(cmap));
48 PetscCall(PetscLayoutGetRange(rmap, &rstart, NULL));
49 PetscCall(PetscLayoutGetRange(cmap, &cstart, NULL));
50 PetscCall(PetscLayoutGetSize(rmap, &M));
51 PetscCall(PetscLayoutGetSize(cmap, &N));
52
53 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
54 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
55
56 /* create fake l2g maps to test the local API */
57 PetscCall(ISCreateStride(PETSC_COMM_WORLD, M - rstart, rstart, 1, &is));
58 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
59 PetscCall(ISDestroy(&is));
60 PetscCall(ISCreateStride(PETSC_COMM_WORLD, N, 0, 1, &is));
61 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
62 PetscCall(ISDestroy(&is));
63 PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g));
64 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
65 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
66
67 PetscCall(MatCreateVecs(A, &x, &y));
68 PetscCall(MatCreateVecs(A, NULL, &z));
69 PetscCall(VecSet(x, 1.));
70 PetscCall(VecSet(z, 2.));
71 if (!localapi)
72 for (i = 0; i < n1; i++) i1[i] += rstart;
73 if (!localapi)
74 for (i = 0; i < n2; i++) i2[i] += rstart;
75 if (loc) {
76 if (locdiag) {
77 for (i = 0; i < n1; i++) j1[i] += cstart;
78 for (i = 0; i < n2; i++) j2[i] += cstart;
79 } else {
80 for (i = 0; i < n1; i++) j1[i] += cstart + m;
81 for (i = 0; i < n2; i++) j2[i] += cstart + m;
82 }
83 }
84 if (neg) {
85 n1 += 2;
86 n2 += 2;
87 }
88 /* MatSetPreallocationCOOLocal maps the indices! */
89 PetscCall(PetscMalloc2(PetscMax(n1, n2), &it, PetscMax(n1, n2), &jt));
90 /* test with repeated entries */
91 if (!localapi) {
92 PetscCall(MatSetPreallocationCOO(A, n1, i1, j1));
93 } else {
94 PetscCall(PetscArraycpy(it, i1, n1));
95 PetscCall(PetscArraycpy(jt, j1, n1));
96 PetscCall(MatSetPreallocationCOOLocal(A, n1, it, jt));
97 }
98 PetscCall(MatSetValuesCOO(A, v1, ADD_VALUES));
99 PetscCall(MatMult(A, x, y));
100 PetscCall(MatView(A, NULL));
101 PetscCall(VecView(y, NULL));
102 PetscCall(MatSetValuesCOO(A, v2, ADD_VALUES));
103 PetscCall(MatMultAdd(A, x, y, y));
104 PetscCall(MatView(A, NULL));
105 PetscCall(VecView(y, NULL));
106 T = A;
107 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
108 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At));
109 if (!ismatis) {
110 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
111 PetscCall(MatView(AAt, NULL));
112 PetscCall(MatDestroy(&AAt));
113 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
114 PetscCall(MatView(AAt, NULL));
115 PetscCall(MatDestroy(&AAt));
116 }
117 PetscCall(MatDestroy(&At));
118 if (ishypre) PetscCall(MatDestroy(&T));
119
120 /* INSERT_VALUES will overwrite matrix entries but
121 still perform the sum of the repeated entries */
122 PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES));
123 PetscCall(MatView(A, NULL));
124
125 /* test with unique entries */
126 PetscCall(PetscArraycpy(it, i2, n2));
127 PetscCall(PetscArraycpy(jt, j2, n2));
128 if (!localapi) {
129 PetscCall(MatSetPreallocationCOO(A, n2, it, jt));
130 } else {
131 PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt));
132 }
133 PetscCall(MatSetValuesCOO(A, v1, ADD_VALUES));
134 PetscCall(MatMult(A, x, y));
135 PetscCall(MatView(A, NULL));
136 PetscCall(VecView(y, NULL));
137 PetscCall(MatSetValuesCOO(A, v2, ADD_VALUES));
138 PetscCall(MatMultAdd(A, x, y, z));
139 PetscCall(MatView(A, NULL));
140 PetscCall(VecView(z, NULL));
141 PetscCall(PetscArraycpy(it, i2, n2));
142 PetscCall(PetscArraycpy(jt, j2, n2));
143 if (!localapi) {
144 PetscCall(MatSetPreallocationCOO(A, n2, it, jt));
145 } else {
146 PetscCall(MatSetPreallocationCOOLocal(A, n2, it, jt));
147 }
148 PetscCall(MatSetValuesCOO(A, v1, INSERT_VALUES));
149 PetscCall(MatMult(A, x, y));
150 PetscCall(MatView(A, NULL));
151 PetscCall(VecView(y, NULL));
152 PetscCall(MatSetValuesCOO(A, v2, INSERT_VALUES));
153 PetscCall(MatMultAdd(A, x, y, z));
154 PetscCall(MatView(A, NULL));
155 PetscCall(VecView(z, NULL));
156 T = A;
157 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
158 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At));
159 if (!ismatis) {
160 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
161 PetscCall(MatView(AAt, NULL));
162 PetscCall(MatDestroy(&AAt));
163 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
164 PetscCall(MatView(AAt, NULL));
165 PetscCall(MatDestroy(&AAt));
166 }
167 PetscCall(MatDestroy(&At));
168 if (ishypre) PetscCall(MatDestroy(&T));
169
170 /* test providing diagonal first, then off-diagonal */
171 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
172 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
173 if ((ismpiaij || ishypre) && size > 1) {
174 Mat lA, lB;
175 const PetscInt *garray, *iA, *jA, *iB, *jB;
176 const PetscScalar *vA, *vB;
177 PetscScalar *coo_v;
178 PetscInt *coo_i, *coo_j;
179 PetscInt i, j, nA, nB, nnz;
180 PetscBool flg;
181
182 T = A;
183 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
184 PetscCall(MatMPIAIJGetSeqAIJ(T, &lA, &lB, &garray));
185 PetscCall(MatSeqAIJGetArrayRead(lA, &vA));
186 PetscCall(MatSeqAIJGetArrayRead(lB, &vB));
187 PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg));
188 PetscCall(MatGetRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg));
189 nnz = iA[nA] + iB[nB];
190 PetscCall(PetscMalloc3(nnz, &coo_i, nnz, &coo_j, nnz, &coo_v));
191 nnz = 0;
192 for (i = 0; i < nA; i++) {
193 for (j = iA[i]; j < iA[i + 1]; j++, nnz++) {
194 coo_i[nnz] = i + rstart;
195 coo_j[nnz] = jA[j] + cstart;
196 coo_v[nnz] = vA[j];
197 }
198 }
199 for (i = 0; i < nB; i++) {
200 for (j = iB[i]; j < iB[i + 1]; j++, nnz++) {
201 coo_i[nnz] = i + rstart;
202 coo_j[nnz] = garray[jB[j]];
203 coo_v[nnz] = vB[j];
204 }
205 }
206 PetscCall(MatRestoreRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg));
207 PetscCall(MatRestoreRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg));
208 PetscCall(MatSeqAIJRestoreArrayRead(lA, &vA));
209 PetscCall(MatSeqAIJRestoreArrayRead(lB, &vB));
210 if (ishypre) PetscCall(MatDestroy(&T));
211
212 PetscCall(MatSetPreallocationCOO(A, nnz, coo_i, coo_j));
213 PetscCall(MatSetValuesCOO(A, coo_v, ADD_VALUES));
214 PetscCall(MatMult(A, x, y));
215 PetscCall(MatView(A, NULL));
216 PetscCall(VecView(y, NULL));
217 PetscCall(MatSetValuesCOO(A, coo_v, INSERT_VALUES));
218 PetscCall(MatMult(A, x, y));
219 PetscCall(MatView(A, NULL));
220 PetscCall(VecView(y, NULL));
221
222 T = A;
223 if (ishypre) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T));
224 PetscCall(MatTranspose(T, MAT_INITIAL_MATRIX, &At));
225 PetscCall(MatMatMult(T, At, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
226 PetscCall(MatView(AAt, NULL));
227 PetscCall(MatDestroy(&AAt));
228 PetscCall(MatMatMult(At, T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AAt));
229 PetscCall(MatView(AAt, NULL));
230 PetscCall(MatDestroy(&AAt));
231 PetscCall(MatDestroy(&At));
232 if (ishypre) PetscCall(MatDestroy(&T));
233
234 PetscCall(PetscFree3(coo_i, coo_j, coo_v));
235 }
236 PetscCall(PetscFree2(it, jt));
237 PetscCall(VecDestroy(&z));
238 PetscCall(VecDestroy(&x));
239 PetscCall(VecDestroy(&y));
240 PetscCall(MatDestroy(&A));
241 PetscCall(PetscFinalize());
242 return 0;
243 }
244
245 /*TEST
246
247 test:
248 suffix: 1
249 filter: grep -v type | grep -v "Mat Object"
250 diff_args: -j
251 args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}}
252
253 test:
254 requires: hypre
255 suffix: 1_hypre
256 filter: grep -v type | grep -v "Mat Object"
257 diff_args: -j
258 args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}}
259 output_file: output/ex123_1.out
260
261 test:
262 requires: cuda
263 suffix: 1_cuda
264 filter: grep -v type | grep -v "Mat Object"
265 diff_args: -j
266 args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}}
267 output_file: output/ex123_1.out
268
269 test:
270 requires: kokkos_kernels
271 suffix: 1_kokkos
272 filter: grep -v type | grep -v "Mat Object"
273 diff_args: -j
274 args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}}
275 output_file: output/ex123_1.out
276
277 test:
278 suffix: 2
279 nsize: 7
280 filter: grep -v type | grep -v "Mat Object"
281 diff_args: -j
282 args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}}
283
284 test:
285 requires: hypre
286 suffix: 2_hypre
287 nsize: 7
288 filter: grep -v type | grep -v "Mat Object"
289 diff_args: -j
290 args: -mat_type hypre -localapi {{0 1}} -neg {{0 1}}
291 output_file: output/ex123_2.out
292
293 test:
294 requires: cuda
295 suffix: 2_cuda
296 nsize: 7
297 filter: grep -v type | grep -v "Mat Object"
298 diff_args: -j
299 args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}}
300 output_file: output/ex123_2.out
301
302 test:
303 requires: kokkos_kernels
304 suffix: 2_kokkos
305 nsize: 7
306 filter: grep -v type | grep -v "Mat Object"
307 diff_args: -j
308 args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}}
309 output_file: output/ex123_2.out
310
311 test:
312 suffix: 3
313 nsize: 3
314 filter: grep -v type | grep -v "Mat Object"
315 diff_args: -j
316 args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}}
317
318 test:
319 requires: hypre
320 suffix: 3_hypre
321 nsize: 3
322 filter: grep -v type | grep -v "Mat Object"
323 diff_args: -j
324 args: -mat_type hypre -loc -localapi {{0 1}} -neg {{0 1}}
325 output_file: output/ex123_3.out
326
327 test:
328 requires: cuda
329 suffix: 3_cuda
330 nsize: 3
331 filter: grep -v type | grep -v "Mat Object"
332 diff_args: -j
333 args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}}
334 output_file: output/ex123_3.out
335
336 test:
337 requires: kokkos_kernels
338 suffix: 3_kokkos
339 nsize: 3
340 filter: grep -v type | grep -v "Mat Object"
341 diff_args: -j
342 args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}}
343 output_file: output/ex123_3.out
344
345 test:
346 suffix: 4
347 nsize: 4
348 filter: grep -v type | grep -v "Mat Object"
349 diff_args: -j
350 args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
351
352 test:
353 requires: hypre
354 suffix: 4_hypre
355 nsize: 4
356 filter: grep -v type | grep -v "Mat Object"
357 diff_args: -j
358 args: -mat_type hypre -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
359 output_file: output/ex123_4.out
360
361 test:
362 requires: cuda
363 suffix: 4_cuda
364 nsize: 4
365 filter: grep -v type | grep -v "Mat Object"
366 diff_args: -j
367 args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
368 output_file: output/ex123_4.out
369
370 test:
371 requires: kokkos_kernels
372 suffix: 4_kokkos
373 nsize: 4
374 filter: grep -v type | grep -v "Mat Object"
375 diff_args: -j
376 args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
377 output_file: output/ex123_4.out
378
379 test:
380 suffix: matis
381 nsize: 3
382 filter: grep -v type | grep -v "Mat Object"
383 diff_args: -j
384 args: -mat_type is -localapi {{0 1}} -neg {{0 1}}
385
386 TEST*/
387