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