1/*
2 * Vector_petsc.hpp
3 *
4 * Created on: Apr 29, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_PETSC_HPP_
9#define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_PETSC_HPP_
10
11#include "Vector/map_vector.hpp"
12#include "Vector/vector_def.hpp"
13#include <boost/mpl/int.hpp>
14#include <petscvec.h>
15#include "util/petsc_util.hpp"
16#include <unordered_map>
17
18#define PETSC_RVAL 2
19
20/*! \brief It store one row value of a vector
21 *
22 * Given a row, store a value
23 *
24 *
25 */
26template<typename T>
27class rval<T,PETSC_RVAL>
28{
29public:
30
31 //! boost fusion that store the point
32 typedef boost::fusion::vector<PetscInt,T> type;
33
34 //! structure that store the data of the point
35 type data;
36
37 //! Property id for the row
38 static const unsigned int row = 0;
39
40 //! Property id for the value
41 static const unsigned int value = 1;
42
43 //! This object has 2 properties
44 static const unsigned int max_prop = 2;
45
46 /*! \brief Get the row
47 *
48 * \return the row
49 *
50 */
51 long int & rw()
52 {
53 return boost::fusion::at_c<row>(data);
54 }
55
56 /*! \brief Get the value
57 *
58 * \return the value
59 *
60 */
61 T & val()
62 {
63 return boost::fusion::at_c<value>(data);
64 }
65
66 /*! \brief Default constructor
67 *
68 */
69 rval() {}
70
71 /*! \brief Constructor from row, column and value
72 *
73 * \param i row
74 * \param val value
75 *
76 */
77 rval(long int i, T val)
78 {
79 rw() = i;
80 val() = val;
81 }
82
83 /*! \brief Indicate that the structure has no pointer
84 *
85 * \return true
86 *
87 */
88 static inline bool noPointers()
89 {
90 return true;
91 }
92};
93
94constexpr unsigned int row_id = 0;
95constexpr unsigned int val_id = 1;
96
97
98/*! \brief PETSC vector for linear algebra
99 *
100 * This vector wrap the PETSC vector for solving linear systems
101 *
102 */
103template<typename T>
104class Vector<T,PETSC_BASE>
105{
106 //! Number of row the petsc vector has
107 size_t n_row;
108
109 //! Number of local rows
110 size_t n_row_local;
111
112 //! Indicate if v has been allocated
113 mutable bool v_created = false;
114
115 //! Mutable vector
116 mutable Vec v;
117
118 //! Mutable row value vector
119 mutable openfpm::vector<rval<PetscScalar,PETSC_RVAL>,HeapMemory, memory_traits_inte > row_val;
120
121 //! Global to local map
122 mutable std::unordered_map<size_t,size_t> map;
123
124 //! invalid
125 T invalid;
126
127 /*! \brief Set the Eigen internal vector
128 *
129 *
130 */
131 void setPetsc() const
132 {
133 if (v_created == false)
134 {PETSC_SAFE_CALL(VecSetType(v,VECMPI));}
135
136 // set the vector
137
138 if (row_val.size() != 0)
139 {PETSC_SAFE_CALL(VecSetValues(v,row_val.size(),&row_val.template get<row_id>(0),&row_val.template get<val_id>(0),INSERT_VALUES))}
140
141 PETSC_SAFE_CALL(VecAssemblyBegin(v));
142 PETSC_SAFE_CALL(VecAssemblyEnd(v));
143
144 v_created = true;
145 }
146
147public:
148
149 /*! \brief Copy the vector
150 *
151 * \param v vector to copy
152 *
153 */
154 Vector(const Vector<T,PETSC_BASE> & v)
155 {
156 this->operator=(v);
157 }
158
159 /*! \brief Copy the vector
160 *
161 * \param v vector to copy
162 *
163 */
164 Vector(Vector<T,PETSC_BASE> && v)
165 :n_row(0),n_row_local(0),invalid(0)
166 {
167 this->operator=(v);
168 }
169
170 /*! \brief Destroy the vector
171 *
172 *
173 */
174 ~Vector()
175 {
176 if (is_openfpm_init() == true)
177 {PETSC_SAFE_CALL(VecDestroy(&v));}
178 }
179
180 /*! \brief Create a vector with n elements
181 *
182 * \param n global number of elements in the vector
183 * \param n_row_local number
184 *
185 */
186 Vector(size_t n, size_t n_row_local)
187 :n_row_local(n_row_local),v(NULL),invalid(0)
188 {
189 // Create the vector
190 PETSC_SAFE_CALL(VecCreate(PETSC_COMM_WORLD,&v));
191
192 resize(n,n_row_local);
193 }
194
195 /*! \brief Create a vector with 0 elements
196 *
197 */
198 Vector()
199 :n_row(0),n_row_local(0),invalid(0)
200 {
201 // Create the vector
202 PETSC_SAFE_CALL(VecCreate(PETSC_COMM_WORLD,&v));
203 }
204
205 /*! \brief Resize the Vector
206 *
207 * \param row numbers of row
208 * \param l_row number of local row
209 *
210 */
211 void resize(size_t row, size_t l_row)
212 {
213 n_row = row;
214 n_row_local = l_row;
215
216 PETSC_SAFE_CALL(VecSetSizes(v,n_row_local,n_row));
217 }
218
219 /*! \brief Return a reference to the vector element
220 *
221 * \param i element
222 * \param val value
223 *
224 */
225 void insert(size_t i, T val)
226 {
227 row_val.add();
228
229 // Map
230 map[i] = row_val.size()-1;
231
232 row_val.last().template get<row_id>() = i;
233 row_val.last().template get<val_id>() = val;
234 }
235
236 /*! \brief Return a reference to the vector element
237 *
238 * \param i element
239 *
240 * \return reference to the element vector
241 *
242 */
243 inline PetscScalar & insert(size_t i)
244 {
245 row_val.add();
246
247 // Map
248 map[i] = row_val.size()-1;
249
250 row_val.last().template get<row_id>() = i;
251 return row_val.last().template get<val_id>();
252 }
253
254 /*! \brief Return a reference to the vector element
255 *
256 * \param i element
257 *
258 * \return reference to the element vector
259 *
260 */
261 inline const PetscScalar & insert(size_t i) const
262 {
263 row_val.add();
264
265 // Map
266 map[i] = row_val.size()-1;
267
268 row_val.last().template get<row_id>() = i;
269 return row_val.last().template get<val_id>();
270 }
271
272 /*! \brief Return a reference to the vector element
273 *
274 * \warning The element must exist
275 *
276 * \param i element
277 *
278 * \return reference to the element vector
279 *
280 */
281 const PetscScalar & operator()(size_t i) const
282 {
283 // Search if exist
284
285 std::unordered_map<size_t,size_t>::iterator it = map.find(i);
286
287 if ( it != map.end() )
288 return row_val.template get<val_id>(it->second);
289
290 return insert(i);
291 }
292
293 /*! \brief Return a reference to the vector element
294 *
295 * \warning The element must exist
296 *
297 * \param i element
298 *
299 * \return reference to the element vector
300 *
301 */
302 PetscScalar & operator()(size_t i)
303 {
304 // Search if exist
305
306 std::unordered_map<size_t,size_t>::iterator it = map.find(i);
307
308 if ( it != map.end() )
309 return row_val.template get<val_id>(it->second);
310
311 return insert(i);
312 }
313
314 /*! \brief Get the PETSC Vector object
315 *
316 * \return the PETSC Vector
317 *
318 */
319 const Vec & getVec() const
320 {
321 setPetsc();
322
323 return v;
324 }
325
326 /*! \brief Get the PETSC Vector object
327 *
328 * \return the PETSC Vector
329 *
330 */
331 Vec & getVec()
332 {
333 setPetsc();
334
335 return v;
336 }
337
338 /*! \brief Update the Vector with the PETSC object
339 *
340 */
341 void update()
342 {
343 PetscInt n_row;
344 PetscInt n_row_local;
345
346 // Get the size of the vector from PETSC
347 VecGetSize(v,&n_row);
348 VecGetLocalSize(v,&n_row_local);
349
350 this->n_row = n_row;
351 this->n_row_local = n_row_local;
352
353 row_val.resize(n_row_local);
354
355 //
356
357 PetscInt low;
358 PetscInt high;
359
360 VecGetOwnershipRange(v,&low,&high);
361
362 // Fill the index and construct the map
363
364 size_t k = 0;
365 for (size_t i = low ; i < (size_t)high ; i++)
366 {
367 row_val.template get<row_id>(k) = i;
368 map[i] = k;
369 k++;
370 }
371
372 PETSC_SAFE_CALL(VecGetValues(v,row_val.size(),&row_val.template get<row_id>(0),&row_val.template get<val_id>(0)))
373 }
374
375 /*! \brief Copy the vector
376 *
377 * \param v vector to copy
378 *
379 */
380 Vector<T,PETSC_BASE> & operator=(const Vector<T,PETSC_BASE> & v)
381 {
382 map = v.map;
383 row_val = v.row_val;
384
385 return *this;
386 }
387
388 /*! \brief Copy the vector
389 *
390 * \param v vector to copy
391 *
392 */
393 Vector<T,PETSC_BASE> & operator=(Vector<T,PETSC_BASE> && v)
394 {
395 map.swap(v.map);
396 row_val.swap(v.row_val);
397
398 return *this;
399 }
400
401 /*! \brief Set to zero all the entries
402 *
403 *
404 */
405 void setZero()
406 {
407 if (v_created == false)
408 {PETSC_SAFE_CALL(VecSetType(v,VECMPI));}
409
410 v_created = true;
411 }
412};
413
414
415#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_EIGEN_HPP_ */
416
417