xref: /petsc/src/mat/tests/ex2.c (revision ebead697dbf761eb322f829370bbe90b3bd93fa3) 
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   PetscFunctionBeginUser;
89   PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
90   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON));
91   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
92   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
93   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
94   n    = m;
95   PetscCall(PetscOptionsHasName(NULL,NULL,"-rectA",&flg));
96   if (flg) {n += 2; rect = 1;}
97   PetscCall(PetscOptionsHasName(NULL,NULL,"-rectB",&flg));
98   if (flg) {n -= 2; rect = 1;}
99 
100   /* ------- Assemble matrix --------- */
101   PetscCall(MatCreate(PETSC_COMM_WORLD,&mat));
102   PetscCall(MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n));
103   PetscCall(MatSetFromOptions(mat));
104   PetscCall(MatSetUp(mat));
105   PetscCall(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       PetscCall(MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES));
110     }
111   }
112   PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY));
113   PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY));
114 
115   /* ----------------- Test MatNorm()  ----------------- */
116   PetscCall(MatNorm(mat,NORM_FROBENIUS,&normf));
117   PetscCall(MatNorm(mat,NORM_1,&norm1));
118   PetscCall(MatNorm(mat,NORM_INFINITY,&normi));
119   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi));
120   PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
121 
122   /* --------------- Test MatTranspose()  -------------- */
123   PetscCall(PetscOptionsHasName(NULL,NULL,"-in_place",&flg));
124   if (!rect && flg) {
125     PetscCall(MatTranspose(mat,MAT_REUSE_MATRIX,&mat));   /* in-place transpose */
126     tmat = mat;
127     mat  = NULL;
128   } else { /* out-of-place transpose */
129     PetscCall(MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat));
130   }
131 
132   /* ----------------- Test MatNorm()  ----------------- */
133   /* Print info about transpose matrix */
134   PetscCall(MatNorm(tmat,NORM_FROBENIUS,&normf));
135   PetscCall(MatNorm(tmat,NORM_1,&norm1));
136   PetscCall(MatNorm(tmat,NORM_INFINITY,&normi));
137   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi));
138   PetscCall(MatView(tmat,PETSC_VIEWER_STDOUT_WORLD));
139 
140   /* ----------------- Test MatAXPY(), MatAYPX()  ----------------- */
141   if (mat && !rect) {
142     alpha = 1.0;
143     PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
144     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A\n"));
145     PetscCall(MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN));
146     PetscCall(MatView(tmat,PETSC_VIEWER_STDOUT_WORLD));
147 
148     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatAYPX:  B = alpha*B + A\n"));
149     PetscCall(MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN));
150     PetscCall(MatView(tmat,PETSC_VIEWER_STDOUT_WORLD));
151   }
152 
153   {
154     Mat C;
155     alpha = 1.0;
156     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n"));
157     PetscCall(MatDuplicate(mat,MAT_COPY_VALUES,&C));
158     PetscCall(MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN));
159     PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD));
160     PetscCall(MatDestroy(&C));
161     PetscCall(TransposeAXPY(C,alpha,mat,MatCreateTranspose));
162     PetscCall(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     PetscCall(MatCreate(PETSC_COMM_WORLD,&matB));
169     PetscCall(MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n));
170     PetscCall(MatSetFromOptions(matB));
171     PetscCall(MatSetUp(matB));
172     PetscCall(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         PetscCall(MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES));
178       }
179     }
180     PetscCall(MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY));
181     PetscCall(MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY));
182     PetscCall(PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n"));
183     PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
184     PetscCall(PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n"));
185     PetscCall(MatView(matB,PETSC_VIEWER_STDOUT_WORLD));
186     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n"));
187     PetscCall(MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN));
188     PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
189     PetscCall(MatDestroy(&matB));
190   }
191 
192   /* Test MatZeroRows */
193   j = rstart - 1;
194   if (j < 0) j = m-1;
195   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatZeroRows:\n"));
196   PetscCall(MatZeroRows(mat,1,&j,0.0,NULL,NULL));
197   PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
198 
199   /* Test MatShift */
200   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatShift: B = B - 2*I\n"));
201   PetscCall(MatShift(mat,-2.0));
202   PetscCall(MatView(mat,PETSC_VIEWER_STDOUT_WORLD));
203 
204   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
205   /* Free data structures */
206   PetscCall(MatDestroy(&mat));
207   PetscCall(MatDestroy(&tmat));
208   PetscCall(PetscFinalize());
209   return 0;
210 }
211 
212 /*TEST
213 
214    test:
215       suffix: 11_A
216       args: -mat_type seqaij -rectA
217       filter: grep -v "Mat Object"
218 
219    test:
220       suffix: 12_A
221       args: -mat_type seqdense -rectA
222       filter: grep -v type | grep -v "Mat Object"
223 
224    test:
225       requires: cuda
226       suffix: 12_A_cuda
227       args: -mat_type seqdensecuda -rectA
228       output_file: output/ex2_12_A.out
229       filter: grep -v type | grep -v "Mat Object"
230 
231    test:
232       requires: kokkos_kernels
233       suffix: 12_A_kokkos
234       args: -mat_type aijkokkos -rectA
235       output_file: output/ex2_12_A.out
236       filter: grep -v type | grep -v "Mat Object"
237 
238    test:
239       suffix: 11_B
240       args: -mat_type seqaij -rectB
241       filter: grep -v "Mat Object"
242 
243    test:
244       suffix: 12_B
245       args: -mat_type seqdense -rectB
246       filter: grep -v type | grep -v "Mat Object"
247 
248    testset:
249       args: -rectB
250       output_file: output/ex2_12_B.out
251       filter: grep -v type | grep -v "Mat Object"
252 
253       test:
254          requires: cuda
255          suffix: 12_B_cuda
256          args: -mat_type {{seqdensecuda seqaijcusparse}}
257 
258       test:
259          requires: kokkos_kernels
260          suffix: 12_B_kokkos
261          args: -mat_type aijkokkos
262 
263       test:
264          suffix: 12_B_aij
265          args: -mat_type aij
266    test:
267       suffix: 21
268       args: -mat_type mpiaij
269       filter: grep -v type | grep -v " MPI process"
270 
271    test:
272       suffix: 22
273       args: -mat_type mpidense
274       filter: grep -v type | grep -v "Mat Object"
275 
276    test:
277       requires: cuda
278       suffix: 22_cuda
279       output_file: output/ex2_22.out
280       args: -mat_type mpidensecuda
281       filter: grep -v type | grep -v "Mat Object"
282 
283    test:
284       requires: kokkos_kernels
285       suffix: 22_kokkos
286       output_file: output/ex2_22.out
287       args: -mat_type aijkokkos
288       filter: grep -v type | grep -v "Mat Object"
289 
290    test:
291       suffix: 23
292       nsize: 3
293       args: -mat_type mpiaij
294       filter: grep -v type | grep -v " MPI process"
295 
296    test:
297       suffix: 24
298       nsize: 3
299       args: -mat_type mpidense
300       filter: grep -v type | grep -v "Mat Object"
301 
302    test:
303       requires: cuda
304       suffix: 24_cuda
305       nsize: 3
306       output_file: output/ex2_24.out
307       args: -mat_type mpidensecuda
308       filter: grep -v type | grep -v "Mat Object"
309 
310    test:
311       suffix: 2_aijcusparse_1
312       args: -mat_type mpiaijcusparse
313       output_file: output/ex2_21.out
314       requires: cuda
315       filter: grep -v type | grep -v " MPI process"
316 
317    test:
318       suffix: 2_aijkokkos_1
319       args: -mat_type aijkokkos
320       output_file: output/ex2_21.out
321       requires: kokkos_kernels
322       filter: grep -v type | grep -v " MPI process"
323 
324    test:
325       suffix: 2_aijcusparse_2
326       nsize: 3
327       args: -mat_type mpiaijcusparse
328       output_file: output/ex2_23.out
329       requires: cuda
330       filter: grep -v type | grep -v " MPI process"
331 
332    test:
333       suffix: 2_aijkokkos_2
334       nsize: 3
335       args: -mat_type aijkokkos
336       output_file: output/ex2_23.out
337       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
338       requires: !sycl !hip kokkos_kernels
339       filter: grep -v type | grep -v " MPI process"
340 
341    test:
342       suffix: 3
343       nsize: 2
344       args: -mat_type mpiaij -rectA
345 
346    test:
347       suffix: 3_aijcusparse
348       nsize: 2
349       args: -mat_type mpiaijcusparse -rectA
350       requires: cuda
351 
352    test:
353       suffix: 4
354       nsize: 2
355       args: -mat_type mpidense -rectA
356       filter: grep -v type | grep -v " MPI process"
357 
358    test:
359       requires: cuda
360       suffix: 4_cuda
361       nsize: 2
362       output_file: output/ex2_4.out
363       args: -mat_type mpidensecuda -rectA
364       filter: grep -v type | grep -v " MPI process"
365 
366    test:
367       suffix: aijcusparse_1
368       args: -mat_type seqaijcusparse -rectA
369       filter: grep -v "Mat Object"
370       output_file: output/ex2_11_A_aijcusparse.out
371       requires: cuda
372 
373 TEST*/
374