Actual source code: ex2.c
2: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
4: #include <petscmat.h>
6: static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
7: {
8: Mat D, E, F, G;
9: MatType mtype;
11: PetscFunctionBegin;
12: PetscCall(MatGetType(mat, &mtype));
13: if (f == MatCreateTranspose) {
14: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
15: } else {
16: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nMatAXPY: (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
17: }
18: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
19: PetscCall(f(C, &D));
20: PetscCall(f(D, &E));
21: PetscCall(MatAXPY(E, alpha, mat, SAME_NONZERO_PATTERN));
22: PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
23: PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
24: PetscCall(MatDestroy(&E));
25: PetscCall(MatDestroy(&D));
26: PetscCall(MatDestroy(&C));
27: if (f == MatCreateTranspose) {
28: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
29: } else {
30: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
31: }
32: if (f == MatCreateTranspose) {
33: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &D));
34: } else {
35: PetscCall(MatHermitianTranspose(mat, MAT_INITIAL_MATRIX, &D));
36: }
37: PetscCall(f(D, &E));
38: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
39: PetscCall(MatAXPY(C, alpha, E, SAME_NONZERO_PATTERN));
40: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
41: PetscCall(MatDestroy(&E));
42: PetscCall(MatDestroy(&D));
43: PetscCall(MatDestroy(&C));
44: if (f == MatCreateTranspose) {
45: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n"));
46: } else {
47: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n"));
48: }
49: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
50: PetscCall(f(C, &D));
51: PetscCall(f(D, &E));
52: PetscCall(f(mat, &F));
53: PetscCall(f(F, &G));
54: PetscCall(MatAXPY(E, alpha, G, SAME_NONZERO_PATTERN));
55: PetscCall(MatConvert(E, mtype, MAT_INPLACE_MATRIX, &E));
56: PetscCall(MatView(E, PETSC_VIEWER_STDOUT_WORLD));
57: PetscCall(MatDestroy(&G));
58: PetscCall(MatDestroy(&F));
59: PetscCall(MatDestroy(&E));
60: PetscCall(MatDestroy(&D));
61: PetscCall(MatDestroy(&C));
63: /* Call f on a matrix that does not implement the transposition */
64: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: Now without the transposition operation\n"));
65: PetscCall(MatConvert(mat, MATSHELL, MAT_INITIAL_MATRIX, &C));
66: PetscCall(f(C, &D));
67: PetscCall(f(D, &E));
68: /* XXX cannot use MAT_INPLACE_MATRIX, it leaks mat */
69: PetscCall(MatConvert(E, mtype, MAT_INITIAL_MATRIX, &F));
70: PetscCall(MatAXPY(F, alpha, mat, SAME_NONZERO_PATTERN));
71: PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
72: PetscCall(MatDestroy(&F));
73: PetscCall(MatDestroy(&E));
74: PetscCall(MatDestroy(&D));
75: PetscCall(MatDestroy(&C));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: int main(int argc, char **argv)
80: {
81: Mat mat, tmat = 0;
82: PetscInt m = 7, n, i, j, rstart, rend, rect = 0;
83: PetscMPIInt size, rank;
84: PetscBool flg;
85: PetscScalar v, alpha;
86: PetscReal normf, normi, norm1;
88: PetscFunctionBeginUser;
89: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
90: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
91: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
92: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
93: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
94: n = m;
95: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
96: if (flg) {
97: n += 2;
98: rect = 1;
99: }
100: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
101: if (flg) {
102: n -= 2;
103: rect = 1;
104: }
106: /* ------- Assemble matrix --------- */
107: PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
108: PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
109: PetscCall(MatSetFromOptions(mat));
110: PetscCall(MatSetUp(mat));
111: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
112: for (i = rstart; i < rend; i++) {
113: for (j = 0; j < n; j++) {
114: v = 10.0 * i + j + 1.0;
115: PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
116: }
117: }
118: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
119: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
121: /* ----------------- Test MatNorm() ----------------- */
122: PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
123: PetscCall(MatNorm(mat, NORM_1, &norm1));
124: PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
125: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
126: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
128: /* --------------- Test MatTranspose() -------------- */
129: PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
130: if (!rect && flg) {
131: PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
132: tmat = mat;
133: mat = NULL;
134: } else { /* out-of-place transpose */
135: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
136: }
138: /* ----------------- Test MatNorm() ----------------- */
139: /* Print info about transpose matrix */
140: PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
141: PetscCall(MatNorm(tmat, NORM_1, &norm1));
142: PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
143: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
144: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
146: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
147: if (mat && !rect) {
148: alpha = 1.0;
149: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A\n"));
151: PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
152: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX: B = alpha*B + A\n"));
155: PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
156: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
157: }
159: {
160: Mat C;
161: alpha = 1.0;
162: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
163: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
164: PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
165: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
166: PetscCall(MatDestroy(&C));
167: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
168: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
169: }
171: {
172: Mat matB;
173: /* get matB that has nonzeros of mat in all even numbers of row and col */
174: PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
175: PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
176: PetscCall(MatSetFromOptions(matB));
177: PetscCall(MatSetUp(matB));
178: PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
179: if (rstart % 2 != 0) rstart++;
180: for (i = rstart; i < rend; i += 2) {
181: for (j = 0; j < n; j += 2) {
182: v = 10.0 * i + j + 1.0;
183: PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
184: }
185: }
186: PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
187: PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
189: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
190: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
191: PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
192: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
193: PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
194: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
195: PetscCall(MatDestroy(&matB));
196: }
198: /* Test MatZeroRows */
199: j = rstart - 1;
200: if (j < 0) j = m - 1;
201: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
202: PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
203: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
205: /* Test MatShift */
206: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
207: PetscCall(MatShift(mat, -2.0));
208: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
210: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
211: /* Free data structures */
212: PetscCall(MatDestroy(&mat));
213: PetscCall(MatDestroy(&tmat));
214: PetscCall(PetscFinalize());
215: return 0;
216: }
218: /*TEST
220: test:
221: suffix: 11_A
222: args: -mat_type seqaij -rectA
223: filter: grep -v "Mat Object"
225: test:
226: suffix: 12_A
227: args: -mat_type seqdense -rectA
228: filter: grep -v type | grep -v "Mat Object"
230: test:
231: requires: cuda
232: suffix: 12_A_cuda
233: args: -mat_type seqdensecuda -rectA
234: output_file: output/ex2_12_A.out
235: filter: grep -v type | grep -v "Mat Object"
237: test:
238: requires: kokkos_kernels
239: suffix: 12_A_kokkos
240: args: -mat_type aijkokkos -rectA
241: output_file: output/ex2_12_A.out
242: filter: grep -v type | grep -v "Mat Object"
244: test:
245: suffix: 11_B
246: args: -mat_type seqaij -rectB
247: filter: grep -v "Mat Object"
249: test:
250: suffix: 12_B
251: args: -mat_type seqdense -rectB
252: filter: grep -v type | grep -v "Mat Object"
254: testset:
255: args: -rectB
256: output_file: output/ex2_12_B.out
257: filter: grep -v type | grep -v "Mat Object"
259: test:
260: requires: cuda
261: suffix: 12_B_cuda
262: args: -mat_type {{seqdensecuda seqaijcusparse}}
264: test:
265: requires: kokkos_kernels
266: suffix: 12_B_kokkos
267: args: -mat_type aijkokkos
269: test:
270: suffix: 12_B_aij
271: args: -mat_type aij
272: test:
273: suffix: 21
274: args: -mat_type mpiaij
275: filter: grep -v type | grep -v " MPI process"
277: test:
278: suffix: 22
279: args: -mat_type mpidense
280: filter: grep -v type | grep -v "Mat Object"
282: test:
283: requires: cuda
284: suffix: 22_cuda
285: output_file: output/ex2_22.out
286: args: -mat_type mpidensecuda
287: filter: grep -v type | grep -v "Mat Object"
289: test:
290: requires: kokkos_kernels
291: suffix: 22_kokkos
292: output_file: output/ex2_22.out
293: args: -mat_type aijkokkos
294: filter: grep -v type | grep -v "Mat Object"
296: test:
297: suffix: 23
298: nsize: 3
299: args: -mat_type mpiaij
300: filter: grep -v type | grep -v " MPI process"
302: test:
303: suffix: 24
304: nsize: 3
305: args: -mat_type mpidense
306: filter: grep -v type | grep -v "Mat Object"
308: test:
309: requires: cuda
310: suffix: 24_cuda
311: nsize: 3
312: output_file: output/ex2_24.out
313: args: -mat_type mpidensecuda
314: filter: grep -v type | grep -v "Mat Object"
316: test:
317: suffix: 2_aijcusparse_1
318: args: -mat_type mpiaijcusparse
319: output_file: output/ex2_21.out
320: requires: cuda
321: filter: grep -v type | grep -v " MPI process"
323: test:
324: suffix: 2_aijkokkos_1
325: args: -mat_type aijkokkos
326: output_file: output/ex2_21.out
327: requires: kokkos_kernels
328: filter: grep -v type | grep -v " MPI process"
330: test:
331: suffix: 2_aijcusparse_2
332: nsize: 3
333: args: -mat_type mpiaijcusparse
334: output_file: output/ex2_23.out
335: requires: cuda
336: filter: grep -v type | grep -v " MPI process"
338: test:
339: suffix: 2_aijkokkos_2
340: nsize: 3
341: args: -mat_type aijkokkos
342: output_file: output/ex2_23.out
343: # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
344: requires: !hip kokkos_kernels
345: filter: grep -v type | grep -v "MPI processes"
347: test:
348: suffix: 3
349: nsize: 2
350: args: -mat_type mpiaij -rectA
352: test:
353: suffix: 3_aijcusparse
354: nsize: 2
355: args: -mat_type mpiaijcusparse -rectA
356: requires: cuda
358: test:
359: suffix: 4
360: nsize: 2
361: args: -mat_type mpidense -rectA
362: filter: grep -v type | grep -v " MPI process"
364: test:
365: requires: cuda
366: suffix: 4_cuda
367: nsize: 2
368: output_file: output/ex2_4.out
369: args: -mat_type mpidensecuda -rectA
370: filter: grep -v type | grep -v " MPI process"
372: test:
373: suffix: aijcusparse_1
374: args: -mat_type seqaijcusparse -rectA
375: filter: grep -v "Mat Object"
376: output_file: output/ex2_11_A_aijcusparse.out
377: requires: cuda
379: TEST*/