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