xref: /petsc/src/mat/tests/ex2.c (revision a69119a591a03a9d906b29c0a4e9802e4d7c9795) 
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   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(0);
76 }
77 
78 int main(int argc, char **argv) {
79   Mat         mat, tmat = 0;
80   PetscInt    m = 7, n, i, j, rstart, rend, rect = 0;
81   PetscMPIInt size, rank;
82   PetscBool   flg;
83   PetscScalar v, alpha;
84   PetscReal   normf, normi, norm1;
85 
86   PetscFunctionBeginUser;
87   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
88   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
89   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
90   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
91   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
92   n = m;
93   PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
94   if (flg) {
95     n += 2;
96     rect = 1;
97   }
98   PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
99   if (flg) {
100     n -= 2;
101     rect = 1;
102   }
103 
104   /* ------- Assemble matrix --------- */
105   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
106   PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
107   PetscCall(MatSetFromOptions(mat));
108   PetscCall(MatSetUp(mat));
109   PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
110   for (i = rstart; i < rend; i++) {
111     for (j = 0; j < n; j++) {
112       v = 10.0 * i + j + 1.0;
113       PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
114     }
115   }
116   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
117   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
118 
119   /* ----------------- Test MatNorm()  ----------------- */
120   PetscCall(MatNorm(mat, NORM_FROBENIUS, &normf));
121   PetscCall(MatNorm(mat, NORM_1, &norm1));
122   PetscCall(MatNorm(mat, NORM_INFINITY, &normi));
123   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
124   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
125 
126   /* --------------- Test MatTranspose()  -------------- */
127   PetscCall(PetscOptionsHasName(NULL, NULL, "-in_place", &flg));
128   if (!rect && flg) {
129     PetscCall(MatTranspose(mat, MAT_REUSE_MATRIX, &mat)); /* in-place transpose */
130     tmat = mat;
131     mat  = NULL;
132   } else { /* out-of-place transpose */
133     PetscCall(MatTranspose(mat, MAT_INITIAL_MATRIX, &tmat));
134   }
135 
136   /* ----------------- Test MatNorm()  ----------------- */
137   /* Print info about transpose matrix */
138   PetscCall(MatNorm(tmat, NORM_FROBENIUS, &normf));
139   PetscCall(MatNorm(tmat, NORM_1, &norm1));
140   PetscCall(MatNorm(tmat, NORM_INFINITY, &normi));
141   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n", (double)normf, (double)norm1, (double)normi));
142   PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
143 
144   /* ----------------- Test MatAXPY(), MatAYPX()  ----------------- */
145   if (mat && !rect) {
146     alpha = 1.0;
147     PetscCall(PetscOptionsGetScalar(NULL, NULL, "-alpha", &alpha, NULL));
148     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  B = B + alpha * A\n"));
149     PetscCall(MatAXPY(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
150     PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
151 
152     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAYPX:  B = alpha*B + A\n"));
153     PetscCall(MatAYPX(tmat, alpha, mat, DIFFERENT_NONZERO_PATTERN));
154     PetscCall(MatView(tmat, PETSC_VIEWER_STDOUT_WORLD));
155   }
156 
157   {
158     Mat C;
159     alpha = 1.0;
160     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
161     PetscCall(MatDuplicate(mat, MAT_COPY_VALUES, &C));
162     PetscCall(MatAXPY(C, alpha, mat, SAME_NONZERO_PATTERN));
163     PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
164     PetscCall(MatDestroy(&C));
165     PetscCall(TransposeAXPY(C, alpha, mat, MatCreateTranspose));
166     PetscCall(TransposeAXPY(C, alpha, mat, MatCreateHermitianTranspose));
167   }
168 
169   {
170     Mat matB;
171     /* get matB that has nonzeros of mat in all even numbers of row and col */
172     PetscCall(MatCreate(PETSC_COMM_WORLD, &matB));
173     PetscCall(MatSetSizes(matB, PETSC_DECIDE, PETSC_DECIDE, m, n));
174     PetscCall(MatSetFromOptions(matB));
175     PetscCall(MatSetUp(matB));
176     PetscCall(MatGetOwnershipRange(matB, &rstart, &rend));
177     if (rstart % 2 != 0) rstart++;
178     for (i = rstart; i < rend; i += 2) {
179       for (j = 0; j < n; j += 2) {
180         v = 10.0 * i + j + 1.0;
181         PetscCall(MatSetValues(matB, 1, &i, 1, &j, &v, INSERT_VALUES));
182       }
183     }
184     PetscCall(MatAssemblyBegin(matB, MAT_FINAL_ASSEMBLY));
185     PetscCall(MatAssemblyEnd(matB, MAT_FINAL_ASSEMBLY));
186     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A: original matrix:\n"));
187     PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
188     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " B(a subset of A):\n"));
189     PetscCall(MatView(matB, PETSC_VIEWER_STDOUT_WORLD));
190     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
191     PetscCall(MatAXPY(mat, alpha, matB, SUBSET_NONZERO_PATTERN));
192     PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
193     PetscCall(MatDestroy(&matB));
194   }
195 
196   /* Test MatZeroRows */
197   j = rstart - 1;
198   if (j < 0) j = m - 1;
199   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatZeroRows:\n"));
200   PetscCall(MatZeroRows(mat, 1, &j, 0.0, NULL, NULL));
201   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
202 
203   /* Test MatShift */
204   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift: B = B - 2*I\n"));
205   PetscCall(MatShift(mat, -2.0));
206   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
207 
208   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
209   /* Free data structures */
210   PetscCall(MatDestroy(&mat));
211   PetscCall(MatDestroy(&tmat));
212   PetscCall(PetscFinalize());
213   return 0;
214 }
215 
216 /*TEST
217 
218    test:
219       suffix: 11_A
220       args: -mat_type seqaij -rectA
221       filter: grep -v "Mat Object"
222 
223    test:
224       suffix: 12_A
225       args: -mat_type seqdense -rectA
226       filter: grep -v type | grep -v "Mat Object"
227 
228    test:
229       requires: cuda
230       suffix: 12_A_cuda
231       args: -mat_type seqdensecuda -rectA
232       output_file: output/ex2_12_A.out
233       filter: grep -v type | grep -v "Mat Object"
234 
235    test:
236       requires: kokkos_kernels
237       suffix: 12_A_kokkos
238       args: -mat_type aijkokkos -rectA
239       output_file: output/ex2_12_A.out
240       filter: grep -v type | grep -v "Mat Object"
241 
242    test:
243       suffix: 11_B
244       args: -mat_type seqaij -rectB
245       filter: grep -v "Mat Object"
246 
247    test:
248       suffix: 12_B
249       args: -mat_type seqdense -rectB
250       filter: grep -v type | grep -v "Mat Object"
251 
252    testset:
253       args: -rectB
254       output_file: output/ex2_12_B.out
255       filter: grep -v type | grep -v "Mat Object"
256 
257       test:
258          requires: cuda
259          suffix: 12_B_cuda
260          args: -mat_type {{seqdensecuda seqaijcusparse}}
261 
262       test:
263          requires: kokkos_kernels
264          suffix: 12_B_kokkos
265          args: -mat_type aijkokkos
266 
267       test:
268          suffix: 12_B_aij
269          args: -mat_type aij
270    test:
271       suffix: 21
272       args: -mat_type mpiaij
273       filter: grep -v type | grep -v " MPI process"
274 
275    test:
276       suffix: 22
277       args: -mat_type mpidense
278       filter: grep -v type | grep -v "Mat Object"
279 
280    test:
281       requires: cuda
282       suffix: 22_cuda
283       output_file: output/ex2_22.out
284       args: -mat_type mpidensecuda
285       filter: grep -v type | grep -v "Mat Object"
286 
287    test:
288       requires: kokkos_kernels
289       suffix: 22_kokkos
290       output_file: output/ex2_22.out
291       args: -mat_type aijkokkos
292       filter: grep -v type | grep -v "Mat Object"
293 
294    test:
295       suffix: 23
296       nsize: 3
297       args: -mat_type mpiaij
298       filter: grep -v type | grep -v " MPI process"
299 
300    test:
301       suffix: 24
302       nsize: 3
303       args: -mat_type mpidense
304       filter: grep -v type | grep -v "Mat Object"
305 
306    test:
307       requires: cuda
308       suffix: 24_cuda
309       nsize: 3
310       output_file: output/ex2_24.out
311       args: -mat_type mpidensecuda
312       filter: grep -v type | grep -v "Mat Object"
313 
314    test:
315       suffix: 2_aijcusparse_1
316       args: -mat_type mpiaijcusparse
317       output_file: output/ex2_21.out
318       requires: cuda
319       filter: grep -v type | grep -v " MPI process"
320 
321    test:
322       suffix: 2_aijkokkos_1
323       args: -mat_type aijkokkos
324       output_file: output/ex2_21.out
325       requires: kokkos_kernels
326       filter: grep -v type | grep -v " MPI process"
327 
328    test:
329       suffix: 2_aijcusparse_2
330       nsize: 3
331       args: -mat_type mpiaijcusparse
332       output_file: output/ex2_23.out
333       requires: cuda
334       filter: grep -v type | grep -v " MPI process"
335 
336    test:
337       suffix: 2_aijkokkos_2
338       nsize: 3
339       args: -mat_type aijkokkos
340       output_file: output/ex2_23.out
341       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
342       requires: !hip kokkos_kernels
343       filter: grep -v type | grep -v "MPI processes"
344 
345    test:
346       suffix: 3
347       nsize: 2
348       args: -mat_type mpiaij -rectA
349 
350    test:
351       suffix: 3_aijcusparse
352       nsize: 2
353       args: -mat_type mpiaijcusparse -rectA
354       requires: cuda
355 
356    test:
357       suffix: 4
358       nsize: 2
359       args: -mat_type mpidense -rectA
360       filter: grep -v type | grep -v " MPI process"
361 
362    test:
363       requires: cuda
364       suffix: 4_cuda
365       nsize: 2
366       output_file: output/ex2_4.out
367       args: -mat_type mpidensecuda -rectA
368       filter: grep -v type | grep -v " MPI process"
369 
370    test:
371       suffix: aijcusparse_1
372       args: -mat_type seqaijcusparse -rectA
373       filter: grep -v "Mat Object"
374       output_file: output/ex2_11_A_aijcusparse.out
375       requires: cuda
376 
377 TEST*/
378