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