1 static char help[] = "Test ViennaCL Matrix Conversions";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 PetscMPIInt rank, size;
8
9 PetscFunctionBeginUser;
10 PetscCall(PetscInitialize(&argc, &args, NULL, help));
11
12 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
13 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
14
15 /* Construct a sequential ViennaCL matrix and vector */
16 if (rank == 0) {
17 Mat A_vcl;
18 Vec v_vcl, r_vcl;
19 PetscInt n = 17, m = 31, nz = 5, i, cnt, j;
20 const PetscScalar val = 1.0;
21
22 PetscCall(MatCreateSeqAIJViennaCL(PETSC_COMM_SELF, m, n, nz, NULL, &A_vcl));
23
24 /* Add nz arbitrary entries per row in arbitrary columns */
25 for (i = 0; i < m; ++i) {
26 for (cnt = 0; cnt < nz; ++cnt) {
27 j = (19 * cnt + (7 * i + 3)) % n;
28 PetscCall(MatSetValue(A_vcl, i, j, (PetscScalar)(0.3 * i + j), INSERT_VALUES));
29 }
30 }
31 PetscCall(MatAssemblyBegin(A_vcl, MAT_FINAL_ASSEMBLY));
32 PetscCall(MatAssemblyEnd(A_vcl, MAT_FINAL_ASSEMBLY));
33
34 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, n, &v_vcl));
35 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, m, &r_vcl));
36 PetscCall(VecSet(v_vcl, val));
37
38 PetscCall(MatMult(A_vcl, v_vcl, r_vcl));
39
40 PetscCall(VecDestroy(&v_vcl));
41 PetscCall(VecDestroy(&r_vcl));
42 PetscCall(MatDestroy(&A_vcl));
43 }
44
45 /* Create a sequential AIJ matrix on rank 0 convert it to a new ViennaCL matrix, and apply it to a ViennaCL vector */
46 if (rank == 0) {
47 Mat A, A_vcl;
48 Vec v, r, v_vcl, r_vcl, d_vcl;
49 PetscInt n = 17, m = 31, nz = 5, i, cnt, j;
50 const PetscScalar val = 1.0;
51 PetscReal dnorm;
52 const PetscReal tol = 1e-5;
53
54 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, nz, NULL, &A));
55
56 /* Add nz arbitrary entries per row in arbitrary columns */
57 for (i = 0; i < m; ++i) {
58 for (cnt = 0; cnt < nz; ++cnt) {
59 j = (19 * cnt + (7 * i + 3)) % n;
60 PetscCall(MatSetValue(A, i, j, (PetscScalar)(0.3 * i + j), INSERT_VALUES));
61 }
62 }
63 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
64 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
65
66 PetscCall(PetscObjectSetName((PetscObject)A, "Sequential CPU Matrix"));
67
68 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v));
69 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &r));
70 PetscCall(PetscObjectSetName((PetscObject)r, "CPU result vector"));
71 PetscCall(VecSet(v, val));
72 PetscCall(MatMult(A, v, r));
73
74 PetscCall(MatConvert(A, MATSEQAIJVIENNACL, MAT_INITIAL_MATRIX, &A_vcl));
75 PetscCall(PetscObjectSetName((PetscObject)A_vcl, "New ViennaCL Matrix"));
76
77 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, n, &v_vcl));
78 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, m, &r_vcl));
79 PetscCall(PetscObjectSetName((PetscObject)r_vcl, "ViennaCL result vector"));
80 PetscCall(VecSet(v_vcl, val));
81 PetscCall(MatMult(A_vcl, v_vcl, r_vcl));
82
83 PetscCall(VecDuplicate(r_vcl, &d_vcl));
84 PetscCall(VecCopy(r_vcl, d_vcl));
85 PetscCall(VecAXPY(d_vcl, -1.0, r_vcl));
86 PetscCall(VecNorm(d_vcl, NORM_INFINITY, &dnorm));
87 PetscCheck(dnorm <= tol, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Sequential CPU and MPI ViennaCL vector results incompatible with inf norm greater than tolerance of %g", tol);
88
89 PetscCall(VecDestroy(&v));
90 PetscCall(VecDestroy(&r));
91 PetscCall(VecDestroy(&v_vcl));
92 PetscCall(VecDestroy(&r_vcl));
93 PetscCall(VecDestroy(&d_vcl));
94 PetscCall(MatDestroy(&A));
95 PetscCall(MatDestroy(&A_vcl));
96 }
97
98 /* Create a sequential AIJ matrix on rank 0 convert it inplace to a new ViennaCL matrix, and apply it to a ViennaCL vector */
99 if (rank == 0) {
100 Mat A;
101 Vec v, r, v_vcl, r_vcl, d_vcl;
102 PetscInt n = 17, m = 31, nz = 5, i, cnt, j;
103 const PetscScalar val = 1.0;
104 PetscReal dnorm;
105 const PetscReal tol = 1e-5;
106
107 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, nz, NULL, &A));
108
109 /* Add nz arbitrary entries per row in arbitrary columns */
110 for (i = 0; i < m; ++i) {
111 for (cnt = 0; cnt < nz; ++cnt) {
112 j = (19 * cnt + (7 * i + 3)) % n;
113 PetscCall(MatSetValue(A, i, j, (PetscScalar)(0.3 * i + j), INSERT_VALUES));
114 }
115 }
116 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
117 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
118
119 PetscCall(PetscObjectSetName((PetscObject)A, "Sequential CPU Matrix"));
120
121 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v));
122 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &r));
123 PetscCall(PetscObjectSetName((PetscObject)r, "CPU result vector"));
124 PetscCall(VecSet(v, val));
125 PetscCall(MatMult(A, v, r));
126
127 PetscCall(MatConvert(A, MATSEQAIJVIENNACL, MAT_INPLACE_MATRIX, &A));
128 PetscCall(PetscObjectSetName((PetscObject)A, "Converted ViennaCL Matrix"));
129
130 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, n, &v_vcl));
131 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, m, &r_vcl));
132 PetscCall(PetscObjectSetName((PetscObject)r_vcl, "ViennaCL result vector"));
133 PetscCall(VecSet(v_vcl, val));
134 PetscCall(MatMult(A, v_vcl, r_vcl));
135
136 PetscCall(VecDuplicate(r_vcl, &d_vcl));
137 PetscCall(VecCopy(r_vcl, d_vcl));
138 PetscCall(VecAXPY(d_vcl, -1.0, r_vcl));
139 PetscCall(VecNorm(d_vcl, NORM_INFINITY, &dnorm));
140 PetscCheck(dnorm <= tol, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g", tol);
141
142 PetscCall(VecDestroy(&v));
143 PetscCall(VecDestroy(&r));
144 PetscCall(VecDestroy(&v_vcl));
145 PetscCall(VecDestroy(&r_vcl));
146 PetscCall(VecDestroy(&d_vcl));
147 PetscCall(MatDestroy(&A));
148 }
149
150 /* Create a parallel AIJ matrix, convert it to a new ViennaCL matrix, and apply it to a ViennaCL vector */
151 if (size > 1) {
152 Mat A, A_vcl;
153 Vec v, r, v_vcl, r_vcl, d_vcl;
154 PetscInt N = 17, M = 31, nz = 5, i, cnt, j, rlow, rhigh;
155 const PetscScalar val = 1.0;
156 PetscReal dnorm;
157 const PetscReal tol = 1e-5;
158
159 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, M, N, nz, NULL, nz, NULL, &A));
160
161 /* Add nz arbitrary entries per row in arbitrary columns */
162 PetscCall(MatGetOwnershipRange(A, &rlow, &rhigh));
163 for (i = rlow; i < rhigh; ++i) {
164 for (cnt = 0; cnt < nz; ++cnt) {
165 j = (19 * cnt + (7 * i + 3)) % N;
166 PetscCall(MatSetValue(A, i, j, (PetscScalar)(0.3 * i + j), INSERT_VALUES));
167 }
168 }
169 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
170 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
171
172 PetscCall(PetscObjectSetName((PetscObject)A, "MPI CPU Matrix"));
173
174 PetscCall(MatCreateVecs(A, &v, &r));
175 PetscCall(PetscObjectSetName((PetscObject)r, "MPI CPU result vector"));
176 PetscCall(VecSet(v, val));
177 PetscCall(MatMult(A, v, r));
178
179 PetscCall(MatConvert(A, MATMPIAIJVIENNACL, MAT_INITIAL_MATRIX, &A_vcl));
180 PetscCall(PetscObjectSetName((PetscObject)A_vcl, "New MPI ViennaCL Matrix"));
181
182 PetscCall(MatCreateVecs(A_vcl, &v_vcl, &r_vcl));
183 PetscCall(PetscObjectSetName((PetscObject)r_vcl, "ViennaCL result vector"));
184 PetscCall(VecSet(v_vcl, val));
185 PetscCall(MatMult(A_vcl, v_vcl, r_vcl));
186
187 PetscCall(VecDuplicate(r_vcl, &d_vcl));
188 PetscCall(VecCopy(r_vcl, d_vcl));
189 PetscCall(VecAXPY(d_vcl, -1.0, r_vcl));
190 PetscCall(VecNorm(d_vcl, NORM_INFINITY, &dnorm));
191 PetscCheck(dnorm <= tol, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g", tol);
192
193 PetscCall(VecDestroy(&v));
194 PetscCall(VecDestroy(&r));
195 PetscCall(VecDestroy(&v_vcl));
196 PetscCall(VecDestroy(&r_vcl));
197 PetscCall(VecDestroy(&d_vcl));
198 PetscCall(MatDestroy(&A));
199 PetscCall(MatDestroy(&A_vcl));
200 }
201
202 /* Create a parallel AIJ matrix, convert it in-place to a ViennaCL matrix, and apply it to a ViennaCL vector */
203 if (size > 1) {
204 Mat A;
205 Vec v, r, v_vcl, r_vcl, d_vcl;
206 PetscInt N = 17, M = 31, nz = 5, i, cnt, j, rlow, rhigh;
207 const PetscScalar val = 1.0;
208 PetscReal dnorm;
209 const PetscReal tol = 1e-5;
210
211 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, M, N, nz, NULL, nz, NULL, &A));
212
213 /* Add nz arbitrary entries per row in arbitrary columns */
214 PetscCall(MatGetOwnershipRange(A, &rlow, &rhigh));
215 for (i = rlow; i < rhigh; ++i) {
216 for (cnt = 0; cnt < nz; ++cnt) {
217 j = (19 * cnt + (7 * i + 3)) % N;
218 PetscCall(MatSetValue(A, i, j, (PetscScalar)(0.3 * i + j), INSERT_VALUES));
219 }
220 }
221 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
222 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
223
224 PetscCall(PetscObjectSetName((PetscObject)A, "MPI CPU Matrix"));
225
226 PetscCall(MatCreateVecs(A, &v, &r));
227 PetscCall(PetscObjectSetName((PetscObject)r, "MPI CPU result vector"));
228 PetscCall(VecSet(v, val));
229 PetscCall(MatMult(A, v, r));
230
231 PetscCall(MatConvert(A, MATMPIAIJVIENNACL, MAT_INPLACE_MATRIX, &A));
232 PetscCall(PetscObjectSetName((PetscObject)A, "Converted MPI ViennaCL Matrix"));
233
234 PetscCall(MatCreateVecs(A, &v_vcl, &r_vcl));
235 PetscCall(PetscObjectSetName((PetscObject)r_vcl, "ViennaCL result vector"));
236 PetscCall(VecSet(v_vcl, val));
237 PetscCall(MatMult(A, v_vcl, r_vcl));
238
239 PetscCall(VecDuplicate(r_vcl, &d_vcl));
240 PetscCall(VecCopy(r_vcl, d_vcl));
241 PetscCall(VecAXPY(d_vcl, -1.0, r_vcl));
242 PetscCall(VecNorm(d_vcl, NORM_INFINITY, &dnorm));
243 PetscCheck(dnorm <= tol, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MPI CPU and MPI ViennaCL Vector results incompatible with inf norm greater than tolerance of %g", tol);
244
245 PetscCall(VecDestroy(&v));
246 PetscCall(VecDestroy(&r));
247 PetscCall(VecDestroy(&v_vcl));
248 PetscCall(VecDestroy(&r_vcl));
249 PetscCall(VecDestroy(&d_vcl));
250 PetscCall(MatDestroy(&A));
251 }
252
253 PetscCall(PetscFinalize());
254 return 0;
255 }
256
257 /*TEST
258
259 build:
260 requires: viennacl
261
262 test:
263 output_file: output/empty.out
264
265 test:
266 suffix: 2
267 nsize: 2
268 output_file: output/empty.out
269
270 TEST*/
271