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