Actual source code: ex2.c
1: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
3: #include <petscmat.h>
5: static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
6: {
7: Mat D, E, F, G;
8: MatType mtype;
10: PetscFunctionBegin;
11: PetscCall(MatGetType(mat, &mtype));
12: if (f == MatCreateTranspose) {
13: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
14: } else {
15: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
16: }
17: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
18: PetscCall(f(C, &D));
19: PetscCall(f(D, &E));
20: PetscCall(MatAXPY(E, alpha, mat, SAME_NONZERO_PATTERN));
21: PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
22: PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
23: PetscCall(MatDestroy(&E));
24: PetscCall(MatDestroy(&D));
25: PetscCall(MatDestroy(&C));
26: if (f == MatCreateTranspose) {
27: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
28: } else {
29: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
30: }
31: if (f == MatCreateTranspose) {
32: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &D));
33: } else {
34: PetscCall(MatHermitianTranspose(mat, MAT_INITIAL_MATRIX, &D));
35: }
36: PetscCall(f(D, &E));
37: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
38: PetscCall(MatAXPY(C, alpha, E, SAME_NONZERO_PATTERN));
39: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
40: PetscCall(MatDestroy(&E));
41: PetscCall(MatDestroy(&D));
42: PetscCall(MatDestroy(&C));
43: if (f == MatCreateTranspose) {
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
45: } else {
46: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
47: }
48: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
49: PetscCall(f(C, &D));
50: PetscCall(f(D, &E));
51: PetscCall(f(mat, &F));
52: PetscCall(f(F, &G));
53: PetscCall(MatAXPY(E, alpha, G, SAME_NONZERO_PATTERN));
54: PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
55: PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
56: PetscCall(MatDestroy(&G));
57: PetscCall(MatDestroy(&F));
58: PetscCall(MatDestroy(&E));
59: PetscCall(MatDestroy(&D));
60: PetscCall(MatDestroy(&C));
62: /* Call f on a matrix that does not implement the transposition */
63: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: Now without the transposition operation\n"));
64: PetscCall(MatConvert(mat, MATSHELL, MAT_INITIAL_MATRIX, &C));
65: PetscCall(f(C, &D));
66: PetscCall(f(D, &E));
67: /* XXX cannot use MAT_INPLACE_MATRIX, it leaks mat */
68: PetscCall(MatConvert(E, mtype, MAT_INITIAL_MATRIX, &F));
69: PetscCall(MatAXPY(F, alpha, mat, SAME_NONZERO_PATTERN));
70: PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
71: PetscCall(MatDestroy(&F));
72: PetscCall(MatDestroy(&E));
73: PetscCall(MatDestroy(&D));
74: PetscCall(MatDestroy(&C));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: int main(int argc, char **argv)
79: {
80: Mat mat, tmat = NULL;
81: PetscInt m = 7, n, i, j, rstart, rend, rect = 0;
82: PetscMPIInt size, rank;
83: PetscBool flg;
84: PetscScalar v, alpha;
85: PetscReal normf, normi, norm1;
87: PetscFunctionBeginUser;
88: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
89: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
90: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
91: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
92: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
93: n = m;
94: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
95: if (flg) {
96: n += 2;
97: rect = 1;
98: }
99: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
100: if (flg) {
101: n -= 2;
102: rect = 1;
103: }
105: /* ------- Assemble matrix --------- */
106: PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
107: PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
108: PetscCall(MatSetFromOptions(mat));
109: PetscCall(MatSetUp(mat));
110: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
111: for (i = rstart; i < rend; i++) {
112: for (j = 0; j < n; j++) {
113: v = 10.0 * i + j + 1.0;
114: PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
115: }
116: }
117: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
118: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
120: /* ----------------- Test MatNorm() ----------------- */
121: PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
122: PetscCall(MatNorm(mat, NORM_1, &norm1));
123: PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
124: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original A: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
125: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
126: {
127: /* The printed norms are masked by petscdiff, so compare the norms of -A with those of a host MATAIJ copy of A.
128: The entries of -A are negative, so the norms of -A only match if absolute values are taken */
129: const NormType types[] = {NORM_FROBENIUS, NORM_1, NORM_INFINITY};
130: Mat C, D;
131: PetscReal nc, nd;
133: PetscCall(MatConvert(mat, MATAIJ, MAT_INITIAL_MATRIX, &C));
134: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &D));
135: PetscCall(MatScale(D, -1.0));
136: for (PetscInt t = 0; t < 3; t++) {
137: PetscCall(MatNorm(C, types[t], &nc));
138: PetscCall(MatNorm(D, types[t], &nd));
139: PetscCheck(PetscIsCloseAtTol(nd, nc, PETSC_SMALL, 0.0), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "NORM_%s of -A is %g, but %g for a MATAIJ copy of A", NormTypes[types[t]], (double)nd, (double)nc);
140: }
141: PetscCall(MatHasOperation(D, MATOP_GET_ROW_SUM_ABS, &flg));
142: if (flg) {
143: Vec vc, vd;
145: PetscCall(MatCreateVecs(C, NULL, &vc));
146: PetscCall(MatCreateVecs(D, NULL, &vd));
147: PetscCall(MatGetRowSumAbs(C, vc));
148: PetscCall(MatGetRowSumAbs(D, vd));
149: PetscCall(VecAXPY(vd, -1.0, vc));
150: PetscCall(VecNorm(vd, NORM_INFINITY, &nd));
151: PetscCheck(nd <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatGetRowSumAbs() of -A differs by %g from that of a MATAIJ copy of A", (double)nd);
152: PetscCall(VecDestroy(&vd));
153: PetscCall(VecDestroy(&vc));
154: }
155: PetscCall(MatDestroy(&D));
156: PetscCall(MatDestroy(&C));
157: }
158: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_comm_self", &flg));
159: if (flg) {
160: /* Each rank builds its own matrix of the requested type on PETSC_COMM_SELF, so a parallel type has a single-rank
161: communicator even when PETSC_COMM_WORLD has several. The entries are negative so that the row sums only match
162: if absolute values are taken, and they depend on the rank to keep the matrices distinct */
163: Mat C, D;
164: Vec vc, vd;
165: PetscReal diff;
167: PetscCall(MatCreate(PETSC_COMM_SELF, &D));
168: PetscCall(MatSetSizes(D, m, n, m, n));
169: PetscCall(MatSetFromOptions(D));
170: PetscCall(MatSetUp(D));
171: for (i = 0; i < m; i++) {
172: for (j = 0; j < n; j++) {
173: v = -(i + j + rank + 1.0);
174: PetscCall(MatSetValues(D, 1, &i, 1, &j, &v, INSERT_VALUES));
175: }
176: }
177: PetscCall(MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY));
178: PetscCall(MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY));
179: PetscCall(MatConvert(D, MATSEQAIJ, MAT_INITIAL_MATRIX, &C));
180: PetscCall(MatCreateVecs(C, NULL, &vc));
181: PetscCall(MatCreateVecs(D, NULL, &vd));
182: PetscCall(MatGetRowSumAbs(C, vc));
183: PetscCall(MatGetRowSumAbs(D, vd));
184: PetscCall(VecAXPY(vd, -1.0, vc));
185: PetscCall(VecNorm(vd, NORM_INFINITY, &diff));
186: PetscCheck(diff <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRowSumAbs() of a PETSC_COMM_SELF matrix differs by %g from that of a MATSEQAIJ copy", (double)diff);
187: PetscCall(VecDestroy(&vd));
188: PetscCall(VecDestroy(&vc));
189: PetscCall(MatDestroy(&D));
190: PetscCall(MatDestroy(&C));
191: }
193: /* --------------- Test MatTranspose() -------------- */
194: PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
195: if (!rect && flg) {
196: PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
197: tmat = mat;
198: mat = NULL;
199: } else { /* out-of-place transpose */
200: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
201: }
203: /* ----------------- Test MatNorm() ----------------- */
204: /* Print info about transpose matrix */
205: PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
206: PetscCall(MatNorm(tmat, NORM_1, &norm1));
207: PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
208: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
209: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
211: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
212: if (mat && !rect) {
213: alpha = 1.0;
214: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
215: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A\n"));
216: PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
217: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
219: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX: B = alpha*B + A\n"));
220: PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
221: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
223: {
224: Mat A, B;
225: PetscBool equal;
227: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &A));
228: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &B));
229: PetscCall(MatAYPX(A, 2.0, A, SAME_NONZERO_PATTERN));
230: PetscCall(MatScale(B, 3.0));
231: PetscCall(MatEqual(A, B, &equal));
232: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatAYPX() failed when Y == X");
233: PetscCall(MatDestroy(&A));
234: PetscCall(MatDestroy(&B));
235: }
236: }
238: {
239: Mat C;
240: alpha = 1.0;
241: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
242: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
243: PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
244: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
245: PetscCall(MatDestroy(&C));
246: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
247: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
248: }
250: {
251: Mat matB;
252: /* get matB that has nonzeros of mat in all even numbers of row and col */
253: PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
254: PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
255: PetscCall(MatSetFromOptions(matB));
256: PetscCall(MatSetUp(matB));
257: PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
258: if (rstart % 2 != 0) rstart++;
259: for (i = rstart; i < rend; i += 2) {
260: for (j = 0; j < n; j += 2) {
261: v = 10.0 * i + j + 1.0;
262: PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
263: }
264: }
265: PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
266: PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
267: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
268: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
269: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
270: PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
271: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
272: PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
273: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
274: PetscCall(MatDestroy(&matB));
275: }
277: /* Test MatZeroRows */
278: j = rstart - 1;
279: if (j < 0) j = m - 1;
280: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
281: PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
282: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
284: /* Test MatShift */
285: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
286: PetscCall(MatShift(mat, -2.0));
287: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
289: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
290: /* Free data structures */
291: PetscCall(MatDestroy(&mat));
292: PetscCall(MatDestroy(&tmat));
293: PetscCall(PetscFinalize());
294: return 0;
295: }
297: /*TEST
299: test:
300: suffix: 11_A
301: args: -mat_type seqaij -rectA
302: filter: grep -v "Mat Object"
304: test:
305: suffix: 12_A
306: args: -mat_type seqdense -rectA
307: filter: grep -v type | grep -v "Mat Object"
309: test:
310: requires: cuda
311: suffix: 12_A_cuda
312: args: -mat_type seqdensecuda -rectA
313: output_file: output/ex2_12_A.out
314: filter: grep -v type | grep -v "Mat Object"
316: test:
317: requires: kokkos_kernels
318: suffix: 12_A_kokkos
319: args: -mat_type aijkokkos -rectA
320: output_file: output/ex2_12_A.out
321: filter: grep -v type | grep -v "Mat Object"
323: test:
324: suffix: 11_B
325: args: -mat_type seqaij -rectB
326: filter: grep -v "Mat Object"
328: test:
329: suffix: 12_B
330: args: -mat_type seqdense -rectB
331: filter: grep -v type | grep -v "Mat Object"
333: testset:
334: args: -rectB
335: output_file: output/ex2_12_B.out
336: filter: grep -v type | grep -v "Mat Object"
338: test:
339: requires: cuda
340: suffix: 12_B_cuda
341: args: -mat_type {{seqdensecuda seqaijcusparse}}
343: test:
344: requires: kokkos_kernels
345: suffix: 12_B_kokkos
346: args: -mat_type aijkokkos
348: test:
349: suffix: 12_B_aij
350: args: -mat_type aij
351: test:
352: suffix: 21
353: args: -mat_type mpiaij
354: filter: grep -v type | grep -v " MPI process"
356: test:
357: suffix: 22
358: args: -mat_type mpidense
359: filter: grep -v type | grep -v "Mat Object"
361: test:
362: requires: cuda
363: suffix: 22_cuda
364: output_file: output/ex2_22.out
365: args: -mat_type mpidensecuda
366: filter: grep -v type | grep -v "Mat Object"
368: test:
369: requires: kokkos_kernels
370: suffix: 22_kokkos
371: output_file: output/ex2_22.out
372: args: -mat_type aijkokkos
373: filter: grep -v type | grep -v "Mat Object"
375: test:
376: suffix: 23
377: nsize: 3
378: args: -mat_type mpiaij
379: filter: grep -v type | grep -v " MPI process"
381: test:
382: suffix: 24
383: nsize: 3
384: args: -mat_type mpidense
385: filter: grep -v type | grep -v "Mat Object"
387: test:
388: requires: cuda
389: suffix: 24_cuda
390: nsize: 3
391: output_file: output/ex2_24.out
392: args: -mat_type mpidensecuda
393: filter: grep -v type | grep -v "Mat Object"
395: test:
396: suffix: 2_aijcusparse_1
397: args: -mat_type mpiaijcusparse
398: output_file: output/ex2_21.out
399: requires: cuda
400: filter: grep -v type | grep -v " MPI process"
402: test:
403: suffix: 2_aijkokkos_1
404: args: -mat_type aijkokkos
405: output_file: output/ex2_21.out
406: requires: kokkos_kernels
407: filter: grep -v type | grep -v " MPI process"
409: test:
410: suffix: 2_mpiaijkokkos_comm_self
411: nsize: 3
412: args: -mat_type mpiaijkokkos -test_comm_self
413: output_file: output/ex2_23.out
414: requires: kokkos_kernels
415: filter: grep -v type | grep -v "MPI processes"
417: test:
418: suffix: 2_aijcusparse_2
419: nsize: 3
420: args: -mat_type mpiaijcusparse
421: output_file: output/ex2_23.out
422: requires: cuda
423: filter: grep -v type | grep -v " MPI process"
425: test:
426: suffix: 2_aijkokkos_2
427: nsize: 3
428: args: -mat_type aijkokkos
429: output_file: output/ex2_23.out
430: # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
431: requires: !hip kokkos_kernels
432: filter: grep -v type | grep -v "MPI processes"
434: test:
435: suffix: 3
436: nsize: 2
437: args: -mat_type mpiaij -rectA
439: test:
440: suffix: 3_aijkokkos
441: nsize: 2
442: args: -mat_type aijkokkos -rectA
443: output_file: output/ex2_3.out
444: requires: kokkos_kernels
445: filter: sed -e "s/mpiaijkokkos/mpiaij/"
447: test:
448: suffix: 3_aijcusparse
449: nsize: 2
450: args: -mat_type mpiaijcusparse -rectA
451: requires: cuda
453: test:
454: suffix: 4
455: nsize: 2
456: args: -mat_type mpidense -rectA
457: filter: grep -v type | grep -v " MPI process"
459: test:
460: requires: cuda
461: suffix: 4_cuda
462: nsize: 2
463: output_file: output/ex2_4.out
464: args: -mat_type mpidensecuda -rectA
465: filter: grep -v type | grep -v " MPI process"
467: test:
468: suffix: aijcusparse_1
469: args: -mat_type seqaijcusparse -rectA
470: filter: grep -v "Mat Object"
471: output_file: output/ex2_11_A_aijcusparse.out
472: requires: cuda
474: TEST*/