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;
17: PetscInitialize(&argc, &args, (char *)0, help);
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) PetscPrintf(PETSC_COMM_WORLD, "Norm of error difference = %g\n", (double)norm);
97: /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */
99: for (i = rstart; i < rend; i++) {
100: v = one * ((PetscReal)i);
101: VecSetValues(x, 1, &i, &v, INSERT_VALUES);
102: }
103: VecAssemblyBegin(x);
104: VecAssemblyEnd(x);
105: PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTranspose()\n");
106: MatMultTranspose(C, x, y);
107: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
109: PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTransposeAdd()\n");
110: MatMultTransposeAdd(C, x, u, s);
111: VecAXPY(y, one, u);
112: VecAXPY(y, negone, s);
113: VecNorm(y, NORM_2, &norm);
114: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "Norm of error difference = %g\n", (double)norm);
116: /* -------------------- Test MatGetDiagonal() ------------------ */
118: PetscPrintf(PETSC_COMM_WORLD, "testing MatGetDiagonal(), MatDiagonalScale()\n");
119: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
120: VecSet(x, one);
121: MatGetDiagonal(C, x);
122: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
123: for (i = vstart; i < vend; i++) {
124: v = one * ((PetscReal)(i + 1));
125: VecSetValues(y, 1, &i, &v, INSERT_VALUES);
126: }
128: /* -------------------- Test () MatDiagonalScale ------------------ */
129: PetscOptionsHasName(NULL, NULL, "-test_diagonalscale", &flg);
130: if (flg) {
131: MatDiagonalScale(C, x, y);
132: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
133: }
134: /* Free data structures */
135: VecDestroy(&u);
136: VecDestroy(&s);
137: VecDestroy(&w);
138: VecDestroy(&x);
139: VecDestroy(&y);
140: 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: 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*/