Actual source code: ex5.c
petsc-3.12.5 2020-03-29
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);
75: MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);
76: CHKMEMQ;
77: /* The Assembly for the second Stage */
78: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
80: VecAssemblyBegin(y);
81: VecAssemblyEnd(y);
82: MatScale(C,alpha);
83: VecAssemblyBegin(u);
84: VecAssemblyEnd(u);
85: PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");
86: CHKMEMQ;
87: MatMult(C,y,x);
88: CHKMEMQ;
89: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
90: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");
91: MatMultAdd(C,y,z,w);
92: VecAXPY(x,one,z);
93: VecAXPY(x,negone,w);
94: VecNorm(x,NORM_2,&norm);
95: if (norm > tol) {
96: PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
97: }
99: /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */
101: for (i=rstart; i<rend; i++) {
102: v = one*((PetscReal)i);
103: VecSetValues(x,1,&i,&v,INSERT_VALUES);
104: }
105: VecAssemblyBegin(x);
106: VecAssemblyEnd(x);
107: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");
108: MatMultTranspose(C,x,y);
109: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
111: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
112: MatMultTransposeAdd(C,x,u,s);
113: VecAXPY(y,one,u);
114: VecAXPY(y,negone,s);
115: VecNorm(y,NORM_2,&norm);
116: if (norm > tol) {
117: PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
118: }
120: /* -------------------- Test MatGetDiagonal() ------------------ */
122: PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");
123: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
124: VecSet(x,one);
125: MatGetDiagonal(C,x);
126: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
127: for (i=vstart; i<vend; i++) {
128: v = one*((PetscReal)(i+1));
129: VecSetValues(y,1,&i,&v,INSERT_VALUES);
130: }
132: /* -------------------- Test () MatDiagonalScale ------------------ */
133: PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg);
134: if (flg) {
135: MatDiagonalScale(C,x,y);
136: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
137: }
138: /* Free data structures */
139: VecDestroy(&u); VecDestroy(&s);
140: VecDestroy(&w); VecDestroy(&x);
141: VecDestroy(&y); VecDestroy(&z);
142: MatDestroy(&C);
144: PetscFinalize();
145: return ierr;
146: }
149: /*TEST
152: test:
153: suffix: 11_A
154: args: -mat_type seqaij -rectA
155: filter: grep -v type
157: test:
158: args: -mat_type seqdense -rectA
159: suffix: 12_A
161: test:
162: args: -mat_type seqaij -rectB
163: suffix: 11_B
164: filter: grep -v type
166: test:
167: args: -mat_type seqdense -rectB
168: suffix: 12_B
170: test:
171: suffix: 21
172: args: -mat_type mpiaij
173: filter: grep -v type
175: test:
176: suffix: 22
177: args: -mat_type mpidense
179: test:
180: suffix: 23
181: nsize: 3
182: args: -mat_type mpiaij
183: filter: grep -v type
185: test:
186: suffix: 24
187: nsize: 3
188: args: -mat_type mpidense
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
197: test:
198: suffix: 2_aijcusparse_2
199: nsize: 3
200: args: -mat_type mpiaijcusparse -vec_type cuda
201: filter: grep -v type
202: output_file: output/ex5_23.out
203: requires: cuda
205: test:
206: suffix: 2_aijcusparse_3
207: nsize: 3
208: args: -mat_type mpiaijcusparse -vec_type cuda -vecscatter_type sf
209: filter: grep -v type
210: output_file: output/ex5_23.out
211: requires: cuda
213: test:
214: suffix: 31
215: args: -mat_type mpiaij -test_diagonalscale
216: filter: grep -v type
218: test:
219: suffix: 32
220: args: -mat_type mpibaij -test_diagonalscale
221: filter: grep -v Mat_
223: test:
224: suffix: 33
225: nsize: 3
226: args: -mat_type mpiaij -test_diagonalscale
227: filter: grep -v type
229: test:
230: suffix: 34
231: nsize: 3
232: args: -mat_type mpibaij -test_diagonalscale
233: filter: grep -v Mat_
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
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
250: test:
251: suffix: aijcusparse_1
252: args: -mat_type seqaijcusparse -vec_type cuda -rectA
253: filter: grep -v type
254: output_file: output/ex5_11_A.out
255: requires: cuda
257: test:
258: suffix: aijcusparse_2
259: args: -mat_type seqaijcusparse -vec_type cuda -rectB
260: filter: grep -v type
261: output_file: output/ex5_11_B.out
262: requires: cuda
264: test:
265: suffix: sell_1
266: args: -mat_type sell
267: output_file: output/ex5_41.out
269: test:
270: suffix: sell_2
271: nsize: 3
272: args: -mat_type sell
273: output_file: output/ex5_43.out
275: test:
276: suffix: sell_3
277: args: -mat_type sell -test_diagonalscale
278: output_file: output/ex5_51.out
280: test:
281: suffix: sell_4
282: nsize: 3
283: args: -mat_type sell -test_diagonalscale
284: output_file: output/ex5_53.out
286: TEST*/