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