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, (char *)0, 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: Frobenious 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: Frobenious 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));
156: }
158: {
159: Mat C;
160: alpha = 1.0;
161: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
162: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
163: PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
164: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
165: PetscCall(MatDestroy(&C));
166: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
167: PetscCall(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: PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
174: PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
175: PetscCall(MatSetFromOptions(matB));
176: PetscCall(MatSetUp(matB));
177: PetscCall(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: PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
183: }
184: }
185: PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
187: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
188: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
189: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
190: PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
191: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
192: PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
193: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
194: PetscCall(MatDestroy(&matB));
195: }
197: /* Test MatZeroRows */
198: j = rstart - 1;
199: if (j < 0) j = m - 1;
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
201: PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
202: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
204: /* Test MatShift */
205: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
206: PetscCall(MatShift(mat, -2.0));
207: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
209: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
210: /* Free data structures */
211: PetscCall(MatDestroy(&mat));
212: PetscCall(MatDestroy(&tmat));
213: PetscCall(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*/