Actual source code: ex5.c
petsc-3.13.6 2020-09-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
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
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)
216: test:
217: suffix: 31
218: args: -mat_type mpiaij -test_diagonalscale
219: filter: grep -v type
221: test:
222: suffix: 32
223: args: -mat_type mpibaij -test_diagonalscale
224: filter: grep -v Mat_
226: test:
227: suffix: 33
228: nsize: 3
229: args: -mat_type mpiaij -test_diagonalscale
230: filter: grep -v type
232: test:
233: suffix: 34
234: nsize: 3
235: args: -mat_type mpibaij -test_diagonalscale
236: filter: grep -v Mat_
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
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
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
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
267: test:
268: suffix: sell_1
269: args: -mat_type sell
270: output_file: output/ex5_41.out
272: test:
273: suffix: sell_2
274: nsize: 3
275: args: -mat_type sell
276: output_file: output/ex5_43.out
278: test:
279: suffix: sell_3
280: args: -mat_type sell -test_diagonalscale
281: output_file: output/ex5_51.out
283: test:
284: suffix: sell_4
285: nsize: 3
286: args: -mat_type sell -test_diagonalscale
287: output_file: output/ex5_53.out
289: TEST*/