Actual source code: ex5.c
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;
11: PetscInt i,j,m = 8,n,rstart,rend,vstart,vend;
12: PetscScalar one = 1.0,negone = -1.0,v,alpha=0.1;
13: PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON;
14: PetscBool flg;
16: PetscInitialize(&argc,&args,(char*)0,help);
17: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
18: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
19: n = m;
20: PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
21: if (flg) n += 2;
22: PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
23: if (flg) n -= 2;
25: /* ---------- Assemble matrix and vectors ----------- */
27: MatCreate(PETSC_COMM_WORLD,&C);
28: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
29: MatSetFromOptions(C);
30: MatSetUp(C);
31: MatGetOwnershipRange(C,&rstart,&rend);
32: VecCreate(PETSC_COMM_WORLD,&x);
33: VecSetSizes(x,PETSC_DECIDE,m);
34: VecSetFromOptions(x);
35: VecDuplicate(x,&z);
36: VecDuplicate(x,&w);
37: VecCreate(PETSC_COMM_WORLD,&y);
38: VecSetSizes(y,PETSC_DECIDE,n);
39: VecSetFromOptions(y);
40: VecDuplicate(y,&u);
41: VecDuplicate(y,&s);
42: VecGetOwnershipRange(y,&vstart,&vend);
44: /* Assembly */
45: for (i=rstart; i<rend; i++) {
46: v = 100*(i+1);
47: VecSetValues(z,1,&i,&v,INSERT_VALUES);
48: for (j=0; j<n; j++) {
49: v = 10*(i+1)+j+1;
50: MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
51: }
52: }
54: /* Flush off proc Vec values and do more assembly */
55: VecAssemblyBegin(z);
56: for (i=vstart; i<vend; i++) {
57: v = one*((PetscReal)i);
58: VecSetValues(y,1,&i,&v,INSERT_VALUES);
59: v = 100.0*i;
60: VecSetValues(u,1,&i,&v,INSERT_VALUES);
61: }
63: /* Flush off proc Mat values and do more assembly */
64: MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);
65: for (i=rstart; i<rend; i++) {
66: for (j=0; j<n; j++) {
67: v = 10*(i+1)+j+1;
68: MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
69: }
70: }
71: /* Try overlap Coomunication with the next stage XXXSetValues */
72: 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);
84: PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");
85: CHKMEMQ;
86: MatMult(C,y,x);
87: CHKMEMQ;
88: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
89: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");
90: MatMultAdd(C,y,z,w);
91: VecAXPY(x,one,z);
92: VecAXPY(x,negone,w);
93: VecNorm(x,NORM_2,&norm);
94: if (norm > tol) {
95: PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
96: }
98: /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */
100: for (i=rstart; i<rend; i++) {
101: v = one*((PetscReal)i);
102: VecSetValues(x,1,&i,&v,INSERT_VALUES);
103: }
104: VecAssemblyBegin(x);
105: VecAssemblyEnd(x);
106: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");
107: MatMultTranspose(C,x,y);
108: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
110: PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
111: MatMultTransposeAdd(C,x,u,s);
112: VecAXPY(y,one,u);
113: VecAXPY(y,negone,s);
114: VecNorm(y,NORM_2,&norm);
115: if (norm > tol) {
116: PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
117: }
119: /* -------------------- Test MatGetDiagonal() ------------------ */
121: PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");
122: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
123: VecSet(x,one);
124: MatGetDiagonal(C,x);
125: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
126: for (i=vstart; i<vend; i++) {
127: v = one*((PetscReal)(i+1));
128: VecSetValues(y,1,&i,&v,INSERT_VALUES);
129: }
131: /* -------------------- Test () MatDiagonalScale ------------------ */
132: PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg);
133: if (flg) {
134: MatDiagonalScale(C,x,y);
135: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
136: }
137: /* Free data structures */
138: VecDestroy(&u)); PetscCall(VecDestroy(&s);
139: VecDestroy(&w)); PetscCall(VecDestroy(&x);
140: VecDestroy(&y)); PetscCall(VecDestroy(&z);
141: MatDestroy(&C);
143: PetscFinalize();
144: return 0;
145: }
147: /*TEST
149: test:
150: suffix: 11_A
151: args: -mat_type seqaij -rectA
152: filter: grep -v type
154: test:
155: args: -mat_type seqdense -rectA
156: suffix: 12_A
158: test:
159: args: -mat_type seqaij -rectB
160: suffix: 11_B
161: filter: grep -v type
163: test:
164: args: -mat_type seqdense -rectB
165: suffix: 12_B
167: test:
168: suffix: 21
169: args: -mat_type mpiaij
170: filter: grep -v type
172: test:
173: suffix: 22
174: args: -mat_type mpidense
176: test:
177: suffix: 23
178: nsize: 3
179: args: -mat_type mpiaij
180: filter: grep -v type
182: test:
183: suffix: 24
184: nsize: 3
185: args: -mat_type mpidense
187: test:
188: suffix: 2_aijcusparse_1
189: args: -mat_type mpiaijcusparse -vec_type cuda
190: filter: grep -v type
191: output_file: output/ex5_21.out
192: requires: cuda
194: test:
195: nsize: 3
196: suffix: 2_aijcusparse_2
197: filter: grep -v type
198: args: -mat_type mpiaijcusparse -vec_type cuda
199: args: -sf_type {{basic neighbor}}
200: output_file: output/ex5_23.out
201: requires: cuda
203: test:
204: nsize: 3
205: suffix: 2_aijcusparse_3
206: filter: grep -v type
207: args: -mat_type mpiaijcusparse -vec_type cuda
208: args: -sf_type {{basic neighbor}}
209: output_file: output/ex5_23.out
210: requires: cuda defined(PETSC_HAVE_MPI_GPU_AWARE)
212: test:
213: suffix: 31
214: args: -mat_type mpiaij -test_diagonalscale
215: filter: grep -v type
217: test:
218: suffix: 32
219: args: -mat_type mpibaij -test_diagonalscale
220: filter: grep -v Mat_
222: test:
223: suffix: 33
224: nsize: 3
225: args: -mat_type mpiaij -test_diagonalscale
226: filter: grep -v type
228: test:
229: suffix: 34
230: nsize: 3
231: args: -mat_type mpibaij -test_diagonalscale
232: filter: grep -v Mat_
234: test:
235: suffix: 3_aijcusparse_1
236: args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
237: filter: grep -v type
238: output_file: output/ex5_31.out
239: requires: cuda
241: test:
242: suffix: 3_aijcusparse_2
243: nsize: 3
244: args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
245: filter: grep -v type
246: output_file: output/ex5_33.out
247: requires: cuda
249: test:
250: suffix: 3_kokkos
251: nsize: 3
252: args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale
253: filter: grep -v type
254: output_file: output/ex5_33.out
255: requires: !sycl kokkos_kernels
257: test:
258: suffix: aijcusparse_1
259: args: -mat_type seqaijcusparse -vec_type cuda -rectA
260: filter: grep -v type
261: output_file: output/ex5_11_A.out
262: requires: cuda
264: test:
265: suffix: aijcusparse_2
266: args: -mat_type seqaijcusparse -vec_type cuda -rectB
267: filter: grep -v type
268: output_file: output/ex5_11_B.out
269: requires: cuda
271: test:
272: suffix: sell_1
273: args: -mat_type sell
274: output_file: output/ex5_41.out
276: test:
277: suffix: sell_2
278: nsize: 3
279: args: -mat_type sell
280: output_file: output/ex5_43.out
282: test:
283: suffix: sell_3
284: args: -mat_type sell -test_diagonalscale
285: output_file: output/ex5_51.out
287: test:
288: suffix: sell_4
289: nsize: 3
290: args: -mat_type sell -test_diagonalscale
291: output_file: output/ex5_53.out
293: TEST*/