1/*
2 * SparseMatrix_unit_tests.hpp
3 *
4 * Created on: Apr 4, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_UNIT_TESTS_HPP_
9#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_UNIT_TESTS_HPP_
10
11#define BOOST_TEST_DYN_LINK
12#include <boost/test/unit_test.hpp>
13
14#include "Matrix/SparseMatrix.hpp"
15#include "Vector/Vector.hpp"
16#include "Solvers/umfpack_solver.hpp"
17#include "Solvers/petsc_solver.hpp"
18
19#ifdef HAVE_PETSC
20#include <petscmat.h>
21#endif
22
23BOOST_AUTO_TEST_SUITE( sparse_matrix_test_suite )
24
25
26
27BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_parallel)
28{
29#if defined(HAVE_EIGEN) && defined(HAVE_SUITESPARSE)
30
31 Vcluster<> & vcl = create_vcluster();
32
33 if (vcl.getProcessingUnits() != 3)
34 return;
35
36 // 3 Processors 9x9 Matrix to invert
37
38 SparseMatrix<double,int> sm(9,9);
39 Vector<double> v(9);
40
41 typedef SparseMatrix<double,int>::triplet_type triplet;
42
43 auto & triplets = sm.getMatrixTriplets();
44
45 if (vcl.getProcessUnitID() == 0)
46 {
47 // row 1
48 triplets.add(triplet(0,0,-2));
49 triplets.add(triplet(0,1,1));
50
51 // row 2
52 triplets.add(triplet(1,0,1));
53 triplets.add(triplet(1,1,-2));
54 triplets.add(triplet(1,2,1));
55
56 // row 3
57 triplets.add(triplet(2,1,1));
58 triplets.add(triplet(2,2,-2));
59 triplets.add(triplet(2,3,1));
60
61 // v row1
62 v.insert(0,1);
63 v.insert(1,1);
64 v.insert(2,1);
65 }
66 else if (vcl.getProcessUnitID() == 1)
67 {
68 // row 4
69 triplets.add(triplet(3,2,1));
70 triplets.add(triplet(3,3,-2));
71 triplets.add(triplet(3,4,1));
72
73 // row 5
74 triplets.add(triplet(4,3,1));
75 triplets.add(triplet(4,4,-2));
76 triplets.add(triplet(4,5,1));
77
78 // row 6
79 triplets.add(triplet(5,4,1));
80 triplets.add(triplet(5,5,-2));
81 triplets.add(triplet(5,6,1));
82
83 v.insert(3,1);
84 v.insert(4,1);
85 v.insert(5,1);
86 }
87 else if (vcl.getProcessUnitID() == 2)
88 {
89 // row 7
90 triplets.add(triplet(6,5,1));
91 triplets.add(triplet(6,6,-2));
92 triplets.add(triplet(6,7,1));
93
94 // row 8
95 triplets.add(triplet(7,6,1));
96 triplets.add(triplet(7,7,-2));
97 triplets.add(triplet(7,8,1));
98
99 // row 9
100 triplets.add(triplet(8,7,1));
101 triplets.add(triplet(8,8,-2));
102
103 v.insert(6,1);
104 v.insert(7,1);
105 v.insert(8,1);
106 }
107
108 // force to sync
109 sm.getMat();
110
111 // Master has the full Matrix
112
113 if (vcl.getProcessUnitID() == 0)
114 {
115 BOOST_REQUIRE_EQUAL(sm(0,0),-2);
116 BOOST_REQUIRE_EQUAL(sm(0,1),1);
117
118 BOOST_REQUIRE_EQUAL(sm(1,0),1);
119 BOOST_REQUIRE_EQUAL(sm(1,1),-2);
120 BOOST_REQUIRE_EQUAL(sm(1,2),1);
121
122 BOOST_REQUIRE_EQUAL(sm(2,1),1);
123 BOOST_REQUIRE_EQUAL(sm(2,2),-2);
124 BOOST_REQUIRE_EQUAL(sm(2,3),1);
125
126 BOOST_REQUIRE_EQUAL(sm(3,2),1);
127 BOOST_REQUIRE_EQUAL(sm(3,3),-2);
128 BOOST_REQUIRE_EQUAL(sm(3,4),1);
129
130 BOOST_REQUIRE_EQUAL(sm(4,3),1);
131 BOOST_REQUIRE_EQUAL(sm(4,4),-2);
132 BOOST_REQUIRE_EQUAL(sm(4,5),1);
133
134 BOOST_REQUIRE_EQUAL(sm(5,4),1);
135 BOOST_REQUIRE_EQUAL(sm(5,5),-2);
136 BOOST_REQUIRE_EQUAL(sm(5,6),1);
137
138 BOOST_REQUIRE_EQUAL(sm(6,5),1);
139 BOOST_REQUIRE_EQUAL(sm(6,6),-2);
140 BOOST_REQUIRE_EQUAL(sm(6,7),1);
141
142 BOOST_REQUIRE_EQUAL(sm(7,6),1);
143 BOOST_REQUIRE_EQUAL(sm(7,7),-2);
144 BOOST_REQUIRE_EQUAL(sm(7,8),1);
145
146 BOOST_REQUIRE_EQUAL(sm(8,7),1);
147 BOOST_REQUIRE_EQUAL(sm(8,8),-2);
148 }
149
150 // try to invert the Matrix with umfpack
151
152 auto x = umfpack_solver<double>::solve(sm,v);
153
154 // we control the solution
155
156 if (vcl.getProcessUnitID() == 0)
157 {
158 BOOST_REQUIRE_CLOSE(x(0), -4.5, 0.001);
159 BOOST_REQUIRE_CLOSE(x(1), -8, 0.001);
160 BOOST_REQUIRE_CLOSE(x(2), -10.5, 0.001);
161 }
162 else if (vcl.getProcessUnitID() == 1)
163 {
164 BOOST_REQUIRE_CLOSE(x(3), -12.0, 0.001);
165 BOOST_REQUIRE_CLOSE(x(4), -12.5, 0.001);
166 BOOST_REQUIRE_CLOSE(x(5), -12.0, 0.001);
167 }
168 else if (vcl.getProcessUnitID() == 2)
169 {
170 BOOST_REQUIRE_CLOSE(x(6), -10.5, 0.001);
171 BOOST_REQUIRE_CLOSE(x(7), -8, 0.001);
172 BOOST_REQUIRE_CLOSE(x(8), -4.5, 0.001);
173 }
174
175#endif
176}
177
178#ifdef HAVE_PETSC
179
180BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc)
181{
182 Vcluster<> & vcl = create_vcluster();
183
184 if (vcl.getProcessingUnits() != 3)
185 return;
186
187 // 3 Processors 9x9 Matrix to invert
188
189 SparseMatrix<double,int,PETSC_BASE> sm(9,9,3);
190 Vector<double,PETSC_BASE> v(9,3);
191
192 typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
193
194 auto & triplets = sm.getMatrixTriplets();
195
196 if (vcl.getProcessUnitID() == 0)
197 {
198 // row 1
199 triplets.add(triplet(0,0,-2));
200 triplets.add(triplet(0,1,1));
201
202 // row 2
203 triplets.add(triplet(1,0,1));
204 triplets.add(triplet(1,1,-2));
205 triplets.add(triplet(1,2,1));
206
207 // row 3
208 triplets.add(triplet(2,1,1));
209 triplets.add(triplet(2,2,-2));
210 triplets.add(triplet(2,3,1));
211
212 // v row1
213 v.insert(0,1);
214 v.insert(1,1);
215 v.insert(2,1);
216 }
217 else if (vcl.getProcessUnitID() == 1)
218 {
219 // row 4
220 triplets.add(triplet(3,2,1));
221 triplets.add(triplet(3,3,-2));
222 triplets.add(triplet(3,4,1));
223
224 // row 5
225 triplets.add(triplet(4,3,1));
226 triplets.add(triplet(4,4,-2));
227 triplets.add(triplet(4,5,1));
228
229 // row 6
230 triplets.add(triplet(5,4,1));
231 triplets.add(triplet(5,5,-2));
232 triplets.add(triplet(5,6,1));
233
234 v.insert(3,1);
235 v.insert(4,1);
236 v.insert(5,1);
237 }
238 else if (vcl.getProcessUnitID() == 2)
239 {
240 // row 7
241 triplets.add(triplet(6,5,1));
242 triplets.add(triplet(6,6,-2));
243 triplets.add(triplet(6,7,1));
244
245 // row 8
246 triplets.add(triplet(7,6,1));
247 triplets.add(triplet(7,7,-2));
248 triplets.add(triplet(7,8,1));
249
250 // row 9
251 triplets.add(triplet(8,7,1));
252 triplets.add(triplet(8,8,-2));
253
254 v.insert(6,1);
255 v.insert(7,1);
256 v.insert(8,1);
257 }
258
259 // Get the petsc Matrix
260 Mat & matp = sm.getMat();
261
262 if (vcl.getProcessUnitID() == 0)
263 {
264 PetscInt r0[1] = {0};
265 PetscInt r1[1] = {1};
266 PetscInt r2[1] = {2};
267 PetscInt r0cols[3] = {0,1};
268 PetscInt r1cols[3] = {0,1,2};
269 PetscInt r2cols[3] = {1,2,3};
270 PetscScalar y[3];
271
272 MatGetValues(matp,1,r0,2,r0cols,y);
273
274 BOOST_REQUIRE_EQUAL(y[0],-2);
275 BOOST_REQUIRE_EQUAL(y[1],1);
276
277 MatGetValues(matp,1,r1,3,r1cols,y);
278
279 BOOST_REQUIRE_EQUAL(y[0],1);
280 BOOST_REQUIRE_EQUAL(y[1],-2);
281 BOOST_REQUIRE_EQUAL(y[2],1);
282
283 MatGetValues(matp,1,r2,3,r2cols,y);
284
285 BOOST_REQUIRE_EQUAL(y[0],1);
286 BOOST_REQUIRE_EQUAL(y[1],-2);
287 BOOST_REQUIRE_EQUAL(y[2],1);
288
289 }
290 else if (vcl.getProcessUnitID() == 1)
291 {
292 PetscInt r0[1] = {3};
293 PetscInt r1[1] = {4};
294 PetscInt r2[1] = {5};
295 PetscInt r0cols[3] = {2,3,4};
296 PetscInt r1cols[3] = {3,4,5};
297 PetscInt r2cols[3] = {4,5,6};
298 PetscScalar y[3];
299
300 MatGetValues(matp,1,r0,3,r0cols,y);
301
302 BOOST_REQUIRE_EQUAL(y[0],1);
303 BOOST_REQUIRE_EQUAL(y[1],-2);
304 BOOST_REQUIRE_EQUAL(y[2],1);
305
306 MatGetValues(matp,1,r1,3,r1cols,y);
307
308 BOOST_REQUIRE_EQUAL(y[0],1);
309 BOOST_REQUIRE_EQUAL(y[1],-2);
310 BOOST_REQUIRE_EQUAL(y[2],1);
311
312 MatGetValues(matp,1,r2,3,r2cols,y);
313
314 BOOST_REQUIRE_EQUAL(y[0],1);
315 BOOST_REQUIRE_EQUAL(y[1],-2);
316 BOOST_REQUIRE_EQUAL(y[2],1);
317 }
318 else if (vcl.getProcessUnitID() == 2)
319 {
320 PetscInt r0[1] = {6};
321 PetscInt r1[1] = {7};
322 PetscInt r2[1] = {8};
323 PetscInt r0cols[3] = {5,6,7};
324 PetscInt r1cols[3] = {6,7,8};
325 PetscInt r2cols[3] = {7,8};
326 PetscScalar y[3];
327
328 MatGetValues(matp,1,r0,3,r0cols,y);
329
330 BOOST_REQUIRE_EQUAL(y[0],1);
331 BOOST_REQUIRE_EQUAL(y[1],-2);
332 BOOST_REQUIRE_EQUAL(y[2],1);
333
334 MatGetValues(matp,1,r1,3,r1cols,y);
335
336 BOOST_REQUIRE_EQUAL(y[0],1);
337 BOOST_REQUIRE_EQUAL(y[1],-2);
338 BOOST_REQUIRE_EQUAL(y[2],1);
339
340 MatGetValues(matp,1,r2,2,r2cols,y);
341
342 BOOST_REQUIRE_EQUAL(y[0],1);
343 BOOST_REQUIRE_EQUAL(y[1],-2);
344
345 }
346
347 petsc_solver<double> solver;
348 solver.solve(sm,v);
349}
350
351BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc_solve)
352{
353 Vcluster<> & vcl = create_vcluster();
354
355 if (vcl.getProcessingUnits() != 3)
356 return;
357
358 const int loc = 200;
359
360 // 3 Processors 9x9 Matrix to invert
361
362 SparseMatrix<double,int,PETSC_BASE> sm(3*loc,3*loc,loc);
363 Vector<double,PETSC_BASE> v(3*loc,loc);
364
365 typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
366
367 auto & triplets = sm.getMatrixTriplets();
368
369 if (vcl.getProcessUnitID() == 0)
370 {
371 // row 1
372 triplets.add(triplet(0,0,-2));
373 triplets.add(triplet(0,1,1));
374 v.insert(0,1);
375
376 // row 2
377 for (size_t i = 1 ; i < loc ; i++)
378 {
379 triplets.add(triplet(i,i-1,1));
380 triplets.add(triplet(i,i,-2));
381 triplets.add(triplet(i,i+1,1));
382
383 v.insert(i,1);
384 }
385 }
386 else if (vcl.getProcessUnitID() == 1)
387 {
388 // row 2
389 for (size_t i = loc ; i < 2*loc ; i++)
390 {
391 triplets.add(triplet(i,i-1,1));
392 triplets.add(triplet(i,i,-2));
393 triplets.add(triplet(i,i+1,1));
394
395 v.insert(i,1);
396 }
397 }
398 else if (vcl.getProcessUnitID() == 2)
399 {
400 // row 2
401 for (size_t i = 2*loc ; i < 3*loc-1 ; i++)
402 {
403 triplets.add(triplet(i,i-1,1));
404 triplets.add(triplet(i,i,-2));
405 triplets.add(triplet(i,i+1,1));
406
407 v.insert(i,1);
408 }
409
410 // row 9
411 triplets.add(triplet(3*loc-1,3*loc-2,1));
412 triplets.add(triplet(3*loc-1,3*loc-1,-2));
413 v.insert(3*loc-1,1);
414 }
415
416 // Get the petsc Matrix
417 sm.getMat();
418
419
420 petsc_solver<double> solver;
421 solver.solve(sm,v);
422}
423
424#endif
425
426BOOST_AUTO_TEST_SUITE_END()
427
428
429#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_UNIT_TESTS_HPP_ */
430