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: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
20: n = m;
21: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
22: if (flg) n += 2;
23: PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
24: if (flg) n -= 2;
26: /* ---------- Assemble matrix and vectors ----------- */
28: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
29: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m, n));
30: PetscCall(MatSetFromOptions(C));
31: PetscCall(MatSetUp(C));
32: PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
33: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
34: PetscCall(VecSetSizes(x, PETSC_DECIDE, m));
35: PetscCall(VecSetFromOptions(x));
36: PetscCall(VecDuplicate(x, &z));
37: PetscCall(VecDuplicate(x, &w));
38: PetscCall(VecCreate(PETSC_COMM_WORLD, &y));
39: PetscCall(VecSetSizes(y, PETSC_DECIDE, n));
40: PetscCall(VecSetFromOptions(y));
41: PetscCall(VecDuplicate(y, &u));
42: PetscCall(VecDuplicate(y, &s));
43: PetscCall(VecGetOwnershipRange(y, &vstart, &vend));
45: /* Assembly */
46: for (i = rstart; i < rend; i++) {
47: v = 100 * (i + 1);
48: PetscCall(VecSetValues(z, 1, &i, &v, INSERT_VALUES));
49: for (j = 0; j < n; j++) {
50: v = 10 * (i + 1) + j + 1;
51: PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES));
52: }
53: }
55: /* Flush off proc Vec values and do more assembly */
56: PetscCall(VecAssemblyBegin(z));
57: for (i = vstart; i < vend; i++) {
58: v = one * ((PetscReal)i);
59: PetscCall(VecSetValues(y, 1, &i, &v, INSERT_VALUES));
60: v = 100.0 * i;
61: PetscCall(VecSetValues(u, 1, &i, &v, INSERT_VALUES));
62: }
64: /* Flush off proc Mat values and do more assembly */
65: PetscCall(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: PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES));
70: }
71: }
72: /* Try overlap Coomunication with the next stage XXXSetValues */
73: PetscCall(VecAssemblyEnd(z));
75: PetscCall(MatAssemblyEnd(C, MAT_FLUSH_ASSEMBLY));
76: CHKMEMQ;
77: /* The Assembly for the second Stage */
78: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
79: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
80: PetscCall(VecAssemblyBegin(y));
81: PetscCall(VecAssemblyEnd(y));
82: PetscCall(MatScale(C, alpha));
83: PetscCall(VecAssemblyBegin(u));
84: PetscCall(VecAssemblyEnd(u));
85: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMult()\n"));
86: CHKMEMQ;
87: PetscCall(MatMult(C, y, x));
88: CHKMEMQ;
89: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
90: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultAdd()\n"));
91: PetscCall(MatMultAdd(C, y, z, w));
92: PetscCall(VecAXPY(x, one, z));
93: PetscCall(VecAXPY(x, negone, w));
94: PetscCall(VecNorm(x, NORM_2, &norm));
95: if (norm > tol) PetscCall(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: PetscCall(VecSetValues(x, 1, &i, &v, INSERT_VALUES));
102: }
103: PetscCall(VecAssemblyBegin(x));
104: PetscCall(VecAssemblyEnd(x));
105: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTranspose()\n"));
106: PetscCall(MatMultTranspose(C, x, y));
107: PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
109: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatMultTransposeAdd()\n"));
110: PetscCall(MatMultTransposeAdd(C, x, u, s));
111: PetscCall(VecAXPY(y, one, u));
112: PetscCall(VecAXPY(y, negone, s));
113: PetscCall(VecNorm(y, NORM_2, &norm));
114: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error difference = %g\n", (double)norm));
116: /* -------------------- Test MatGetDiagonal() ------------------ */
118: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "testing MatGetDiagonal(), MatDiagonalScale()\n"));
119: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
120: PetscCall(VecSet(x, one));
121: PetscCall(MatGetDiagonal(C, x));
122: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
123: for (i = vstart; i < vend; i++) {
124: v = one * ((PetscReal)(i + 1));
125: PetscCall(VecSetValues(y, 1, &i, &v, INSERT_VALUES));
126: }
128: /* -------------------- Test () MatDiagonalScale ------------------ */
129: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_diagonalscale", &flg));
130: if (flg) {
131: PetscCall(MatDiagonalScale(C, x, y));
132: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
133: }
134: /* Free data structures */
135: PetscCall(VecDestroy(&u));
136: PetscCall(VecDestroy(&s));
137: PetscCall(VecDestroy(&w));
138: PetscCall(VecDestroy(&x));
139: PetscCall(VecDestroy(&y));
140: PetscCall(VecDestroy(&z));
141: PetscCall(MatDestroy(&C));
143: PetscCall(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
221: test:
222: suffix: 33
223: nsize: 3
224: args: -mat_type mpiaij -test_diagonalscale
225: filter: grep -v type
227: test:
228: suffix: 34
229: nsize: 3
230: args: -mat_type mpibaij -test_diagonalscale
232: test:
233: suffix: 3_aijcusparse_1
234: args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
235: filter: grep -v type
236: output_file: output/ex5_31.out
237: requires: cuda
239: test:
240: suffix: 3_aijcusparse_2
241: nsize: 3
242: args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
243: filter: grep -v type
244: output_file: output/ex5_33.out
245: requires: cuda
247: test:
248: suffix: 3_kokkos
249: nsize: 3
250: args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale
251: filter: grep -v type
252: output_file: output/ex5_33.out
253: requires: kokkos_kernels
255: test:
256: suffix: aijcusparse_1
257: args: -mat_type seqaijcusparse -vec_type cuda -rectA
258: filter: grep -v type
259: output_file: output/ex5_11_A.out
260: requires: cuda
262: test:
263: suffix: aijcusparse_2
264: args: -mat_type seqaijcusparse -vec_type cuda -rectB
265: filter: grep -v type
266: output_file: output/ex5_11_B.out
267: requires: cuda
269: test:
270: suffix: sell_1
271: args: -mat_type sell
272: output_file: output/ex5_41.out
274: test:
275: suffix: sell_2
276: nsize: 3
277: args: -mat_type sell
278: output_file: output/ex5_43.out
280: test:
281: suffix: sell_3
282: args: -mat_type sell -test_diagonalscale
283: output_file: output/ex5_51.out
285: test:
286: suffix: sell_4
287: nsize: 3
288: args: -mat_type sell -test_diagonalscale
289: output_file: output/ex5_53.out
291: TEST*/