xref: /petsc/src/mat/tests/ex5.c (revision 909dfd52d29915732d44f9443658cd10cc8ebe94) 
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 
149 /*TEST
150 
151 
152    test:
153       suffix: 11_A
154       args: -mat_type seqaij -rectA
155       filter: grep -v type
156 
157    test:
158       args: -mat_type seqdense -rectA
159       suffix: 12_A
160 
161    test:
162       args: -mat_type seqaij -rectB
163       suffix: 11_B
164       filter: grep -v type
165 
166    test:
167       args: -mat_type seqdense -rectB
168       suffix: 12_B
169 
170    test:
171       suffix: 21
172       args: -mat_type mpiaij
173       filter: grep -v type
174 
175    test:
176       suffix: 22
177       args: -mat_type mpidense
178 
179    test:
180       suffix: 23
181       nsize: 3
182       args: -mat_type mpiaij
183       filter: grep -v type
184 
185    test:
186       suffix: 24
187       nsize: 3
188       args: -mat_type mpidense
189 
190    test:
191       suffix: 2_aijcusparse_1
192       args: -mat_type mpiaijcusparse -vec_type cuda
193       filter: grep -v type
194       output_file: output/ex5_21.out
195       requires: cuda
196 
197 
198    test:
199       nsize: 3
200       suffix: 2_aijcusparse_2
201       filter: grep -v type
202       args: -mat_type mpiaijcusparse -vec_type cuda
203       args: -sf_type {{basic neighbor}} -vecscatter_packongpu {{0 1}}
204       output_file: output/ex5_23.out
205       requires: cuda
206 
207 
208    test:
209       nsize: 3
210       suffix: 2_aijcusparse_3
211       filter: grep -v type
212       args: -mat_type mpiaijcusparse -vec_type cuda
213       args: -sf_type {{basic neighbor}}
214       output_file: output/ex5_23.out
215       requires: cuda define(PETSC_HAVE_MPI_GPU_AWARE)
216 
217    test:
218       suffix: 31
219       args: -mat_type mpiaij -test_diagonalscale
220       filter: grep -v type
221 
222    test:
223       suffix: 32
224       args: -mat_type mpibaij -test_diagonalscale
225       filter: grep -v Mat_
226 
227    test:
228       suffix: 33
229       nsize: 3
230       args: -mat_type mpiaij -test_diagonalscale
231       filter: grep -v type
232 
233    test:
234       suffix: 34
235       nsize: 3
236       args: -mat_type mpibaij -test_diagonalscale
237       filter: grep -v Mat_
238 
239    test:
240       suffix: 3_aijcusparse_1
241       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
242       filter: grep -v type
243       output_file: output/ex5_31.out
244       requires: cuda
245 
246    test:
247       suffix: 3_aijcusparse_2
248       nsize: 3
249       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
250       filter: grep -v type
251       output_file: output/ex5_33.out
252       requires: cuda
253 
254    test:
255       suffix: 3_kokkos
256       nsize: 3
257       args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale
258       filter: grep -v type
259       output_file: output/ex5_33.out
260       requires: kokkos_kernels
261 
262    test:
263       suffix: aijcusparse_1
264       args: -mat_type seqaijcusparse -vec_type cuda -rectA
265       filter: grep -v type
266       output_file: output/ex5_11_A.out
267       requires: cuda
268 
269    test:
270       suffix: aijcusparse_2
271       args: -mat_type seqaijcusparse -vec_type cuda -rectB
272       filter: grep -v type
273       output_file: output/ex5_11_B.out
274       requires: cuda
275 
276    test:
277       suffix: sell_1
278       args: -mat_type sell
279       output_file: output/ex5_41.out
280 
281    test:
282       suffix: sell_2
283       nsize: 3
284       args: -mat_type sell
285       output_file: output/ex5_43.out
286 
287    test:
288       suffix: sell_3
289       args: -mat_type sell -test_diagonalscale
290       output_file: output/ex5_51.out
291 
292    test:
293       suffix: sell_4
294       nsize: 3
295       args: -mat_type sell -test_diagonalscale
296       output_file: output/ex5_53.out
297 
298 TEST*/
299