xref: /petsc/src/mat/tests/ex5.c (revision d5c9c0c4eebc2f2a01a1bd0c86fca87e2acd2a03) 
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    test:
208       nsize: 3
209       suffix: 2_aijcusparse_3
210       filter: grep -v type
211       args: -mat_type mpiaijcusparse -vec_type cuda
212       args: -sf_type {{basic neighbor}}
213       output_file: output/ex5_23.out
214       requires: cuda define(PETSC_HAVE_MPI_GPU_AWARE)
215 
216    test:
217       suffix: 31
218       args: -mat_type mpiaij -test_diagonalscale
219       filter: grep -v type
220 
221    test:
222       suffix: 32
223       args: -mat_type mpibaij -test_diagonalscale
224       filter: grep -v Mat_
225 
226    test:
227       suffix: 33
228       nsize: 3
229       args: -mat_type mpiaij -test_diagonalscale
230       filter: grep -v type
231 
232    test:
233       suffix: 34
234       nsize: 3
235       args: -mat_type mpibaij -test_diagonalscale
236       filter: grep -v Mat_
237 
238    test:
239       suffix: 3_aijcusparse_1
240       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
241       filter: grep -v type
242       output_file: output/ex5_31.out
243       requires: cuda
244 
245    test:
246       suffix: 3_aijcusparse_2
247       nsize: 3
248       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
249       filter: grep -v type
250       output_file: output/ex5_33.out
251       requires: cuda
252 
253    test:
254       suffix: aijcusparse_1
255       args: -mat_type seqaijcusparse -vec_type cuda -rectA
256       filter: grep -v type
257       output_file: output/ex5_11_A.out
258       requires: cuda
259 
260    test:
261       suffix: aijcusparse_2
262       args: -mat_type seqaijcusparse -vec_type cuda -rectB
263       filter: grep -v type
264       output_file: output/ex5_11_B.out
265       requires: cuda
266 
267    test:
268       suffix: sell_1
269       args: -mat_type sell
270       output_file: output/ex5_41.out
271 
272    test:
273       suffix: sell_2
274       nsize: 3
275       args: -mat_type sell
276       output_file: output/ex5_43.out
277 
278    test:
279       suffix: sell_3
280       args: -mat_type sell -test_diagonalscale
281       output_file: output/ex5_51.out
282 
283    test:
284       suffix: sell_4
285       nsize: 3
286       args: -mat_type sell -test_diagonalscale
287       output_file: output/ex5_53.out
288 
289 TEST*/
290