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;
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;
29: PetscInitialize(&argc, &args, (char *)0, help);
30: PetscOptionsGetBool(NULL, NULL, "-neg", &neg, NULL);
31: PetscOptionsGetBool(NULL, NULL, "-loc", &loc, NULL);
32: PetscOptionsGetBool(NULL, NULL, "-locdiag", &locdiag, NULL);
33: PetscOptionsGetBool(NULL, NULL, "-localapi", &localapi, NULL);
34: MatCreate(PETSC_COMM_WORLD, &A);
35: if (loc) {
36: if (locdiag) {
37: MatSetSizes(A, m, N, PETSC_DECIDE, PETSC_DECIDE);
38: } else {
39: MatSetSizes(A, m, m + N, PETSC_DECIDE, PETSC_DECIDE);
40: }
41: } else {
42: MatSetSizes(A, m, PETSC_DECIDE, PETSC_DECIDE, N);
43: }
44: MatSetFromOptions(A);
45: MatGetLayouts(A, &rmap, &cmap);
46: PetscLayoutSetUp(rmap);
47: PetscLayoutSetUp(cmap);
48: PetscLayoutGetRange(rmap, &rstart, NULL);
49: PetscLayoutGetRange(cmap, &cstart, NULL);
50: PetscLayoutGetSize(rmap, &M);
51: PetscLayoutGetSize(cmap, &N);
53: PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis);
55: /* create fake l2g maps to test the local API */
56: ISCreateStride(PETSC_COMM_WORLD, M - rstart, rstart, 1, &is);
57: ISLocalToGlobalMappingCreateIS(is, &rl2g);
58: ISDestroy(&is);
59: ISCreateStride(PETSC_COMM_WORLD, N, 0, 1, &is);
60: ISLocalToGlobalMappingCreateIS(is, &cl2g);
61: ISDestroy(&is);
62: MatSetLocalToGlobalMapping(A, rl2g, cl2g);
63: ISLocalToGlobalMappingDestroy(&rl2g);
64: ISLocalToGlobalMappingDestroy(&cl2g);
66: MatCreateVecs(A, &x, &y);
67: MatCreateVecs(A, NULL, &z);
68: VecSet(x, 1.);
69: VecSet(z, 2.);
70: if (!localapi)
71: for (i = 0; i < n1; i++) i1[i] += rstart;
72: if (!localapi)
73: for (i = 0; i < n2; i++) i2[i] += rstart;
74: if (loc) {
75: if (locdiag) {
76: for (i = 0; i < n1; i++) j1[i] += cstart;
77: for (i = 0; i < n2; i++) j2[i] += cstart;
78: } else {
79: for (i = 0; i < n1; i++) j1[i] += cstart + m;
80: for (i = 0; i < n2; i++) j2[i] += cstart + m;
81: }
82: }
83: if (neg) {
84: n1 += 2;
85: n2 += 2;
86: }
87: /* MatSetPreallocationCOOLocal maps the indices! */
88: PetscMalloc2(PetscMax(n1, n2), &it, PetscMax(n1, n2), &jt);
89: /* test with repeated entries */
90: if (!localapi) {
91: MatSetPreallocationCOO(A, n1, i1, j1);
92: } else {
93: PetscArraycpy(it, i1, n1);
94: PetscArraycpy(jt, j1, n1);
95: MatSetPreallocationCOOLocal(A, n1, it, jt);
96: }
97: MatSetValuesCOO(A, v1, ADD_VALUES);
98: MatMult(A, x, y);
99: MatView(A, NULL);
100: VecView(y, NULL);
101: MatSetValuesCOO(A, v2, ADD_VALUES);
102: MatMultAdd(A, x, y, y);
103: MatView(A, NULL);
104: VecView(y, NULL);
105: MatTranspose(A, MAT_INITIAL_MATRIX, &At);
106: if (!ismatis) {
107: MatMatMult(A, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
108: MatView(AAt, NULL);
109: MatDestroy(&AAt);
110: MatMatMult(At, A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
111: MatView(AAt, NULL);
112: MatDestroy(&AAt);
113: }
114: MatDestroy(&At);
116: /* INSERT_VALUES will overwrite matrix entries but
117: still perform the sum of the repeated entries */
118: MatSetValuesCOO(A, v2, INSERT_VALUES);
119: MatView(A, NULL);
121: /* test with unique entries */
122: PetscArraycpy(it, i2, n2);
123: PetscArraycpy(jt, j2, n2);
124: if (!localapi) {
125: MatSetPreallocationCOO(A, n2, it, jt);
126: } else {
127: MatSetPreallocationCOOLocal(A, n2, it, jt);
128: }
129: MatSetValuesCOO(A, v1, ADD_VALUES);
130: MatMult(A, x, y);
131: MatView(A, NULL);
132: VecView(y, NULL);
133: MatSetValuesCOO(A, v2, ADD_VALUES);
134: MatMultAdd(A, x, y, z);
135: MatView(A, NULL);
136: VecView(z, NULL);
137: PetscArraycpy(it, i2, n2);
138: PetscArraycpy(jt, j2, n2);
139: if (!localapi) {
140: MatSetPreallocationCOO(A, n2, it, jt);
141: } else {
142: MatSetPreallocationCOOLocal(A, n2, it, jt);
143: }
144: MatSetValuesCOO(A, v1, INSERT_VALUES);
145: MatMult(A, x, y);
146: MatView(A, NULL);
147: VecView(y, NULL);
148: MatSetValuesCOO(A, v2, INSERT_VALUES);
149: MatMultAdd(A, x, y, z);
150: MatView(A, NULL);
151: VecView(z, NULL);
152: MatTranspose(A, MAT_INITIAL_MATRIX, &At);
153: if (!ismatis) {
154: MatMatMult(A, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
155: MatView(AAt, NULL);
156: MatDestroy(&AAt);
157: MatMatMult(At, A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
158: MatView(AAt, NULL);
159: MatDestroy(&AAt);
160: }
161: MatDestroy(&At);
163: /* test providing diagonal first, then offdiagonal */
164: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
165: PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij);
166: if (ismpiaij && size > 1) {
167: Mat lA, lB;
168: const PetscInt *garray, *iA, *jA, *iB, *jB;
169: const PetscScalar *vA, *vB;
170: PetscScalar *coo_v;
171: PetscInt *coo_i, *coo_j;
172: PetscInt i, j, nA, nB, nnz;
173: PetscBool flg;
175: MatMPIAIJGetSeqAIJ(A, &lA, &lB, &garray);
176: MatSeqAIJGetArrayRead(lA, &vA);
177: MatSeqAIJGetArrayRead(lB, &vB);
178: MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg);
179: MatGetRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg);
180: nnz = iA[nA] + iB[nB];
181: PetscMalloc3(nnz, &coo_i, nnz, &coo_j, nnz, &coo_v);
182: nnz = 0;
183: for (i = 0; i < nA; i++) {
184: for (j = iA[i]; j < iA[i + 1]; j++, nnz++) {
185: coo_i[nnz] = i + rstart;
186: coo_j[nnz] = jA[j] + cstart;
187: coo_v[nnz] = vA[j];
188: }
189: }
190: for (i = 0; i < nB; i++) {
191: for (j = iB[i]; j < iB[i + 1]; j++, nnz++) {
192: coo_i[nnz] = i + rstart;
193: coo_j[nnz] = garray[jB[j]];
194: coo_v[nnz] = vB[j];
195: }
196: }
197: MatRestoreRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nA, &iA, &jA, &flg);
198: MatRestoreRowIJ(lB, 0, PETSC_FALSE, PETSC_FALSE, &nB, &iB, &jB, &flg);
199: MatSeqAIJRestoreArrayRead(lA, &vA);
200: MatSeqAIJRestoreArrayRead(lB, &vB);
202: MatSetPreallocationCOO(A, nnz, coo_i, coo_j);
203: MatSetValuesCOO(A, coo_v, ADD_VALUES);
204: MatMult(A, x, y);
205: MatView(A, NULL);
206: VecView(y, NULL);
207: MatSetValuesCOO(A, coo_v, INSERT_VALUES);
208: MatMult(A, x, y);
209: MatView(A, NULL);
210: VecView(y, NULL);
211: MatTranspose(A, MAT_INITIAL_MATRIX, &At);
212: MatMatMult(A, At, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
213: MatView(AAt, NULL);
214: MatDestroy(&AAt);
215: MatMatMult(At, A, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AAt);
216: MatView(AAt, NULL);
217: MatDestroy(&AAt);
218: MatDestroy(&At);
220: PetscFree3(coo_i, coo_j, coo_v);
221: }
222: PetscFree2(it, jt);
223: VecDestroy(&z);
224: VecDestroy(&x);
225: VecDestroy(&y);
226: MatDestroy(&A);
227: PetscFinalize();
228: return 0;
229: }
231: /*TEST
233: test:
234: suffix: 1
235: filter: grep -v type
236: diff_args: -j
237: args: -mat_type {{seqaij mpiaij}} -localapi {{0 1}} -neg {{0 1}}
239: test:
240: requires: cuda
241: suffix: 1_cuda
242: filter: grep -v type
243: diff_args: -j
244: args: -mat_type {{seqaijcusparse mpiaijcusparse}} -localapi {{0 1}} -neg {{0 1}}
245: output_file: output/ex123_1.out
247: test:
248: requires: kokkos_kernels
249: suffix: 1_kokkos
250: filter: grep -v type
251: diff_args: -j
252: args: -mat_type {{seqaijkokkos mpiaijkokkos}} -localapi {{0 1}} -neg {{0 1}}
253: output_file: output/ex123_1.out
255: test:
256: suffix: 2
257: nsize: 7
258: filter: grep -v type
259: diff_args: -j
260: args: -mat_type mpiaij -localapi {{0 1}} -neg {{0 1}}
262: test:
263: requires: cuda
264: suffix: 2_cuda
265: nsize: 7
266: filter: grep -v type
267: diff_args: -j
268: args: -mat_type mpiaijcusparse -localapi {{0 1}} -neg {{0 1}}
269: output_file: output/ex123_2.out
271: test:
272: requires: kokkos_kernels
273: suffix: 2_kokkos
274: nsize: 7
275: filter: grep -v type
276: diff_args: -j
277: args: -mat_type mpiaijkokkos -localapi {{0 1}} -neg {{0 1}}
278: output_file: output/ex123_2.out
280: test:
281: suffix: 3
282: nsize: 3
283: filter: grep -v type
284: diff_args: -j
285: args: -mat_type mpiaij -loc -localapi {{0 1}} -neg {{0 1}}
287: test:
288: requires: cuda
289: suffix: 3_cuda
290: nsize: 3
291: filter: grep -v type
292: diff_args: -j
293: args: -mat_type mpiaijcusparse -loc -localapi {{0 1}} -neg {{0 1}}
294: output_file: output/ex123_3.out
296: test:
297: requires: kokkos_kernels
298: suffix: 3_kokkos
299: nsize: 3
300: filter: grep -v type
301: diff_args: -j
302: args: -mat_type aijkokkos -loc -localapi {{0 1}} -neg {{0 1}}
303: output_file: output/ex123_3.out
305: test:
306: suffix: 4
307: nsize: 4
308: filter: grep -v type
309: diff_args: -j
310: args: -mat_type mpiaij -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
312: test:
313: requires: cuda
314: suffix: 4_cuda
315: nsize: 4
316: filter: grep -v type
317: diff_args: -j
318: args: -mat_type mpiaijcusparse -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
319: output_file: output/ex123_4.out
321: test:
322: requires: kokkos_kernels
323: suffix: 4_kokkos
324: nsize: 4
325: filter: grep -v type
326: diff_args: -j
327: args: -mat_type aijkokkos -loc -locdiag 0 -localapi {{0 1}} -neg {{0 1}}
328: output_file: output/ex123_4.out
330: test:
331: suffix: matis
332: nsize: 3
333: filter: grep -v type
334: diff_args: -j
335: args: -mat_type is -localapi {{0 1}} -neg {{0 1}}
337: TEST*/