OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
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 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
181BOOST_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
352BOOST_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
427BOOST_AUTO_TEST_SUITE_END()
428
429
430#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_UNIT_TESTS_HPP_ */
SparseMatrix
Sparse Matrix implementation.
Definition SparseMatrix.hpp:79
Vcluster_base::getProcessUnitID
size_t getProcessUnitID()
Get the process unit id.
Definition VCluster_base.hpp:535
Vcluster_base::getProcessingUnits
size_t getProcessingUnits()
Get the total number of processors.
Definition VCluster_base.hpp:479
Vcluster
Implementation of VCluster class.
Definition VCluster.hpp:59
Vector
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition Vector.hpp:40
petsc_solver
In case T does not match the PETSC precision compilation create a stub structure.
Definition petsc_solver.hpp:56
petsc_solver::solve
static Vector< T > solve(const SparseMatrix< T, impl > &A, const Vector< T > &b)
Solve the linear system.
Definition petsc_solver.hpp:69
umfpack_solver
stub when library compiled without eigen
Definition umfpack_solver.hpp:194
umfpack_solver::solve
static Vector< T > solve(const SparseMatrix< T, id_type, impl > &A, const Vector< T, impl > &b)
stub solve
Definition umfpack_solver.hpp:198
triplet
It store the non zero elements of the matrix.
Definition SparseMatrix.hpp:38