xref: /petsc/src/mat/tests/ex2.c (revision ffa8c5705e8ab2cf85ee1d14dbe507a6e2eb5283) 
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(MatZeroRows(mat,1,&j,0.0,NULL,NULL));
195   PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
196 
197   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
198   /* Free data structures */
199   PetscCall(MatDestroy(&mat));
200   PetscCall(MatDestroy(&tmat));
201   PetscCall(PetscFinalize());
202   return 0;
203 }
204 
205 /*TEST
206 
207    test:
208       suffix: 11_A
209       args: -mat_type seqaij -rectA
210       filter: grep -v "Mat Object"
211 
212    test:
213       suffix: 12_A
214       args: -mat_type seqdense -rectA
215       filter: grep -v type | grep -v "Mat Object"
216 
217    test:
218       requires: cuda
219       suffix: 12_A_cuda
220       args: -mat_type seqdensecuda -rectA
221       output_file: output/ex2_12_A.out
222       filter: grep -v type | grep -v "Mat Object"
223 
224    test:
225       requires: kokkos_kernels
226       suffix: 12_A_kokkos
227       args: -mat_type aijkokkos -rectA
228       output_file: output/ex2_12_A.out
229       filter: grep -v type | grep -v "Mat Object"
230 
231    test:
232       suffix: 11_B
233       args: -mat_type seqaij -rectB
234       filter: grep -v "Mat Object"
235 
236    test:
237       suffix: 12_B
238       args: -mat_type seqdense -rectB
239       filter: grep -v type | grep -v "Mat Object"
240 
241    test:
242       requires: cuda
243       suffix: 12_B_cuda
244       args: -mat_type seqdensecuda -rectB
245       output_file: output/ex2_12_B.out
246       filter: grep -v type | grep -v "Mat Object"
247 
248    test:
249       requires: kokkos_kernels
250       suffix: 12_B_kokkos
251       args: -mat_type aijkokkos -rectB
252       output_file: output/ex2_12_B.out
253       filter: grep -v type | grep -v "Mat Object"
254 
255    test:
256       suffix: 21
257       args: -mat_type mpiaij
258       filter: grep -v type | grep -v "MPI processes"
259 
260    test:
261       suffix: 22
262       args: -mat_type mpidense
263       filter: grep -v type | grep -v "Mat Object"
264 
265    test:
266       requires: cuda
267       suffix: 22_cuda
268       output_file: output/ex2_22.out
269       args: -mat_type mpidensecuda
270       filter: grep -v type | grep -v "Mat Object"
271 
272    test:
273       requires: kokkos_kernels
274       suffix: 22_kokkos
275       output_file: output/ex2_22.out
276       args: -mat_type aijkokkos
277       filter: grep -v type | grep -v "Mat Object"
278 
279    test:
280       suffix: 23
281       nsize: 3
282       args: -mat_type mpiaij
283       filter: grep -v type | grep -v "MPI processes"
284 
285    test:
286       suffix: 24
287       nsize: 3
288       args: -mat_type mpidense
289       filter: grep -v type | grep -v "Mat Object"
290 
291    test:
292       requires: cuda
293       suffix: 24_cuda
294       nsize: 3
295       output_file: output/ex2_24.out
296       args: -mat_type mpidensecuda
297       filter: grep -v type | grep -v "Mat Object"
298 
299    test:
300       suffix: 2_aijcusparse_1
301       args: -mat_type mpiaijcusparse
302       output_file: output/ex2_21.out
303       requires: cuda
304       filter: grep -v type | grep -v "MPI processes"
305 
306    test:
307       suffix: 2_aijkokkos_1
308       args: -mat_type aijkokkos
309       output_file: output/ex2_21.out
310       requires: kokkos_kernels
311       filter: grep -v type | grep -v "MPI processes"
312 
313    test:
314       suffix: 2_aijcusparse_2
315       nsize: 3
316       args: -mat_type mpiaijcusparse
317       output_file: output/ex2_23.out
318       requires: cuda
319       filter: grep -v type | grep -v "MPI processes"
320 
321    test:
322       suffix: 2_aijkokkos_2
323       nsize: 3
324       args: -mat_type aijkokkos
325       output_file: output/ex2_23.out
326       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
327       requires: !sycl !hip kokkos_kernels
328       filter: grep -v type | grep -v "MPI processes"
329 
330    test:
331       suffix: 3
332       nsize: 2
333       args: -mat_type mpiaij -rectA
334 
335    test:
336       suffix: 3_aijcusparse
337       nsize: 2
338       args: -mat_type mpiaijcusparse -rectA
339       requires: cuda
340 
341    test:
342       suffix: 4
343       nsize: 2
344       args: -mat_type mpidense -rectA
345       filter: grep -v type | grep -v "MPI processes"
346 
347    test:
348       requires: cuda
349       suffix: 4_cuda
350       nsize: 2
351       output_file: output/ex2_4.out
352       args: -mat_type mpidensecuda -rectA
353       filter: grep -v type | grep -v "MPI processes"
354 
355    test:
356       suffix: aijcusparse_1
357       args: -mat_type seqaijcusparse -rectA
358       filter: grep -v "Mat Object"
359       output_file: output/ex2_11_A_aijcusparse.out
360       requires: cuda
361 
362    test:
363       suffix: aijcusparse_2
364       args: -mat_type seqaijcusparse -rectB
365       filter: grep -v "Mat Object"
366       output_file: output/ex2_11_B_aijcusparse.out
367       requires: cuda
368 
369 TEST*/
370