Actual source code: ex123.c
1: static char help[] = "Test MatSetPreallocationCOO and MatSetValuesCOO\n\n";
3: #include <petscmat.h>
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;
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));
53: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
54: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
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));
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));
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));
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));
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;
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));
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));
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));
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: }
245: /*TEST
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}}
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
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
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
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}}
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
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
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
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}}
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
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
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
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}}
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
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
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
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}}
386: TEST*/