1/*
2 * SparseMatrix_petsc.hpp
3 *
4 * Created on: Apr 26, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
9#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_
10
11#include "util/petsc_util.hpp"
12#include "Vector/map_vector.hpp"
13#include <boost/mpl/int.hpp>
14#include <petscmat.h>
15
16#define PETSC_BASE 2
17
18/*! \brief It store one non-zero element in the sparse matrix
19 *
20 * Given a row, and a column, store a value
21 *
22 *
23 */
24template<typename T>
25class triplet<T,PETSC_BASE>
26{
27 //! Row of the sparse matrix
28 PetscInt row_;
29
30 //! Colum of the sparse matrix
31 PetscInt col_;
32
33 //! Value of the Matrix
34 PetscScalar val_;
35
36public:
37
38 /*! \brief Return the row of the triplet
39 *
40 * \return the row index
41 *
42 */
43 PetscInt & row()
44 {
45 return row_;
46 }
47
48 /*! \brief Return the colum of the triplet
49 *
50 * \return the colum index
51 *
52 */
53 PetscInt & col()
54 {
55 return col_;
56 }
57
58 /*! \brief Return the value of the triplet
59 *
60 * \return the value
61 *
62 */
63 PetscScalar & value()
64 {
65 return val_;
66 }
67
68 /*! \brief Constructor from row, colum and value
69 *
70 * \param i row
71 * \param j colum
72 * \param val value
73 *
74 */
75 triplet(long int i, long int j, T val)
76 {
77 row_ = i;
78 col_ = j;
79 val_ = val;
80 }
81
82 // Default constructor
83 triplet()
84 :row_(0),col_(0),val_(0)
85 {};
86};
87
88/*! \brief Sparse Matrix implementation, that map over Eigen
89 *
90 * \tparam T Type of the sparse Matrix store on each row,colums
91 * \tparam id_t type of id
92 * \tparam impl implementation
93 *
94 */
95template<typename T, typename id_t>
96class SparseMatrix<T,id_t, PETSC_BASE>
97{
98public:
99
100 //! Triplet implementation id
101 typedef boost::mpl::int_<PETSC_BASE> triplet_impl;
102
103 //! Triplet type
104 typedef triplet<T,PETSC_BASE> triplet_type;
105
106private:
107
108 //! Number of matrix row (global)
109 size_t g_row;
110 //! Number of matrix colums (global)
111 size_t g_col;
112
113 //! Number of matrix row (local)
114 size_t l_row;
115 //! Number of matrix colums (local)
116 size_t l_col;
117
118 //! starting row for this processor
119 size_t start_row;
120
121 //! indicate if the matrix has been created
122 bool m_created = false;
123
124 //! PETSC Matrix
125 Mat mat;
126
127 //! Triplets of the matrix
128 openfpm::vector<triplet_type> trpl;
129
130
131 //! temporary list of values
132 mutable openfpm::vector<PetscScalar> vals;
133 //! temporary list of colums
134 mutable openfpm::vector<PetscInt> cols;
135 //! PETSC d_nnz
136 mutable openfpm::vector<PetscInt> d_nnz;
137 //! PETSC o_nnz
138 mutable openfpm::vector<PetscInt> o_nnz;
139
140 /*! \brief Fill the petsc Matrix
141 *
142 *
143 */
144 void fill_petsc()
145 {
146 d_nnz.resize(l_row);
147 o_nnz.resize(l_row);
148
149 d_nnz.fill(0);
150 o_nnz.fill(0);
151
152 // Here we explore every row to count how many non zero we have in the diagonal matrix part,
153 // and the non diagonal matrix part, needed by MatMPIAIJSetPreallocation
154
155 size_t i = 0;
156
157 // Set the Matrix from triplet
158 while (i < trpl.size())
159 {
160 PetscInt row = trpl.get(i).row();
161
162 while(i < trpl.size() && row == trpl.get(i).row())
163 {
164 if ((size_t)trpl.get(i).col() >= start_row && (size_t)trpl.get(i).col() < start_row + l_row)
165 d_nnz.get(row - start_row)++;
166 else
167 o_nnz.get(row - start_row)++;
168 i++;
169 }
170 }
171
172 PETSC_SAFE_CALL(MatMPIAIJSetPreallocation(mat,0,static_cast<const PetscInt*>(d_nnz.getPointer()),0,
173 static_cast<const PetscInt*>(o_nnz.getPointer())));
174
175 // Counter i is zero
176 i = 0;
177
178 // Set the Matrix from triplet
179 while (i < trpl.size())
180 {
181 vals.clear();
182 cols.clear();
183
184 PetscInt row = trpl.get(i).row();
185
186 while(i < trpl.size() && row == trpl.get(i).row())
187 {
188 vals.add(trpl.get(i).value());
189 cols.add(trpl.get(i).col());
190 i++;
191 }
192 PETSC_SAFE_CALL(MatSetValues(mat,1,&row,cols.size(),static_cast<const PetscInt*>(cols.getPointer()),
193 static_cast<const PetscScalar *>(vals.getPointer()),
194 INSERT_VALUES));
195 }
196
197 PETSC_SAFE_CALL(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY));
198 PETSC_SAFE_CALL(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY));
199
200 m_created = true;
201 }
202
203
204public:
205
206 /*! \brief Create an empty Matrix
207 *
208 * \param N1 number of row
209 * \param N2 number of colums
210 * \param N1_loc number of local row
211 *
212 */
213 SparseMatrix(size_t N1, size_t N2, size_t n_row_local)
214 :g_row(N1),g_col(N2),l_row(n_row_local),l_col(n_row_local)
215 {
216 PETSC_SAFE_CALL(MatCreate(PETSC_COMM_WORLD,&mat));
217 PETSC_SAFE_CALL(MatSetType(mat,MATMPIAIJ));
218 PETSC_SAFE_CALL(MatSetSizes(mat,n_row_local,n_row_local,N1,N2));
219
220 Vcluster<> & v_cl = create_vcluster();
221
222 openfpm::vector<size_t> vn_row_local;
223 v_cl.allGather(l_row,vn_row_local);
224 v_cl.execute();
225
226 // Calculate the starting row for this processor
227
228 start_row = 0;
229 for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
230 start_row += vn_row_local.get(i);
231 }
232
233 /*! \brief Create an empty Matrix
234 *
235 */
236 SparseMatrix()
237 :g_row(0),g_col(0),l_row(0l),l_col(0),start_row(0)
238 {
239 PETSC_SAFE_CALL(MatCreate(PETSC_COMM_WORLD,&mat));
240 PETSC_SAFE_CALL(MatSetType(mat,MATMPIAIJ));
241 }
242
243 ~SparseMatrix()
244 {
245 // Destroy the matrix
246 if (is_openfpm_init() == true)
247 {PETSC_SAFE_CALL(MatDestroy(&mat));}
248 }
249
250 /*! \brief Get the Matrix triplets buffer
251 *
252 * It return a buffer that can be filled with triplets
253 *
254 * \return Petsc Matrix
255 *
256 */
257 openfpm::vector<triplet_type> & getMatrixTriplets()
258 {
259 m_created = false;
260
261 return this->trpl;
262 }
263
264 /*! \brief Get the Patsc Matrix object
265 *
266 * \return the Eigen Matrix
267 *
268 */
269 const Mat & getMat() const
270 {
271 if (m_created == false)
272 {fill_petsc();}
273
274 return mat;
275 }
276
277 /*! \brief Get the Petsc Matrix object
278 *
279 * \return the Petsc Matrix
280 *
281 */
282 Mat & getMat()
283 {
284 if (m_created == false)
285 {fill_petsc();}
286
287 return mat;
288 }
289
290 /*! \brief Resize the Sparse Matrix
291 *
292 * \param row number for row
293 * \param col number of colums
294 * \param local number of row
295 * \param local number of colums
296 *
297 */
298 void resize(size_t row, size_t col, size_t l_row, size_t l_col)
299 {
300 this->g_row = row;
301 this->g_col = col;
302
303 this->l_row = l_row;
304 this->l_col = l_col;
305
306 PETSC_SAFE_CALL(MatSetSizes(mat,l_row,l_col,g_row,g_col));
307
308 Vcluster<> & v_cl = create_vcluster();
309
310 openfpm::vector<size_t> vn_row_local;
311 v_cl.allGather(l_row,vn_row_local);
312 v_cl.execute();
313
314 // Calculate the starting row for this processor
315
316 start_row = 0;
317 for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
318 start_row += vn_row_local.get(i);
319 }
320
321 /*! \brief Get the row i and the colum j of the Matrix
322 *
323 * \warning it is slow, consider to get blocks of the matrix
324 *
325 * \return the value of the matrix at row i colum j
326 *
327 */
328 T operator()(id_t i, id_t j)
329 {
330 T v;
331
332 MatGetValues(mat,1,&i,1,&j,&v);
333
334 return v;
335 }
336
337 /*! \brief Get the value from triplet
338 *
339 * \warning It is extremly slow because it do a full search across the triplets elements
340 *
341 * \param r row
342 * \param c colum
343 *
344 */
345 T getValue(size_t r, size_t c)
346 {
347 for (size_t i = 0 ; i < trpl.size() ; i++)
348 {
349 if (r == (size_t)trpl.get(i).row() && c == (size_t)trpl.get(i).col())
350 return trpl.get(i).value();
351 }
352
353 return 0;
354 }
355};
356
357
358#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_PETSC_HPP_ */
359