xref: /petsc/src/mat/tests/ex2.c (revision 98d129c30f3ee9fdddc40fdbc5a989b7be64f888) 
1 static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
2 
3 #include <petscmat.h>
4 
5 static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
6 {
7   Mat     D, E, F, G;
8   MatType mtype;
9 
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));
61 
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 }
77 
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;
86 
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   }
104 
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));
119 
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   /* --------------- 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   }
136 
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));
144 
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));
152 
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   }
157 
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   }
169 
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   }
196 
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));
203 
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));
208 
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 }
216 
217 /*TEST
218 
219    test:
220       suffix: 11_A
221       args: -mat_type seqaij -rectA
222       filter: grep -v "Mat Object"
223 
224    test:
225       suffix: 12_A
226       args: -mat_type seqdense -rectA
227       filter: grep -v type | grep -v "Mat Object"
228 
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"
235 
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"
242 
243    test:
244       suffix: 11_B
245       args: -mat_type seqaij -rectB
246       filter: grep -v "Mat Object"
247 
248    test:
249       suffix: 12_B
250       args: -mat_type seqdense -rectB
251       filter: grep -v type | grep -v "Mat Object"
252 
253    testset:
254       args: -rectB
255       output_file: output/ex2_12_B.out
256       filter: grep -v type | grep -v "Mat Object"
257 
258       test:
259          requires: cuda
260          suffix: 12_B_cuda
261          args: -mat_type {{seqdensecuda seqaijcusparse}}
262 
263       test:
264          requires: kokkos_kernels
265          suffix: 12_B_kokkos
266          args: -mat_type aijkokkos
267 
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"
275 
276    test:
277       suffix: 22
278       args: -mat_type mpidense
279       filter: grep -v type | grep -v "Mat Object"
280 
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"
287 
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"
294 
295    test:
296       suffix: 23
297       nsize: 3
298       args: -mat_type mpiaij
299       filter: grep -v type | grep -v " MPI process"
300 
301    test:
302       suffix: 24
303       nsize: 3
304       args: -mat_type mpidense
305       filter: grep -v type | grep -v "Mat Object"
306 
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"
314 
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"
321 
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"
328 
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"
336 
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"
345 
346    test:
347       suffix: 3
348       nsize: 2
349       args: -mat_type mpiaij -rectA
350 
351    test:
352       suffix: 3_aijcusparse
353       nsize: 2
354       args: -mat_type mpiaijcusparse -rectA
355       requires: cuda
356 
357    test:
358       suffix: 4
359       nsize: 2
360       args: -mat_type mpidense -rectA
361       filter: grep -v type | grep -v " MPI process"
362 
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"
370 
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
377 
378 TEST*/
379