OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
SparseMatrix_unit_tests.cpp
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 
23 BOOST_AUTO_TEST_SUITE( sparse_matrix_test_suite )
24 
25 
26 
27 BOOST_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  umfpack_solver<double> solver;
153  auto x = solver.solve(sm,v);
154 
155  // we control the solution
156 
157  if (vcl.getProcessUnitID() == 0)
158  {
159  BOOST_REQUIRE_CLOSE(x(0), -4.5, 0.001);
160  BOOST_REQUIRE_CLOSE(x(1), -8, 0.001);
161  BOOST_REQUIRE_CLOSE(x(2), -10.5, 0.001);
162  }
163  else if (vcl.getProcessUnitID() == 1)
164  {
165  BOOST_REQUIRE_CLOSE(x(3), -12.0, 0.001);
166  BOOST_REQUIRE_CLOSE(x(4), -12.5, 0.001);
167  BOOST_REQUIRE_CLOSE(x(5), -12.0, 0.001);
168  }
169  else if (vcl.getProcessUnitID() == 2)
170  {
171  BOOST_REQUIRE_CLOSE(x(6), -10.5, 0.001);
172  BOOST_REQUIRE_CLOSE(x(7), -8, 0.001);
173  BOOST_REQUIRE_CLOSE(x(8), -4.5, 0.001);
174  }
175 
176 #endif
177 }
178 
179 #ifdef HAVE_PETSC
180 
181 BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc)
182 {
183  Vcluster<> & vcl = create_vcluster();
184 
185  if (vcl.getProcessingUnits() != 3)
186  return;
187 
188  // 3 Processors 9x9 Matrix to invert
189 
190  SparseMatrix<double,int,PETSC_BASE> sm(9,9,3);
191  Vector<double,PETSC_BASE> v(9,3);
192 
193  typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
194 
195  auto & triplets = sm.getMatrixTriplets();
196 
197  if (vcl.getProcessUnitID() == 0)
198  {
199  // row 1
200  triplets.add(triplet(0,0,-2));
201  triplets.add(triplet(0,1,1));
202 
203  // row 2
204  triplets.add(triplet(1,0,1));
205  triplets.add(triplet(1,1,-2));
206  triplets.add(triplet(1,2,1));
207 
208  // row 3
209  triplets.add(triplet(2,1,1));
210  triplets.add(triplet(2,2,-2));
211  triplets.add(triplet(2,3,1));
212 
213  // v row1
214  v.insert(0,1);
215  v.insert(1,1);
216  v.insert(2,1);
217  }
218  else if (vcl.getProcessUnitID() == 1)
219  {
220  // row 4
221  triplets.add(triplet(3,2,1));
222  triplets.add(triplet(3,3,-2));
223  triplets.add(triplet(3,4,1));
224 
225  // row 5
226  triplets.add(triplet(4,3,1));
227  triplets.add(triplet(4,4,-2));
228  triplets.add(triplet(4,5,1));
229 
230  // row 6
231  triplets.add(triplet(5,4,1));
232  triplets.add(triplet(5,5,-2));
233  triplets.add(triplet(5,6,1));
234 
235  v.insert(3,1);
236  v.insert(4,1);
237  v.insert(5,1);
238  }
239  else if (vcl.getProcessUnitID() == 2)
240  {
241  // row 7
242  triplets.add(triplet(6,5,1));
243  triplets.add(triplet(6,6,-2));
244  triplets.add(triplet(6,7,1));
245 
246  // row 8
247  triplets.add(triplet(7,6,1));
248  triplets.add(triplet(7,7,-2));
249  triplets.add(triplet(7,8,1));
250 
251  // row 9
252  triplets.add(triplet(8,7,1));
253  triplets.add(triplet(8,8,-2));
254 
255  v.insert(6,1);
256  v.insert(7,1);
257  v.insert(8,1);
258  }
259 
260  // Get the petsc Matrix
261  Mat & matp = sm.getMat();
262 
263  if (vcl.getProcessUnitID() == 0)
264  {
265  PetscInt r0[1] = {0};
266  PetscInt r1[1] = {1};
267  PetscInt r2[1] = {2};
268  PetscInt r0cols[3] = {0,1};
269  PetscInt r1cols[3] = {0,1,2};
270  PetscInt r2cols[3] = {1,2,3};
271  PetscScalar y[3];
272 
273  MatGetValues(matp,1,r0,2,r0cols,y);
274 
275  BOOST_REQUIRE_EQUAL(y[0],-2);
276  BOOST_REQUIRE_EQUAL(y[1],1);
277 
278  MatGetValues(matp,1,r1,3,r1cols,y);
279 
280  BOOST_REQUIRE_EQUAL(y[0],1);
281  BOOST_REQUIRE_EQUAL(y[1],-2);
282  BOOST_REQUIRE_EQUAL(y[2],1);
283 
284  MatGetValues(matp,1,r2,3,r2cols,y);
285 
286  BOOST_REQUIRE_EQUAL(y[0],1);
287  BOOST_REQUIRE_EQUAL(y[1],-2);
288  BOOST_REQUIRE_EQUAL(y[2],1);
289 
290  }
291  else if (vcl.getProcessUnitID() == 1)
292  {
293  PetscInt r0[1] = {3};
294  PetscInt r1[1] = {4};
295  PetscInt r2[1] = {5};
296  PetscInt r0cols[3] = {2,3,4};
297  PetscInt r1cols[3] = {3,4,5};
298  PetscInt r2cols[3] = {4,5,6};
299  PetscScalar y[3];
300 
301  MatGetValues(matp,1,r0,3,r0cols,y);
302 
303  BOOST_REQUIRE_EQUAL(y[0],1);
304  BOOST_REQUIRE_EQUAL(y[1],-2);
305  BOOST_REQUIRE_EQUAL(y[2],1);
306 
307  MatGetValues(matp,1,r1,3,r1cols,y);
308 
309  BOOST_REQUIRE_EQUAL(y[0],1);
310  BOOST_REQUIRE_EQUAL(y[1],-2);
311  BOOST_REQUIRE_EQUAL(y[2],1);
312 
313  MatGetValues(matp,1,r2,3,r2cols,y);
314 
315  BOOST_REQUIRE_EQUAL(y[0],1);
316  BOOST_REQUIRE_EQUAL(y[1],-2);
317  BOOST_REQUIRE_EQUAL(y[2],1);
318  }
319  else if (vcl.getProcessUnitID() == 2)
320  {
321  PetscInt r0[1] = {6};
322  PetscInt r1[1] = {7};
323  PetscInt r2[1] = {8};
324  PetscInt r0cols[3] = {5,6,7};
325  PetscInt r1cols[3] = {6,7,8};
326  PetscInt r2cols[3] = {7,8};
327  PetscScalar y[3];
328 
329  MatGetValues(matp,1,r0,3,r0cols,y);
330 
331  BOOST_REQUIRE_EQUAL(y[0],1);
332  BOOST_REQUIRE_EQUAL(y[1],-2);
333  BOOST_REQUIRE_EQUAL(y[2],1);
334 
335  MatGetValues(matp,1,r1,3,r1cols,y);
336 
337  BOOST_REQUIRE_EQUAL(y[0],1);
338  BOOST_REQUIRE_EQUAL(y[1],-2);
339  BOOST_REQUIRE_EQUAL(y[2],1);
340 
341  MatGetValues(matp,1,r2,2,r2cols,y);
342 
343  BOOST_REQUIRE_EQUAL(y[0],1);
344  BOOST_REQUIRE_EQUAL(y[1],-2);
345 
346  }
347 
348  petsc_solver<double> solver;
349  solver.solve(sm,v);
350 }
351 
352 BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc_solve)
353 {
354  Vcluster<> & vcl = create_vcluster();
355 
356  if (vcl.getProcessingUnits() != 3)
357  return;
358 
359  const int loc = 200;
360 
361  // 3 Processors 9x9 Matrix to invert
362 
363  SparseMatrix<double,int,PETSC_BASE> sm(3*loc,3*loc,loc);
364  Vector<double,PETSC_BASE> v(3*loc,loc);
365 
366  typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
367 
368  auto & triplets = sm.getMatrixTriplets();
369 
370  if (vcl.getProcessUnitID() == 0)
371  {
372  // row 1
373  triplets.add(triplet(0,0,-2));
374  triplets.add(triplet(0,1,1));
375  v.insert(0,1);
376 
377  // row 2
378  for (size_t i = 1 ; i < loc ; i++)
379  {
380  triplets.add(triplet(i,i-1,1));
381  triplets.add(triplet(i,i,-2));
382  triplets.add(triplet(i,i+1,1));
383 
384  v.insert(i,1);
385  }
386  }
387  else if (vcl.getProcessUnitID() == 1)
388  {
389  // row 2
390  for (size_t i = loc ; i < 2*loc ; i++)
391  {
392  triplets.add(triplet(i,i-1,1));
393  triplets.add(triplet(i,i,-2));
394  triplets.add(triplet(i,i+1,1));
395 
396  v.insert(i,1);
397  }
398  }
399  else if (vcl.getProcessUnitID() == 2)
400  {
401  // row 2
402  for (size_t i = 2*loc ; i < 3*loc-1 ; i++)
403  {
404  triplets.add(triplet(i,i-1,1));
405  triplets.add(triplet(i,i,-2));
406  triplets.add(triplet(i,i+1,1));
407 
408  v.insert(i,1);
409  }
410 
411  // row 9
412  triplets.add(triplet(3*loc-1,3*loc-2,1));
413  triplets.add(triplet(3*loc-1,3*loc-1,-2));
414  v.insert(3*loc-1,1);
415  }
416 
417  // Get the petsc Matrix
418  sm.getMat();
419 
420 
421  petsc_solver<double> solver;
422  solver.solve(sm,v);
423 }
424 
425 #endif
426 
427 BOOST_AUTO_TEST_SUITE_END()
428 
429 
430 #endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_UNIT_TESTS_HPP_ */
Vcluster_base::getProcessUnitID
size_t getProcessUnitID()
Get the process unit id.
Definition: VCluster_base.hpp:532
Vector
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition: Vector.hpp:39
triplet
It store the non zero elements of the matrix.
Definition: SparseMatrix.hpp:37
petsc_solver< double >::solve
Vector< double, PETSC_BASE > solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, bool initial_guess=false)
Here we invert the matrix and solve the system.
Definition: petsc_solver.hpp:1306
umfpack_solver< double >
stub when library compiled without eigen
Definition: umfpack_solver.hpp:218
Vcluster
Implementation of VCluster class.
Definition: VCluster.hpp:58
SparseMatrix
Sparse Matrix implementation.
Definition: SparseMatrix.hpp:78
petsc_solver< double >
This class is able to do Matrix inversion in parallel with PETSC solvers.
Definition: petsc_solver.hpp:101
Vcluster_base::getProcessingUnits
size_t getProcessingUnits()
Get the total number of processors.
Definition: VCluster_base.hpp:476
umfpack_solver< double >::solve
static Vector< double > solve(SparseMatrix< double, id_type, impl > &A, const Vector< double > &b, size_t opt=UMFPACK_NONE)
stub solve
Definition: umfpack_solver.hpp:224