xref: /petsc/src/mat/tests/ex2.c (revision 2f605a994c568b341eedf7d10867052a20765ec2) 
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    test:
248       requires: cuda
249       suffix: 12_B_cuda
250       args: -mat_type seqdensecuda -rectB
251       output_file: output/ex2_12_B.out
252       filter: grep -v type | grep -v "Mat Object"
253 
254    test:
255       requires: kokkos_kernels
256       suffix: 12_B_kokkos
257       args: -mat_type aijkokkos -rectB
258       output_file: output/ex2_12_B.out
259       filter: grep -v type | grep -v "Mat Object"
260 
261    test:
262       suffix: 21
263       args: -mat_type mpiaij
264       filter: grep -v type | grep -v " MPI process"
265 
266    test:
267       suffix: 22
268       args: -mat_type mpidense
269       filter: grep -v type | grep -v "Mat Object"
270 
271    test:
272       requires: cuda
273       suffix: 22_cuda
274       output_file: output/ex2_22.out
275       args: -mat_type mpidensecuda
276       filter: grep -v type | grep -v "Mat Object"
277 
278    test:
279       requires: kokkos_kernels
280       suffix: 22_kokkos
281       output_file: output/ex2_22.out
282       args: -mat_type aijkokkos
283       filter: grep -v type | grep -v "Mat Object"
284 
285    test:
286       suffix: 23
287       nsize: 3
288       args: -mat_type mpiaij
289       filter: grep -v type | grep -v " MPI process"
290 
291    test:
292       suffix: 24
293       nsize: 3
294       args: -mat_type mpidense
295       filter: grep -v type | grep -v "Mat Object"
296 
297    test:
298       requires: cuda
299       suffix: 24_cuda
300       nsize: 3
301       output_file: output/ex2_24.out
302       args: -mat_type mpidensecuda
303       filter: grep -v type | grep -v "Mat Object"
304 
305    test:
306       suffix: 2_aijcusparse_1
307       args: -mat_type mpiaijcusparse
308       output_file: output/ex2_21.out
309       requires: cuda
310       filter: grep -v type | grep -v " MPI process"
311 
312    test:
313       suffix: 2_aijkokkos_1
314       args: -mat_type aijkokkos
315       output_file: output/ex2_21.out
316       requires: kokkos_kernels
317       filter: grep -v type | grep -v " MPI process"
318 
319    test:
320       suffix: 2_aijcusparse_2
321       nsize: 3
322       args: -mat_type mpiaijcusparse
323       output_file: output/ex2_23.out
324       requires: cuda
325       filter: grep -v type | grep -v " MPI process"
326 
327    test:
328       suffix: 2_aijkokkos_2
329       nsize: 3
330       args: -mat_type aijkokkos
331       output_file: output/ex2_23.out
332       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
333       requires: !sycl !hip kokkos_kernels
334       filter: grep -v type | grep -v " MPI process"
335 
336    test:
337       suffix: 3
338       nsize: 2
339       args: -mat_type mpiaij -rectA
340 
341    test:
342       suffix: 3_aijcusparse
343       nsize: 2
344       args: -mat_type mpiaijcusparse -rectA
345       requires: cuda
346 
347    test:
348       suffix: 4
349       nsize: 2
350       args: -mat_type mpidense -rectA
351       filter: grep -v type | grep -v " MPI process"
352 
353    test:
354       requires: cuda
355       suffix: 4_cuda
356       nsize: 2
357       output_file: output/ex2_4.out
358       args: -mat_type mpidensecuda -rectA
359       filter: grep -v type | grep -v " MPI process"
360 
361    test:
362       suffix: aijcusparse_1
363       args: -mat_type seqaijcusparse -rectA
364       filter: grep -v "Mat Object"
365       output_file: output/ex2_11_A_aijcusparse.out
366       requires: cuda
367 
368    test:
369       suffix: aijcusparse_2
370       args: -mat_type seqaijcusparse -rectB
371       filter: grep -v "Mat Object"
372       output_file: output/ex2_11_B_aijcusparse.out
373       requires: cuda
374 
375 TEST*/
376