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 = 0;
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));
127: /* --------------- Test MatTranspose() -------------- */
128: PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
129: if (!rect && flg) {
130: PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
131: tmat = mat;
132: mat = NULL;
133: } else { /* out-of-place transpose */
134: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
135: }
137: /* ----------------- Test MatNorm() ----------------- */
138: /* Print info about transpose matrix */
139: PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
140: PetscCall(MatNorm(tmat, NORM_1, &norm1));
141: PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
143: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
145: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
146: if (mat && !rect) {
147: alpha = 1.0;
148: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
149: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A\n"));
150: PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
151: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
153: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX: B = alpha*B + A\n"));
154: PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
155: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
157: {
158: Mat A, B;
159: PetscBool equal;
161: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &A));
162: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &B));
163: PetscCall(MatAYPX(A, 2.0, A, SAME_NONZERO_PATTERN));
164: PetscCall(MatScale(B, 3.0));
165: PetscCall(MatEqual(A, B, &equal));
166: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatAYPX() failed when Y == X");
167: PetscCall(MatDestroy(&A));
168: PetscCall(MatDestroy(&B));
169: }
170: }
172: {
173: Mat C;
174: alpha = 1.0;
175: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
176: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
177: PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
178: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
179: PetscCall(MatDestroy(&C));
180: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
181: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
182: }
184: {
185: Mat matB;
186: /* get matB that has nonzeros of mat in all even numbers of row and col */
187: PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
188: PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
189: PetscCall(MatSetFromOptions(matB));
190: PetscCall(MatSetUp(matB));
191: PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
192: if (rstart % 2 != 0) rstart++;
193: for (i = rstart; i < rend; i += 2) {
194: for (j = 0; j < n; j += 2) {
195: v = 10.0 * i + j + 1.0;
196: PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
197: }
198: }
199: PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
200: PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
201: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
202: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
203: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
204: PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
205: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
206: PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
207: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
208: PetscCall(MatDestroy(&matB));
209: }
211: /* Test MatZeroRows */
212: j = rstart - 1;
213: if (j < 0) j = m - 1;
214: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
215: PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
216: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
218: /* Test MatShift */
219: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
220: PetscCall(MatShift(mat, -2.0));
221: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
223: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
224: /* Free data structures */
225: PetscCall(MatDestroy(&mat));
226: PetscCall(MatDestroy(&tmat));
227: PetscCall(PetscFinalize());
228: return 0;
229: }
231: /*TEST
233: test:
234: suffix: 11_A
235: args: -mat_type seqaij -rectA
236: filter: grep -v "Mat Object"
238: test:
239: suffix: 12_A
240: args: -mat_type seqdense -rectA
241: filter: grep -v type | grep -v "Mat Object"
243: test:
244: requires: cuda
245: suffix: 12_A_cuda
246: args: -mat_type seqdensecuda -rectA
247: output_file: output/ex2_12_A.out
248: filter: grep -v type | grep -v "Mat Object"
250: test:
251: requires: kokkos_kernels
252: suffix: 12_A_kokkos
253: args: -mat_type aijkokkos -rectA
254: output_file: output/ex2_12_A.out
255: filter: grep -v type | grep -v "Mat Object"
257: test:
258: suffix: 11_B
259: args: -mat_type seqaij -rectB
260: filter: grep -v "Mat Object"
262: test:
263: suffix: 12_B
264: args: -mat_type seqdense -rectB
265: filter: grep -v type | grep -v "Mat Object"
267: testset:
268: args: -rectB
269: output_file: output/ex2_12_B.out
270: filter: grep -v type | grep -v "Mat Object"
272: test:
273: requires: cuda
274: suffix: 12_B_cuda
275: args: -mat_type {{seqdensecuda seqaijcusparse}}
277: test:
278: requires: kokkos_kernels
279: suffix: 12_B_kokkos
280: args: -mat_type aijkokkos
282: test:
283: suffix: 12_B_aij
284: args: -mat_type aij
285: test:
286: suffix: 21
287: args: -mat_type mpiaij
288: filter: grep -v type | grep -v " MPI process"
290: test:
291: suffix: 22
292: args: -mat_type mpidense
293: filter: grep -v type | grep -v "Mat Object"
295: test:
296: requires: cuda
297: suffix: 22_cuda
298: output_file: output/ex2_22.out
299: args: -mat_type mpidensecuda
300: filter: grep -v type | grep -v "Mat Object"
302: test:
303: requires: kokkos_kernels
304: suffix: 22_kokkos
305: output_file: output/ex2_22.out
306: args: -mat_type aijkokkos
307: filter: grep -v type | grep -v "Mat Object"
309: test:
310: suffix: 23
311: nsize: 3
312: args: -mat_type mpiaij
313: filter: grep -v type | grep -v " MPI process"
315: test:
316: suffix: 24
317: nsize: 3
318: args: -mat_type mpidense
319: filter: grep -v type | grep -v "Mat Object"
321: test:
322: requires: cuda
323: suffix: 24_cuda
324: nsize: 3
325: output_file: output/ex2_24.out
326: args: -mat_type mpidensecuda
327: filter: grep -v type | grep -v "Mat Object"
329: test:
330: suffix: 2_aijcusparse_1
331: args: -mat_type mpiaijcusparse
332: output_file: output/ex2_21.out
333: requires: cuda
334: filter: grep -v type | grep -v " MPI process"
336: test:
337: suffix: 2_aijkokkos_1
338: args: -mat_type aijkokkos
339: output_file: output/ex2_21.out
340: requires: kokkos_kernels
341: filter: grep -v type | grep -v " MPI process"
343: test:
344: suffix: 2_aijcusparse_2
345: nsize: 3
346: args: -mat_type mpiaijcusparse
347: output_file: output/ex2_23.out
348: requires: cuda
349: filter: grep -v type | grep -v " MPI process"
351: test:
352: suffix: 2_aijkokkos_2
353: nsize: 3
354: args: -mat_type aijkokkos
355: output_file: output/ex2_23.out
356: # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
357: requires: !hip kokkos_kernels
358: filter: grep -v type | grep -v "MPI processes"
360: test:
361: suffix: 3
362: nsize: 2
363: args: -mat_type mpiaij -rectA
365: test:
366: suffix: 3_aijcusparse
367: nsize: 2
368: args: -mat_type mpiaijcusparse -rectA
369: requires: cuda
371: test:
372: suffix: 4
373: nsize: 2
374: args: -mat_type mpidense -rectA
375: filter: grep -v type | grep -v " MPI process"
377: test:
378: requires: cuda
379: suffix: 4_cuda
380: nsize: 2
381: output_file: output/ex2_4.out
382: args: -mat_type mpidensecuda -rectA
383: filter: grep -v type | grep -v " MPI process"
385: test:
386: suffix: aijcusparse_1
387: args: -mat_type seqaijcusparse -rectA
388: filter: grep -v "Mat Object"
389: output_file: output/ex2_11_A_aijcusparse.out
390: requires: cuda
392: TEST*/