27#include <unordered_map>
36 template<
typename NumericT,
typename VectorType>
38 template<
typename NumericT,
typename VectorType>
41 template<
typename NumericT,
typename VectorType>
47 comm (MPI_COMM_WORLD), self (PETSC_COMM_SELF)
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);
56 static std::unique_ptr<PETSC_handler> singleton;
59 PetscMPIInt rankW, sizeW, rank, size;
66 if (singleton ==
nullptr)
74 if (singleton ==
nullptr)
124 template<
typename NumericT,
typename VectorType>
125 std::unique_ptr<PETSC_handler<NumericT, VectorType>> PETSC_handler<
126 NumericT, VectorType>::singleton =
nullptr;
130 template<
typename NumericT,
typename VectorType>
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)
150 typedef typename RowType::const_iterator const_iterator2;
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 ;
161 long unsigned int nnz = 0;
162 for (
unsigned long int k = 0; k < copymat.
size1(); k++)
164 nnz = nnz < copymat.
row(k).size()? copymat.
row(k).size() :nnz;
167 PetscInt *col =
new PetscInt[nnz];
169 double *value =
new double[nnz];
171 ierr = MatCreate (PETSC_COMM_WORLD, &Apetsc);
175 if(PETSC_matrix::ninstan == 2){
176 long unsigned int initial;
177 long unsigned int total;
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;
186 m = n=size2_ = total;
187 rstart = initial-local;
189 ierr = MatSetSizes (Apetsc,local, local, PETSC_DECIDE, PETSC_DECIDE);
193 ierr = MatSetSizes (Apetsc, PETSC_DECIDE, PETSC_DECIDE, m, n);
198 ierr = MatSetFromOptions (Apetsc);
200 ierr = MatSeqAIJSetPreallocation (Apetsc, nnz, PETSC_NULL);
202 ierr = MatMPIAIJSetPreallocation (Apetsc, nnz, PETSC_NULL, nnz,
206 ierr = MatSetUp (Apetsc);
208 ierr = MatGetLocalSize (Apetsc, &nlocal1, &nlocal2);
210 ierr = MatSetOption (Apetsc, MAT_NEW_NONZERO_ALLOCATION_ERR,
213 for (PetscInt k = 0; k < (PetscInt)(copymat.
size1()); k++)
218 for (const_iterator2 it = Row.begin (); it != Row.end (); it++)
220 col[nnzt] = it->first+rstart;
221 value[nnzt] = it->second;
225 ierr = MatSetValues (Apetsc, 1, &(nlocal), nnzt, col, value, INSERT_VALUES);
231 ierr = MatAssemblyBegin (Apetsc, MAT_FINAL_ASSEMBLY);
233 ierr = MatAssemblyEnd (Apetsc, MAT_FINAL_ASSEMBLY);
236 KSPCreate (comm, &ksp);
237 KSPSetOperators (ksp, Apetsc, Apetsc);
238 PETSC_matrix::ninstan++;
243 MatDestroy (&Apetsc);
254 KSPSetFromOptions (ksp);
259 KSPGetIterationNumber (ksp, &its);
282 ierr = MatSetValue (Apetsc, i, j, other, INSERT_VALUES);
294 KSPGetConvergedReason (ksp, &Creason);
301 KSPGetIterationNumber (ksp, &its);
361 this->nlocal1 = nlocal1;
373 this->nlocal2 = nlocal2;
409 this->rstart = rstart;
413 PetscMPIInt rank, size;
414 PetscInt rstart , rend, offset ;
417 MatPartitioning part;
418 PetscInt its, nlocal1, nlocal2;
424 KSPConvergedReason Creason;
425 PetscErrorCode ierr = 0;
430 static int nconsiter;
435 template<
typename NumericT,
typename VectorType>
448 i(0),nlocal(0), residualb (false), globalb (false),size_n (V.
size ()), dad (&M)
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);
460 for (
size_t k = 0; k < V.size(); k++)
466 VecSetValues (vec, 1, &(nlocal), &V[k], INSERT_VALUES);
468 VecAssemblyBegin (vec);
469 VecAssemblyEnd (vec);
470 VecAssemblyBegin (solution);
471 VecAssemblyEnd (solution);
472 VecAssemblyBegin (residual);
473 VecAssemblyEnd (residual);
478 VecDestroy (&solution);
479 VecDestroy (&residual);
485 operator std::vector<NumericT> ()
487 std::vector<NumericT> result (size_n);
492 VecGetLocalSize (solution, &
size);
493 if (residualb ==
false)
495 VecGetArray (solution, &array);
496 std::vector<NumericT> result (array, array + (
size));
497 VecRestoreArray (solution, &array);
498 VecDestroy (&residual);
499 VecDestroy (&solution);
505 VecGetLocalSize (residual, &
size);
506 MatMult (dad->getA (), vec, residual);
507 VecAXPY (residual, -1., vec);
508 VecGetArray (residual, &array);
509 std::vector<NumericT> result (array, array + (
size));
510 VecRestoreArray (residual, &array);
516 VecGetLocalSize (global, &
size);
519 if ((residualb ==
false) && (
size != 0))
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);
533 VecGetLocalSize (residual, &
size);
534 MatMult (dad->getA (), vec, residual);
535 VecAXPY (residual, -1., vec);
536 VecGetArray (residual, &array);
537 std::vector<NumericT> result (array, array + (
size));
538 VecRestoreArray (residual, &array);
552 VecScatterCreateToAll (solution, &scat, &global);
553 ierr = PetscObjectSetName ((PetscObject) global,
"Global");
554 VecScatterBegin (scat, solution, global, INSERT_VALUES,
556 VecScatterEnd (scat, solution, global, INSERT_VALUES,
569 this->solution = solution;
571 operator Vec ()
const
603 bool residualb, globalb;
604 PetscInt rstart , rend ;
606 Vec vec, solution, residual, global;
610 ISLocalToGlobalMapping ltog;
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;
Sparse matrix class based on a vector of binary trees for holding the entries.
MapType const & row(size_t i) const
PetscMPIInt getSizeW() const
static PETSC_handler & get_handler()
PETSC_handler(PETSC_handler const &)=delete
PetscMPIInt getRankW() const
PetscMPIInt getRank() const
static PETSC_handler & get_handler(int argc, char **argv)
void setComm(MPI_Comm comm)
PetscMPIInt getSize() const
Wrapper class to the PETSC matrix.
size_type getSize2() const
PETSC_vector< NumericT, VectorType > row_type
self_type & operator=(double other)
PetscInt getNlocal1() const
PetscInt getNlocal2() const
void setNlocal1(PetscInt nlocal1)
PETSC_vector< NumericT, VectorType > RowType
void setRend(PetscInt rend)
void setNlocal2(PetscInt nlocal2)
KSPConvergedReason get_Creason()
PetscInt getRstart() const
void setSize2(size_type size2)
void setNorm(PetscReal norm)
void setSize1(size_type size1)
void setRstart(PetscInt rstart)
PetscReal getNorm() const
PETSC_matrix(const viennashe::math::sparse_matrix< NumericT > ©mat)
Wrapper class to the PETSC matrix.
void solve(PETSC_vector< NumericT, VectorType > &b)
Wrapper class to the PETSC matrix.
size_type getSize1() const
Class Encapsulate The PETSc COMM WORLD -Should be singleton.
size_type sizeLocal() const
void setSolution(Vec solution)
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.
The main ViennaSHE namespace. All functionality resides inside this namespace.