ViennaSHE 1.3.0
Free open-source semiconductor device simulator using spherical harmonics expansions techniques.
petsc.hpp
Go to the documentation of this file.
1/* ============================================================================
2 Copyright (c) 2011-2019, Institute for Microelectronics,
3 Institute for Analysis and Scientific Computing,
4 TU Wien.
5
6 -----------------
7 ViennaSHE - The Vienna Spherical Harmonics Expansion Boltzmann Solver
8 -----------------
9
10 http://viennashe.sourceforge.net/
11
12 Lead Developer: Karl Rupp rupp@iue.tuwien.ac.at
13 Developer: Felipe Senra Riberio ribeiro@iue.tuwien.ac.at
14
15 License: MIT (X11), see file LICENSE in the base directory
16 =============================================================================== */
17
18#ifndef PETSC_HPP_
19#define PETSC_HPP_
20
21#include <petscksp.h>
22#include <mpi.h>
23#include <iterator>
24#include <vector>
25#include <numeric>
26#include <memory>
27#include <unordered_map>
28
29namespace viennashe
30{
31 namespace solvers
32 {
36 template<typename NumericT, typename VectorType>
37 class PETSC_vector;
38 template<typename NumericT, typename VectorType>
39 class PETSC_matrix;
40
41 template<typename NumericT, typename VectorType>
43 {
44 private:
45
46 PETSC_handler (int argc, char **argv) :
47 comm (MPI_COMM_WORLD), self (PETSC_COMM_SELF)
48 {
49 MPI_Comm_size (PETSC_COMM_WORLD, &sizeW);
50 MPI_Comm_rank (PETSC_COMM_WORLD, &rankW);
51 MPI_Comm_size (PETSC_COMM_SELF, &size);
52 MPI_Comm_rank (PETSC_COMM_SELF, &rank);
53 }
54 void
55 operator= (PETSC_handler const&); // Don't implement
56 static std::unique_ptr<PETSC_handler> singleton;
57 MPI_Comm comm, self;
58 PetscViewer viewer;
59 PetscMPIInt rankW, sizeW, rank, size;
60
61 public:
62 // Implement methods for setting the PETSC env.
63 static PETSC_handler&
65 {
66 if (singleton == nullptr)
67 singleton.reset (new PETSC_handler);
68 return *singleton;
69 }
70
71 static PETSC_handler&
72 get_handler (int argc, char **argv)
73 {
74 if (singleton == nullptr)
75 singleton.reset (new PETSC_handler (argc, argv));
76 return *singleton;
77 }
78 PETSC_handler (PETSC_handler const&) = delete;
79
80 MPI_Comm
81 getComm () const
82 {
83 return comm;
84 }
85
86 void
87 setComm (MPI_Comm comm)
88 {
89 this->comm = comm;
90 }
91
92 PetscMPIInt
93 getRank () const
94 {
95 return rank;
96 }
97
98 PetscMPIInt
99 getRankW () const
100 {
101 return rankW;
102 }
103
104 MPI_Comm
105 getSelf () const
106 {
107 return self;
108 }
109
110 PetscMPIInt
111 getSize () const
112 {
113 return size;
114 }
115
116 PetscMPIInt
117 getSizeW () const
118 {
119 return sizeW;
120 }
121 }
122 ;
123
124 template<typename NumericT, typename VectorType>
125 std::unique_ptr<PETSC_handler<NumericT, VectorType>> PETSC_handler<
126 NumericT, VectorType>::singleton = nullptr;
127
130 template<typename NumericT, typename VectorType>
132 {
134 public:
135
136 typedef std::size_t size_type;
138
139 public:
140
145 comm (PETSC_COMM_WORLD), self (
146 PETSC_COMM_SELF),its (0),i (0), j (0), index (0), Creason (KSP_CONVERGED_RTOL_NORMAL), iter (0)
147 {
149 // typedef typename RowType::iterator iterator2;
150 typedef typename RowType::const_iterator const_iterator2;
151
152 //* Get some information from the matrix*/
153 MPI_Comm_size (PETSC_COMM_WORLD, &size);
154 MPI_Comm_rank (PETSC_COMM_WORLD, &rank);
155 int n = copymat.size1 ();
156 int m = copymat.size2 ();
157 long unsigned int local = copymat.size2();
158 rstart = 0, rend =m ;
159 size1_ = n;
160 size2_ = m;
161 long unsigned int nnz = 0;
162 for (unsigned long int k = 0; k < copymat.size1(); k++)
163 {
164 nnz = nnz < copymat.row(k).size()? copymat.row(k).size() :nnz;
165 }
166 nnz++;
167 PetscInt *col = new PetscInt[nnz];
168 nghosts = 0;
169 double *value = new double[nnz];
170 this->index = 0;
171 ierr = MatCreate (PETSC_COMM_WORLD, &Apetsc);
172 CHKERRV (ierr);
173
174 //* Split matrix for SHE only//
175 if(PETSC_matrix::ninstan == 2){
176 long unsigned int initial;
177 long unsigned int total;
178 RowType const &Row = copymat.row (local-1);
179 offset = 0;
180 if(rank != size-1) offset = Row.rbegin()->first -Row.begin()->first ;
181 local = copymat.size2();
182 MPI_Scan(&local, &initial, 1, MPI_UNSIGNED_LONG , MPI_SUM, PETSC_COMM_WORLD);
183 MPI_Allreduce(&local, &total, 1, MPI_UNSIGNED_LONG , MPI_SUM, PETSC_COMM_WORLD);
184 PETSC_matrix::ninstan = -1;
185 PetscBarrier (NULL);
186 m = n=size2_ = total;
187 rstart = initial-local;
188 rend = initial;
189 ierr = MatSetSizes (Apetsc,local, local, PETSC_DECIDE, PETSC_DECIDE);
190 CHKERRV (ierr);
191
192 }else{
193 ierr = MatSetSizes (Apetsc, PETSC_DECIDE, PETSC_DECIDE, m, n);
194 offset = 0;
195 CHKERRV (ierr);
196 }
197
198 ierr = MatSetFromOptions (Apetsc);
199 CHKERRV (ierr);
200 ierr = MatSeqAIJSetPreallocation (Apetsc, nnz, PETSC_NULL);
201 CHKERRV (ierr);
202 ierr = MatMPIAIJSetPreallocation (Apetsc, nnz, PETSC_NULL, nnz,
203 PETSC_NULL);
204 CHKERRV (ierr);
205
206 ierr = MatSetUp (Apetsc);
207 CHKERRV (ierr);
208 ierr = MatGetLocalSize (Apetsc, &nlocal1, &nlocal2);
209 CHKERRV (ierr);
210 ierr = MatSetOption (Apetsc, MAT_NEW_NONZERO_ALLOCATION_ERR,
211 PETSC_FALSE);
212 PetscInt nlocal;
213 for (PetscInt k = 0; k < (PetscInt)(copymat.size1()); k++)
214 {
215 RowType const &Row = copymat.row (k);
216 long nnzt = 0;
217
218 for (const_iterator2 it = Row.begin (); it != Row.end (); it++)
219 {
220 col[nnzt] = it->first+rstart;
221 value[nnzt] = it->second;
222 nlocal = k+rstart;
223 nnzt++;
224 }
225 ierr = MatSetValues (Apetsc, 1, &(nlocal), nnzt, col, value, INSERT_VALUES);
226 CHKERRV (ierr);
227 }
228
229 delete[] col;
230 delete[] value;
231 ierr = MatAssemblyBegin (Apetsc, MAT_FINAL_ASSEMBLY);
232 CHKERRV (ierr);
233 ierr = MatAssemblyEnd (Apetsc, MAT_FINAL_ASSEMBLY);
234 CHKERRV (ierr);
235 nconsiter++;
236 KSPCreate (comm, &ksp);
237 KSPSetOperators (ksp, Apetsc, Apetsc);
238 PETSC_matrix::ninstan++;
239 }
240
242 {
243 MatDestroy (&Apetsc);
244 KSPDestroy (&ksp);
245 }
246
250 void
252 {
253 PetscErrorCode ierr;
254 KSPSetFromOptions (ksp);
256
257 ierr = KSPSolve (ksp, b.getVec (), b.getSolution ());CHKERRV (ierr);
258 iter++;
259 KSPGetIterationNumber (ksp, &its);
260 b.GenerateGlobal ();
261 }
262
264 size1 () const
265 {
266 return size1_;
267 }
269 size2 () const
270 {
271 return size2_;
272 }
273
274
275
278
279 self_type&
280 operator = (double other) // note: argument passed by value!
281 {
282 ierr = MatSetValue (Apetsc, i, j, other, INSERT_VALUES);
283 return *this;
284 }
285
286 PetscReal
287 getNorm () const
288 {
289 return norm;
290 }
291 KSPConvergedReason
293 {
294 KSPGetConvergedReason (ksp, &Creason);
295 return Creason;
296 }
297
298 int
300 {
301 KSPGetIterationNumber (ksp, &its);
302 return its;
303 }
304 void
305 setNorm (PetscReal norm)
306 {
307 this->norm = norm;
308 }
309
311 getSize1 () const
312 {
313 return size1_;
314 }
315
316 void
318 {
319 size1_ = size1;
320 }
321
323 getSize2 () const
324 {
325 return size2_;
326 }
327
328 void
330 {
331 size2_ = size2;
332 }
333
334 const Mat&
336 {
337 return Apetsc;
338 }
339
340 int
341 getIndex () const
342 {
343 return index;
344 }
345
346 void
347 setIndex (int index)
348 {
349 this->index = index;
350 }
351
352 PetscInt
353 getNlocal1 () const
354 {
355 return nlocal1;
356 }
357
358 void
359 setNlocal1 (PetscInt nlocal1)
360 {
361 this->nlocal1 = nlocal1;
362 }
363
364 PetscInt
365 getNlocal2 () const
366 {
367 return nlocal2;
368 }
369
370 void
371 setNlocal2 (PetscInt nlocal2)
372 {
373 this->nlocal2 = nlocal2;
374 }
375
376 int
377 getIter () const
378 {
379 return iter;
380 }
381
382 void
383 setIter (int iter)
384 {
385 this->iter = iter;
386 }
387
388 PetscInt
389 getRend () const
390 {
391 return rend;
392 }
393
394 void
395 setRend (PetscInt rend)
396 {
397 this->rend = rend;
398 }
399
400 PetscInt
401 getRstart () const
402 {
403 return rstart;
404 }
405
406 void
407 setRstart (PetscInt rstart)
408 {
409 this->rstart = rstart;
410 }
411
412 private:
413 PetscMPIInt rank, size;
414 PetscInt rstart , rend, offset ;
415 MPI_Comm comm, self;
416 IS is;
417 MatPartitioning part;
418 PetscInt its, nlocal1, nlocal2;
419 long i, j, index;
420 Mat Apetsc;
421 KSP ksp;
422 // PC pc;
423 PetscReal norm;
424 KSPConvergedReason Creason;
425 PetscErrorCode ierr = 0;
426 size_type size1_;
427 size_type size2_;
428 int nghosts, iter;
429 static int ninstan;
430 static int nconsiter;
431
432 };
435 template<typename NumericT, typename VectorType>
437 {
439 public:
440 typedef std::size_t size_type;
441 typedef NumericT value_type;
442
446 PETSC_vector (const std::vector<NumericT> &V,
448 i(0),nlocal(0), residualb (false), globalb (false),size_n (V.size ()), dad (&M)
449 {
450 PetscInt nlocal;
451 ierr = MatCreateVecs ((M.getA ()), &vec, &solution);
452 VecDuplicate (vec, &residual);
453 ierr = PetscObjectSetName ((PetscObject) vec, "RHS");
454 ierr = PetscObjectSetName ((PetscObject) solution, "Solution");
455 VecGetOwnershipRange (vec, &rstart, &rend);
456 VecGetLocalSize (vec, &nlocal);
457 sizeg = M.size2();
458 rstart = M.getRstart();
459 rend = M.getRend();
460 for (size_t k = 0; k < V.size(); k++)
461 {
462 if (V[k] != 0.0){
463 nnz++;
464 }
465 nlocal = k+rstart;
466 VecSetValues (vec, 1, &(nlocal), &V[k], INSERT_VALUES);
467 }
468 VecAssemblyBegin (vec);
469 VecAssemblyEnd (vec);
470 VecAssemblyBegin (solution);
471 VecAssemblyEnd (solution);
472 VecAssemblyBegin (residual);
473 VecAssemblyEnd (residual);
474 }
476 {
477 VecDestroy (&vec);
478 VecDestroy (&solution);
479 VecDestroy (&residual);
480 }
481
485 operator std::vector<NumericT> () // const
486 {
487 std::vector<NumericT> result (size_n); // static or not
488 PetscScalar *array;
489 PetscInt size;
490 if (!globalb)
491 {
492 VecGetLocalSize (solution, &size);
493 if (residualb == false)
494 {
495 VecGetArray (solution, &array);
496 std::vector<NumericT> result (array, array + (size));
497 VecRestoreArray (solution, &array);
498 VecDestroy (&residual);
499 VecDestroy (&solution);
500 VecDestroy (&vec);
501 return result;
502 }
503 else
504 {
505 VecGetLocalSize (residual, &size);
506 MatMult (dad->getA (), vec, residual); // y <- Ax
507 VecAXPY (residual, -1., vec); // r <- r - b
508 VecGetArray (residual, &array);
509 std::vector<NumericT> result (array, array + (size));
510 VecRestoreArray (residual, &array);
511 return result;
512 }
513 }
514 else
515 {
516 VecGetLocalSize (global, &size);
517 if (size != 0)
518 {
519 if ((residualb == false) && (size != 0))
520 {
521 VecGetArray (global, &array);
522 std::vector<NumericT> result (array+rstart,array+rend);
523 VecRestoreArray (global, &array);
524 VecScatterDestroy (&scat);
525 VecDestroy (&global);
526 VecDestroy (&residual);
527 VecDestroy (&solution);
528 VecDestroy (&vec);
529 return result;
530 }
531 else
532 {
533 VecGetLocalSize (residual, &size);
534 MatMult (dad->getA (), vec, residual); // y <- Ax
535 VecAXPY (residual, -1., vec); // r <- r - b
536 VecGetArray (residual, &array);
537 std::vector<NumericT> result (array, array + (size));
538 VecRestoreArray (residual, &array);
539 return result;
540 }
541 }
542
543 }
544 return result;
545 }
546// /** @brief Send the solution to the first process( process 0)
547// * **/
548 void
550 {
551 globalb = true;
552 VecScatterCreateToAll (solution, &scat, &global);
553 ierr = PetscObjectSetName ((PetscObject) global, "Global");
554 VecScatterBegin (scat, solution, global, INSERT_VALUES,
555 SCATTER_FORWARD);
556 VecScatterEnd (scat, solution, global, INSERT_VALUES,
557 SCATTER_FORWARD);
558 }
559
560 Vec
561 getSolution () const
562 {
563 return solution;
564 }
565
566 void
567 setSolution (Vec solution)
568 {
569 this->solution = solution;
570 }
571 operator Vec () const
572 {
573 return vec;
574 }
575
577 size () const
578 {
579 return this->size_n;
580 }
581
583 sizeLocal () const
584 {
585 return this->nlocal;
586 }
587
588 Vec
589 getVec () const
590 {
591 return vec;
592 }
593
594 void
595 setVec (Vec vec)
596 {
597 this->vec = vec;
598 }
599
600 private:
601 int i;
602 int nnz, nlocal;
603 bool residualb, globalb;
604 PetscInt rstart , rend ;
605 IS indexG, indexL;
606 Vec vec, solution, residual, global;
607 VecScatter scatter;
608 size_type size_n, sizeg;
609 VecScatter scat;
610 ISLocalToGlobalMapping ltog;
612 PetscErrorCode ierr;
613
614 };
615
616 template<typename NumericT, typename VectorType>
617 int PETSC_matrix<NumericT, VectorType>::ninstan = 0;
618 template<typename NumericT, typename VectorType>
619 int PETSC_matrix<NumericT, VectorType>::nconsiter = 0;
620 }
621}
622#endif /* PETSC_HPP_ */
Sparse matrix class based on a vector of binary trees for holding the entries.
Definition: linalg_util.hpp:69
MapType const & row(size_t i) const
PetscMPIInt getSizeW() const
Definition: petsc.hpp:117
static PETSC_handler & get_handler()
Definition: petsc.hpp:64
PETSC_handler(PETSC_handler const &)=delete
MPI_Comm getComm() const
Definition: petsc.hpp:81
PetscMPIInt getRankW() const
Definition: petsc.hpp:99
PetscMPIInt getRank() const
Definition: petsc.hpp:93
static PETSC_handler & get_handler(int argc, char **argv)
Definition: petsc.hpp:72
void setComm(MPI_Comm comm)
Definition: petsc.hpp:87
PetscMPIInt getSize() const
Definition: petsc.hpp:111
Wrapper class to the PETSC matrix.
Definition: petsc.hpp:132
size_type getSize2() const
Definition: petsc.hpp:323
PETSC_vector< NumericT, VectorType > row_type
Definition: petsc.hpp:276
self_type & operator=(double other)
Definition: petsc.hpp:280
PetscInt getNlocal1() const
Definition: petsc.hpp:353
size_type size1() const
Definition: petsc.hpp:264
PetscInt getNlocal2() const
Definition: petsc.hpp:365
size_type size2() const
Definition: petsc.hpp:269
PetscInt getRend() const
Definition: petsc.hpp:389
void setNlocal1(PetscInt nlocal1)
Definition: petsc.hpp:359
PETSC_vector< NumericT, VectorType > RowType
Definition: petsc.hpp:277
void setRend(PetscInt rend)
Definition: petsc.hpp:395
void setNlocal2(PetscInt nlocal2)
Definition: petsc.hpp:371
KSPConvergedReason get_Creason()
Definition: petsc.hpp:292
PetscInt getRstart() const
Definition: petsc.hpp:401
void setSize2(size_type size2)
Definition: petsc.hpp:329
void setNorm(PetscReal norm)
Definition: petsc.hpp:305
void setSize1(size_type size1)
Definition: petsc.hpp:317
void setRstart(PetscInt rstart)
Definition: petsc.hpp:407
PetscReal getNorm() const
Definition: petsc.hpp:287
PETSC_matrix(const viennashe::math::sparse_matrix< NumericT > &copymat)
Wrapper class to the PETSC matrix.
Definition: petsc.hpp:144
void solve(PETSC_vector< NumericT, VectorType > &b)
Wrapper class to the PETSC matrix.
Definition: petsc.hpp:251
size_type getSize1() const
Definition: petsc.hpp:311
Class Encapsulate The PETSc COMM WORLD -Should be singleton.
Definition: petsc.hpp:437
size_type sizeLocal() const
Definition: petsc.hpp:583
size_type size() const
Definition: petsc.hpp:577
void setSolution(Vec solution)
Definition: petsc.hpp:567
PETSC_vector(const std::vector< NumericT > &V, PETSC_matrix< NumericT, VectorType > &M)
Constructor to the PETSC vector Right Hand and Solution M Matrix of the Linear System.
Definition: petsc.hpp:446
The main ViennaSHE namespace. All functionality resides inside this namespace.
Definition: accessors.hpp:40