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(MatDestroy(&D));
142: PetscCall(MatDestroy(&C));
143: }
145: /* --------------- Test MatTranspose() -------------- */
146: PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
147: if (!rect && flg) {
148: PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
149: tmat = mat;
150: mat = NULL;
151: } else { /* out-of-place transpose */
152: PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
153: }
155: /* ----------------- Test MatNorm() ----------------- */
156: /* Print info about transpose matrix */
157: PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
158: PetscCall(MatNorm(tmat, NORM_1, &norm1));
159: PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
160: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenius norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
161: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
163: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
164: if (mat && !rect) {
165: alpha = 1.0;
166: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
167: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A\n"));
168: PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
169: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
171: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX: B = alpha*B + A\n"));
172: PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
173: PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
175: {
176: Mat A, B;
177: PetscBool equal;
179: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &A));
180: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &B));
181: PetscCall(MatAYPX(A, 2.0, A, SAME_NONZERO_PATTERN));
182: PetscCall(MatScale(B, 3.0));
183: PetscCall(MatEqual(A, B, &equal));
184: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatAYPX() failed when Y == X");
185: PetscCall(MatDestroy(&A));
186: PetscCall(MatDestroy(&B));
187: }
188: }
190: {
191: Mat C;
192: alpha = 1.0;
193: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
194: PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
195: PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
196: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
197: PetscCall(MatDestroy(&C));
198: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
199: PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
200: }
202: {
203: Mat matB;
204: /* get matB that has nonzeros of mat in all even numbers of row and col */
205: PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
206: PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
207: PetscCall(MatSetFromOptions(matB));
208: PetscCall(MatSetUp(matB));
209: PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
210: if (rstart % 2 != 0) rstart++;
211: for (i = rstart; i < rend; i += 2) {
212: for (j = 0; j < n; j += 2) {
213: v = 10.0 * i + j + 1.0;
214: PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
215: }
216: }
217: PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
218: PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
219: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
220: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
221: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
222: PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
223: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
224: PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
225: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
226: PetscCall(MatDestroy(&matB));
227: }
229: /* Test MatZeroRows */
230: j = rstart - 1;
231: if (j < 0) j = m - 1;
232: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
233: PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
234: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
236: /* Test MatShift */
237: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
238: PetscCall(MatShift(mat, -2.0));
239: PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
241: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
242: /* Free data structures */
243: PetscCall(MatDestroy(&mat));
244: PetscCall(MatDestroy(&tmat));
245: PetscCall(PetscFinalize());
246: return 0;
247: }
249: /*TEST
251: test:
252: suffix: 11_A
253: args: -mat_type seqaij -rectA
254: filter: grep -v "Mat Object"
256: test:
257: suffix: 12_A
258: args: -mat_type seqdense -rectA
259: filter: grep -v type | grep -v "Mat Object"
261: test:
262: requires: cuda
263: suffix: 12_A_cuda
264: args: -mat_type seqdensecuda -rectA
265: output_file: output/ex2_12_A.out
266: filter: grep -v type | grep -v "Mat Object"
268: test:
269: requires: kokkos_kernels
270: suffix: 12_A_kokkos
271: args: -mat_type aijkokkos -rectA
272: output_file: output/ex2_12_A.out
273: filter: grep -v type | grep -v "Mat Object"
275: test:
276: suffix: 11_B
277: args: -mat_type seqaij -rectB
278: filter: grep -v "Mat Object"
280: test:
281: suffix: 12_B
282: args: -mat_type seqdense -rectB
283: filter: grep -v type | grep -v "Mat Object"
285: testset:
286: args: -rectB
287: output_file: output/ex2_12_B.out
288: filter: grep -v type | grep -v "Mat Object"
290: test:
291: requires: cuda
292: suffix: 12_B_cuda
293: args: -mat_type {{seqdensecuda seqaijcusparse}}
295: test:
296: requires: kokkos_kernels
297: suffix: 12_B_kokkos
298: args: -mat_type aijkokkos
300: test:
301: suffix: 12_B_aij
302: args: -mat_type aij
303: test:
304: suffix: 21
305: args: -mat_type mpiaij
306: filter: grep -v type | grep -v " MPI process"
308: test:
309: suffix: 22
310: args: -mat_type mpidense
311: filter: grep -v type | grep -v "Mat Object"
313: test:
314: requires: cuda
315: suffix: 22_cuda
316: output_file: output/ex2_22.out
317: args: -mat_type mpidensecuda
318: filter: grep -v type | grep -v "Mat Object"
320: test:
321: requires: kokkos_kernels
322: suffix: 22_kokkos
323: output_file: output/ex2_22.out
324: args: -mat_type aijkokkos
325: filter: grep -v type | grep -v "Mat Object"
327: test:
328: suffix: 23
329: nsize: 3
330: args: -mat_type mpiaij
331: filter: grep -v type | grep -v " MPI process"
333: test:
334: suffix: 24
335: nsize: 3
336: args: -mat_type mpidense
337: filter: grep -v type | grep -v "Mat Object"
339: test:
340: requires: cuda
341: suffix: 24_cuda
342: nsize: 3
343: output_file: output/ex2_24.out
344: args: -mat_type mpidensecuda
345: filter: grep -v type | grep -v "Mat Object"
347: test:
348: suffix: 2_aijcusparse_1
349: args: -mat_type mpiaijcusparse
350: output_file: output/ex2_21.out
351: requires: cuda
352: filter: grep -v type | grep -v " MPI process"
354: test:
355: suffix: 2_aijkokkos_1
356: args: -mat_type aijkokkos
357: output_file: output/ex2_21.out
358: requires: kokkos_kernels
359: filter: grep -v type | grep -v " MPI process"
361: test:
362: suffix: 2_aijcusparse_2
363: nsize: 3
364: args: -mat_type mpiaijcusparse
365: output_file: output/ex2_23.out
366: requires: cuda
367: filter: grep -v type | grep -v " MPI process"
369: test:
370: suffix: 2_aijkokkos_2
371: nsize: 3
372: args: -mat_type aijkokkos
373: output_file: output/ex2_23.out
374: # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
375: requires: !hip kokkos_kernels
376: filter: grep -v type | grep -v "MPI processes"
378: test:
379: suffix: 3
380: nsize: 2
381: args: -mat_type mpiaij -rectA
383: test:
384: suffix: 3_aijkokkos
385: nsize: 2
386: args: -mat_type aijkokkos -rectA
387: output_file: output/ex2_3.out
388: requires: kokkos_kernels
389: filter: sed -e "s/mpiaijkokkos/mpiaij/"
391: test:
392: suffix: 3_aijcusparse
393: nsize: 2
394: args: -mat_type mpiaijcusparse -rectA
395: requires: cuda
397: test:
398: suffix: 4
399: nsize: 2
400: args: -mat_type mpidense -rectA
401: filter: grep -v type | grep -v " MPI process"
403: test:
404: requires: cuda
405: suffix: 4_cuda
406: nsize: 2
407: output_file: output/ex2_4.out
408: args: -mat_type mpidensecuda -rectA
409: filter: grep -v type | grep -v " MPI process"
411: test:
412: suffix: aijcusparse_1
413: args: -mat_type seqaijcusparse -rectA
414: filter: grep -v "Mat Object"
415: output_file: output/ex2_11_A_aijcusparse.out
416: requires: cuda
418: TEST*/