xref: /petsc/src/mat/tests/ex5.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
1 
2 static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\
3 Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), and MatDiagonalScale().\n\n";
4 
5 #include <petscmat.h>
6 
7 int main(int argc,char **args)
8 {
9   Mat            C;
10   Vec            s,u,w,x,y,z;
11   PetscErrorCode ierr;
12   PetscInt       i,j,m = 8,n,rstart,rend,vstart,vend;
13   PetscScalar    one = 1.0,negone = -1.0,v,alpha=0.1;
14   PetscReal      norm, tol = PETSC_SQRT_MACHINE_EPSILON;
15   PetscBool      flg;
16 
17   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
18   ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);CHKERRQ(ierr);
19   ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);
20   n    = m;
21   ierr = PetscOptionsHasName(NULL,NULL,"-rectA",&flg);CHKERRQ(ierr);
22   if (flg) n += 2;
23   ierr = PetscOptionsHasName(NULL,NULL,"-rectB",&flg);CHKERRQ(ierr);
24   if (flg) n -= 2;
25 
26   /* ---------- Assemble matrix and vectors ----------- */
27 
28   ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
29   ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
30   ierr = MatSetFromOptions(C);CHKERRQ(ierr);
31   ierr = MatSetUp(C);CHKERRQ(ierr);
32   ierr = MatGetOwnershipRange(C,&rstart,&rend);CHKERRQ(ierr);
33   ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
34   ierr = VecSetSizes(x,PETSC_DECIDE,m);CHKERRQ(ierr);
35   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
36   ierr = VecDuplicate(x,&z);CHKERRQ(ierr);
37   ierr = VecDuplicate(x,&w);CHKERRQ(ierr);
38   ierr = VecCreate(PETSC_COMM_WORLD,&y);CHKERRQ(ierr);
39   ierr = VecSetSizes(y,PETSC_DECIDE,n);CHKERRQ(ierr);
40   ierr = VecSetFromOptions(y);CHKERRQ(ierr);
41   ierr = VecDuplicate(y,&u);CHKERRQ(ierr);
42   ierr = VecDuplicate(y,&s);CHKERRQ(ierr);
43   ierr = VecGetOwnershipRange(y,&vstart,&vend);CHKERRQ(ierr);
44 
45   /* Assembly */
46   for (i=rstart; i<rend; i++) {
47     v    = 100*(i+1);
48     ierr = VecSetValues(z,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
49     for (j=0; j<n; j++) {
50       v    = 10*(i+1)+j+1;
51       ierr = MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);
52     }
53   }
54 
55   /* Flush off proc Vec values and do more assembly */
56   ierr = VecAssemblyBegin(z);CHKERRQ(ierr);
57   for (i=vstart; i<vend; i++) {
58     v    = one*((PetscReal)i);
59     ierr = VecSetValues(y,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
60     v    = 100.0*i;
61     ierr = VecSetValues(u,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
62   }
63 
64   /* Flush off proc Mat values and do more assembly */
65   ierr = MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
66   for (i=rstart; i<rend; i++) {
67     for (j=0; j<n; j++) {
68       v    = 10*(i+1)+j+1;
69       ierr = MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);
70     }
71   }
72   /* Try overlap Coomunication with the next stage XXXSetValues */
73   ierr = VecAssemblyEnd(z);CHKERRQ(ierr);
74 
75   ierr = MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
76   CHKMEMQ;
77   /* The Assembly for the second Stage */
78   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
79   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
80   ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
81   ierr = VecAssemblyEnd(y);CHKERRQ(ierr);
82   ierr = MatScale(C,alpha);CHKERRQ(ierr);
83   ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
84   ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
85   ierr = PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");CHKERRQ(ierr);
86   CHKMEMQ;
87   ierr = MatMult(C,y,x);CHKERRQ(ierr);
88   CHKMEMQ;
89   ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
90   ierr = PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");CHKERRQ(ierr);
91   ierr = MatMultAdd(C,y,z,w);CHKERRQ(ierr);
92   ierr = VecAXPY(x,one,z);CHKERRQ(ierr);
93   ierr = VecAXPY(x,negone,w);CHKERRQ(ierr);
94   ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
95   if (norm > tol) {
96     ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);CHKERRQ(ierr);
97   }
98 
99   /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */
100 
101   for (i=rstart; i<rend; i++) {
102     v    = one*((PetscReal)i);
103     ierr = VecSetValues(x,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
104   }
105   ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
106   ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
107   ierr = PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");CHKERRQ(ierr);
108   ierr = MatMultTranspose(C,x,y);CHKERRQ(ierr);
109   ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
110 
111   ierr = PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");CHKERRQ(ierr);
112   ierr = MatMultTransposeAdd(C,x,u,s);CHKERRQ(ierr);
113   ierr = VecAXPY(y,one,u);CHKERRQ(ierr);
114   ierr = VecAXPY(y,negone,s);CHKERRQ(ierr);
115   ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
116   if (norm > tol) {
117     ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);CHKERRQ(ierr);
118   }
119 
120   /* -------------------- Test MatGetDiagonal() ------------------ */
121 
122   ierr = PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");CHKERRQ(ierr);
123   ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
124   ierr = VecSet(x,one);CHKERRQ(ierr);
125   ierr = MatGetDiagonal(C,x);CHKERRQ(ierr);
126   ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
127   for (i=vstart; i<vend; i++) {
128     v    = one*((PetscReal)(i+1));
129     ierr = VecSetValues(y,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
130   }
131 
132   /* -------------------- Test () MatDiagonalScale ------------------ */
133   ierr = PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg);CHKERRQ(ierr);
134   if (flg) {
135     ierr = MatDiagonalScale(C,x,y);CHKERRQ(ierr);
136     ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
137   }
138   /* Free data structures */
139   ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&s);CHKERRQ(ierr);
140   ierr = VecDestroy(&w);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
141   ierr = VecDestroy(&y);CHKERRQ(ierr); ierr = VecDestroy(&z);CHKERRQ(ierr);
142   ierr = MatDestroy(&C);CHKERRQ(ierr);
143 
144   ierr = PetscFinalize();
145   return ierr;
146 }
147 
148 /*TEST
149 
150    test:
151       suffix: 11_A
152       args: -mat_type seqaij -rectA
153       filter: grep -v type
154 
155    test:
156       args: -mat_type seqdense -rectA
157       suffix: 12_A
158 
159    test:
160       args: -mat_type seqaij -rectB
161       suffix: 11_B
162       filter: grep -v type
163 
164    test:
165       args: -mat_type seqdense -rectB
166       suffix: 12_B
167 
168    test:
169       suffix: 21
170       args: -mat_type mpiaij
171       filter: grep -v type
172 
173    test:
174       suffix: 22
175       args: -mat_type mpidense
176 
177    test:
178       suffix: 23
179       nsize: 3
180       args: -mat_type mpiaij
181       filter: grep -v type
182 
183    test:
184       suffix: 24
185       nsize: 3
186       args: -mat_type mpidense
187 
188    test:
189       suffix: 2_aijcusparse_1
190       args: -mat_type mpiaijcusparse -vec_type cuda
191       filter: grep -v type
192       output_file: output/ex5_21.out
193       requires: cuda
194 
195    test:
196       nsize: 3
197       suffix: 2_aijcusparse_2
198       filter: grep -v type
199       args: -mat_type mpiaijcusparse -vec_type cuda
200       args: -sf_type {{basic neighbor}} -vecscatter_packongpu {{0 1}}
201       output_file: output/ex5_23.out
202       requires: cuda
203 
204    test:
205       nsize: 3
206       suffix: 2_aijcusparse_3
207       filter: grep -v type
208       args: -mat_type mpiaijcusparse -vec_type cuda
209       args: -sf_type {{basic neighbor}}
210       output_file: output/ex5_23.out
211       requires: cuda define(PETSC_HAVE_MPI_GPU_AWARE)
212 
213    test:
214       suffix: 31
215       args: -mat_type mpiaij -test_diagonalscale
216       filter: grep -v type
217 
218    test:
219       suffix: 32
220       args: -mat_type mpibaij -test_diagonalscale
221       filter: grep -v Mat_
222 
223    test:
224       suffix: 33
225       nsize: 3
226       args: -mat_type mpiaij -test_diagonalscale
227       filter: grep -v type
228 
229    test:
230       suffix: 34
231       nsize: 3
232       args: -mat_type mpibaij -test_diagonalscale
233       filter: grep -v Mat_
234 
235    test:
236       suffix: 3_aijcusparse_1
237       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
238       filter: grep -v type
239       output_file: output/ex5_31.out
240       requires: cuda
241 
242    test:
243       suffix: 3_aijcusparse_2
244       nsize: 3
245       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
246       filter: grep -v type
247       output_file: output/ex5_33.out
248       requires: cuda
249 
250    test:
251       suffix: 3_kokkos
252       nsize: 3
253       args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale
254       filter: grep -v type
255       output_file: output/ex5_33.out
256       requires: kokkos_kernels
257 
258    test:
259       suffix: aijcusparse_1
260       args: -mat_type seqaijcusparse -vec_type cuda -rectA
261       filter: grep -v type
262       output_file: output/ex5_11_A.out
263       requires: cuda
264 
265    test:
266       suffix: aijcusparse_2
267       args: -mat_type seqaijcusparse -vec_type cuda -rectB
268       filter: grep -v type
269       output_file: output/ex5_11_B.out
270       requires: cuda
271 
272    test:
273       suffix: sell_1
274       args: -mat_type sell
275       output_file: output/ex5_41.out
276 
277    test:
278       suffix: sell_2
279       nsize: 3
280       args: -mat_type sell
281       output_file: output/ex5_43.out
282 
283    test:
284       suffix: sell_3
285       args: -mat_type sell -test_diagonalscale
286       output_file: output/ex5_51.out
287 
288    test:
289       suffix: sell_4
290       nsize: 3
291       args: -mat_type sell -test_diagonalscale
292       output_file: output/ex5_53.out
293 
294 TEST*/
295