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