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