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