xref: /petsc/src/mat/tests/ex2.c (revision fbf9dbe564678ed6eff1806adbc4c4f01b9743f4) 
1 
2 static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
3 
4 #include <petscmat.h>
5 
6 static PetscErrorCode TransposeAXPY(Mat C, PetscScalar alpha, Mat mat, PetscErrorCode (*f)(Mat, Mat *))
7 {
8   Mat     D, E, F, G;
9   MatType mtype;
10 
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));
62 
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 }
78 
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;
87 
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   }
105 
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));
120 
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));
127 
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   }
137 
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));
145 
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));
153 
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   }
158 
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   }
170 
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   }
197 
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));
204 
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));
209 
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 }
217 
218 /*TEST
219 
220    test:
221       suffix: 11_A
222       args: -mat_type seqaij -rectA
223       filter: grep -v "Mat Object"
224 
225    test:
226       suffix: 12_A
227       args: -mat_type seqdense -rectA
228       filter: grep -v type | grep -v "Mat Object"
229 
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"
236 
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"
243 
244    test:
245       suffix: 11_B
246       args: -mat_type seqaij -rectB
247       filter: grep -v "Mat Object"
248 
249    test:
250       suffix: 12_B
251       args: -mat_type seqdense -rectB
252       filter: grep -v type | grep -v "Mat Object"
253 
254    testset:
255       args: -rectB
256       output_file: output/ex2_12_B.out
257       filter: grep -v type | grep -v "Mat Object"
258 
259       test:
260          requires: cuda
261          suffix: 12_B_cuda
262          args: -mat_type {{seqdensecuda seqaijcusparse}}
263 
264       test:
265          requires: kokkos_kernels
266          suffix: 12_B_kokkos
267          args: -mat_type aijkokkos
268 
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"
276 
277    test:
278       suffix: 22
279       args: -mat_type mpidense
280       filter: grep -v type | grep -v "Mat Object"
281 
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"
288 
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"
295 
296    test:
297       suffix: 23
298       nsize: 3
299       args: -mat_type mpiaij
300       filter: grep -v type | grep -v " MPI process"
301 
302    test:
303       suffix: 24
304       nsize: 3
305       args: -mat_type mpidense
306       filter: grep -v type | grep -v "Mat Object"
307 
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"
315 
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"
322 
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"
329 
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"
337 
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"
346 
347    test:
348       suffix: 3
349       nsize: 2
350       args: -mat_type mpiaij -rectA
351 
352    test:
353       suffix: 3_aijcusparse
354       nsize: 2
355       args: -mat_type mpiaijcusparse -rectA
356       requires: cuda
357 
358    test:
359       suffix: 4
360       nsize: 2
361       args: -mat_type mpidense -rectA
362       filter: grep -v type | grep -v " MPI process"
363 
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"
371 
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
378 
379 TEST*/
380