PetscVectorInline.cxx
1 // Copyright (C) 2011-2012, INRIA
2 // Author(s): Marc Fragu
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_VECTOR_PETSCVECTOR_INLINE_CXX
22 
23 
24 #include "PetscVector.hxx"
25 
26 
27 namespace Seldon
28 {
29 
30 
31  /****************
32  * CONSTRUCTORS *
33  ****************/
34 
35 
37 
40  template <class T, class Allocator>
41  inline PETScVector<T, Allocator>::PETScVector():
42  Vector_Base<T, Allocator>()
43  {
44  this->m_ = 0;
45  petsc_vector_deallocated_ = false;
46  }
47 
48 
50 
54  template <class T, class Allocator>
55  inline PETScVector<T, Allocator>::PETScVector(int i, MPI_Comm mpi_communicator)
56  :Vector_Base<T, Allocator>(i)
57  {
58  this->m_ = i;
59  petsc_vector_deallocated_ = false;
60  }
61 
62 
64 
67  template <class T, class Allocator>
68  inline PETScVector<T, Allocator>::
69  PETScVector(Vec& petsc_vector)
70  {
71  petsc_vector_deallocated_ = true;
72  Copy(petsc_vector);
73  }
74 
75 
77 
80  template <class T, class Allocator>
81  inline PETScVector<T, Allocator>::
82  PETScVector(const PETScVector<T, Allocator>& V)
83  {
84  Copy(V);
85  }
86 
87 
88  /**************
89  * DESTRUCTOR *
90  **************/
91 
92 
94  template <class T, class Allocator>
95  inline PETScVector<T, Allocator>::~PETScVector()
96  {
97  Clear();
98  }
99 
100 
102 
105  template <class T, class Allocator>
106  inline Vec& PETScVector<T, Allocator>::GetPetscVector()
107  {
108  return petsc_vector_;
109  }
110 
111 
113 
116  template <class T, class Allocator>
117  inline const Vec& PETScVector<T, Allocator>::GetPetscVector() const
118  {
119  return petsc_vector_;
120  }
121 
122 
124 
127  template <class T, class Allocator>
128  inline void PETScVector<T, Allocator>::SetCommunicator(MPI_Comm mpi_communicator)
129  {
130  mpi_communicator_ = mpi_communicator;
131  }
132 
133 
134  /*********************
135  * MEMORY MANAGEMENT *
136  *********************/
137 
138 
140 
144  template <class T, class Allocator>
145  inline void PETScVector<T, Allocator>::Clear()
146  {
147  if (petsc_vector_deallocated_)
148  return;
149  int ierr;
150  ierr = VecDestroy(&petsc_vector_);
151  CHKERRABORT(mpi_communicator_, ierr);
152  petsc_vector_deallocated_ = true;
153  }
154 
155 
157 
161  template <class T, class Allocator>
162  inline void PETScVector<T, Allocator>::Resize(int n)
163  {
164  throw Undefined("PETScVector<T, Allocator>::Resize(int n)");
165  }
166 
167 
182  template <class T, class Allocator>
183  inline void PETScVector<T, Allocator>
184  ::SetData(int i, typename PETScVector<T, Allocator>::pointer data)
185  {
186  throw Undefined("PETScVector<T, Allocator>::SetData(int i, "
187  "typename PETScVector<T, Allocator>::pointer data)");
188  }
189 
190 
192 
197  template <class T, class Allocator>
198  inline void PETScVector<T, Allocator>::Nullify()
199  {
200  throw Undefined("PETScVector<T, Allocator>::Nullify()");
201  }
202 
203 
204  /**********************************
205  * ELEMENT ACCESS AND AFFECTATION *
206  **********************************/
207 
208 
210 
214  template <class T, class Allocator>
215  inline typename PETScVector<T, Allocator>::value_type
216  PETScVector<T, Allocator>::operator() (int i) const
217  {
218 #ifdef SELDON_CHECK_BOUNDS
219  if (i < 0 || i >= this->m_)
220  throw WrongIndex("Vector<PETSc>::operator()",
221  string("Index should be in [0, ") + to_str(this->m_-1)
222  + "], but is equal to " + to_str(i) + ".");
223 #endif
224  int ierr;
225  value_type ret[1];
226  int index[1];
227  index[0] = i;
228  ierr = VecGetValues(petsc_vector_, 1, index, ret);
229  CHKERRABORT(mpi_communicator_, ierr);
230  return ret[0];
231  }
232 
233 
235 
243  template <class T, class Allocator>
244  inline void PETScVector<T, Allocator>
245  ::SetBuffer(int i, T value, InsertMode insert_mode)
246  {
247  int ierr;
248  int ix[1] = {i};
249  T data[1] = {value};
250  ierr = VecSetValues(petsc_vector_, 1, ix,
251  data, INSERT_VALUES);
252  CHKERRABORT(mpi_communicator_, ierr);
253  }
254 
255 
257  template <class T, class Allocator>
258  inline void PETScVector<T, Allocator>
259  ::Flush()
260  {
261  int ierr;
262  ierr = VecAssemblyBegin(petsc_vector_);
263  CHKERRABORT(mpi_communicator_, ierr);
264  ierr = VecAssemblyEnd(petsc_vector_);
265  CHKERRABORT(mpi_communicator_, ierr);
266  }
267 
268 
270 
278  template <class T, class Allocator>
279  inline void PETScVector<T, Allocator>
280  ::GetProcessorRange(int& i, int& j) const
281  {
282  int ierr;
283  ierr = VecGetOwnershipRange(petsc_vector_, &i, &j);
284  CHKERRABORT(mpi_communicator_, ierr);
285  }
286 
287 
289 
294  template <class T, class Allocator>
295  inline void PETScVector<T, Allocator>
296  ::Copy(const PETScVector<T, Allocator>& X)
297  {
298  Copy(X.GetPetscVector());
299  }
300 
301 
303 
308  template <class T, class Allocator>
309  inline void PETScVector<T, Allocator>::Append(const T& x)
310  {
311  throw Undefined("PETScVector<T, Allocator>::Append(const T& x)");
312  }
313 
314 
315  /*******************
316  * BASIC FUNCTIONS *
317  *******************/
318 
319 
321 
324  template <class T, class Allocator>
325  inline int PETScVector<T, Allocator>::GetDataSize() const
326  {
327  return this->m_;
328  }
329 
330 
332 
335  template <class T, class Allocator>
336  inline int PETScVector<T, Allocator>::GetLocalM() const
337  {
338  int size;
339  VecGetLocalSize(petsc_vector_, &size);
340  return size;
341  }
342 
343 
345  // SEQUENTIAL PETSC VECTOR //
347 
348 
349  /****************
350  * CONSTRUCTORS *
351  ****************/
352 
353 
355 
358  template <class T, class Allocator>
359  inline Vector<T, PETScSeq, Allocator>::Vector():
360  PETScVector<T, Allocator>()
361  {
362  this->mpi_communicator_ = MPI_COMM_WORLD;
363  this->m_ = 0;
364  int ierr;
365  ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &this->petsc_vector_);
366  this->petsc_vector_deallocated_ = false;
367  }
368 
369 
371 
375  template <class T, class Allocator>
376  inline Vector<T, PETScSeq, Allocator>::Vector(int i, MPI_Comm mpi_communicator)
377  :PETScVector<T, Allocator>(i)
378  {
379  int ierr;
380  this->mpi_communicator_ = mpi_communicator;
381  this->m_ = i;
382  ierr = VecCreateSeq(PETSC_COMM_SELF, i, &this->petsc_vector_);
383  CHKERRABORT(this->mpi_communicator_, ierr);
384  ierr = VecSet(this->petsc_vector_, 0.);
385  CHKERRABORT(this->mpi_communicator_, ierr);
386  this->Flush();
387  }
388 
389 
391 
394  template <class T, class Allocator>
395  inline Vector<T, PETScSeq, Allocator>::
396  Vector(Vec& petsc_vector): PETScVector<T, Allocator>(petsc_vector)
397  {
398  }
399 
400 
402 
405  template <class T, class Allocator>
406  inline Vector<T, PETScSeq, Allocator>::
407  Vector(const Vector<T, PETScSeq, Allocator>& V)
408  {
409  Copy(V);
410  }
411 
412 
413  /**************
414  * DESTRUCTOR *
415  **************/
416 
417 
419  template <class T, class Allocator>
420  inline Vector<T, PETScSeq, Allocator>::~Vector()
421  {
422  }
423 
424 
426 
431  template <class T, class Allocator>
432  inline void Vector<T, PETScSeq, Allocator>
433  ::Copy(const Vector<T, PETScSeq, Allocator>& X)
434  {
435  Copy(X.GetPetscVector());
436  }
437 
438 
440 
445  template <class T, class Allocator>
446  inline void Vector<T, PETScSeq, Allocator>
447  ::Copy(const Vec& X)
448  {
449  PETScVector<T, Allocator>::Copy(X);
450  }
451 
452 
454 
460  template <class T, class Allocator>
461  inline void Vector<T, PETScSeq, Allocator>
462  ::Reallocate(int i)
463  {
464  this->Clear();
465  int ierr;
466  this->m_ = i;
467  ierr = VecCreateSeq(PETSC_COMM_SELF, i, &this->petsc_vector_);
468  CHKERRABORT(this->mpi_communicator_, ierr);
469  Fill(T(0));
470  this->Flush();
471  this->petsc_vector_deallocated_ = false;
472  }
473 
474 
476 
481  template <class T, class Allocator>
482  inline Vector<T, PETScSeq, Allocator>& Vector<T, PETScSeq, Allocator>
483  ::operator= (const Vector<T, PETScSeq, Allocator>& X)
484  {
485  this->Copy(X);
486  return *this;
487  }
488 
489 
491 
494  template <class T, class Allocator>
495  template <class T0>
496  inline Vector<T, PETScSeq, Allocator>&
497  Vector<T, PETScSeq, Allocator>::operator= (const T0& x)
498  {
499  this->Fill(x);
500  return *this;
501  }
502 
503 
505 
508  template <class T, class Allocator>
509  template<class T0>
510  inline Vector<T, PETScSeq, Allocator>& Vector<T, PETScSeq, Allocator>
511  ::operator*= (const T0& alpha)
512  {
513  int ierr;
514  ierr = VecScale(this->petsc_vector_, alpha);
515  CHKERRABORT(this->mpi_communicator_, ierr);
516  return *this;
517  }
518 
519 
521  template <class T, class Allocator>
522  inline void Vector<T, PETScSeq, Allocator>::Print() const
523  {
524  int ierr;
525  ierr = VecView(this->petsc_vector_, PETSC_VIEWER_STDOUT_SELF);
526  }
527 
528 
530  // PARALLEL PETSC VECTOR //
532 
533 
534  /****************
535  * CONSTRUCTORS *
536  ****************/
537 
538 
540 
543  template <class T, class Allocator>
544  inline Vector<T, PETScPar, Allocator>::Vector():
545  PETScVector<T, Allocator>()
546  {
547  this->mpi_communicator_ = MPI_COMM_WORLD;
548  int ierr;
549  ierr = VecCreateMPI(this->mpi_communicator_,
550  PETSC_DECIDE, 0, &this->petsc_vector_);
551  CHKERRABORT(this->mpi_communicator_, ierr);
552  }
553 
554 
556 
560  template <class T, class Allocator>
561  inline Vector<T, PETScPar, Allocator>::Vector(int i, MPI_Comm mpi_communicator):
562  PETScVector<T, Allocator>(i)
563  {
564  int ierr;
565  this->mpi_communicator_ = mpi_communicator;
566  ierr = VecCreateMPI(this->mpi_communicator_, PETSC_DECIDE, i,
567  &this->petsc_vector_);
568  CHKERRABORT(this->mpi_communicator_, ierr);
569  ierr = VecSet(this->petsc_vector_, 0.);
570  CHKERRABORT(this->mpi_communicator_, ierr);
571  this->Flush();
572  }
573 
574 
576 
581  template <class T, class Allocator>
582  inline Vector<T, PETScPar, Allocator>::Vector(int i, int Nlocal,
583  MPI_Comm mpi_communicator):
584  PETScVector<T, Allocator>(i)
585  {
586  int ierr;
587  this->mpi_communicator_ = mpi_communicator;
588  ierr = VecCreateMPI(this->mpi_communicator_, Nlocal,
589  i, &this->petsc_vector_);
590  CHKERRABORT(this->mpi_communicator_, ierr);
591  ierr = VecSet(this->petsc_vector_, 0.);
592  CHKERRABORT(this->mpi_communicator_, ierr);
593  this->Flush();
594  }
595 
596 
598 
601  template <class T, class Allocator>
602  inline Vector<T, PETScPar, Allocator>::
603  Vector(Vec& petsc_vector): PETScVector<T, Allocator>(petsc_vector)
604  {
605  }
606 
607 
609 
612  template <class T, class Allocator>
613  inline Vector<T, PETScPar, Allocator>::
614  Vector(const Vector<T, PETScPar, Allocator>& V)
615  {
616  Copy(V);
617  }
618 
619 
620  /**************
621  * DESTRUCTOR *
622  **************/
623 
624 
626  template <class T, class Allocator>
627  inline Vector<T, PETScPar, Allocator>::~Vector()
628  {
629  }
630 
631 
633 
638  template <class T, class Allocator>
639  inline void Vector<T, PETScPar, Allocator>
640  ::Copy(const Vector<T, PETScPar, Allocator>& X)
641  {
642  Copy(X.GetPetscVector());
643  this->mpi_communicator_ = X.mpi_communicator_;
644  }
645 
646 
648 
653  template <class T, class Allocator>
654  inline void Vector<T, PETScPar, Allocator>
655  ::Copy(const Vec& X)
656  {
657  PETScVector<T, Allocator>::Copy(X);
658  }
659 
660 
662 
667  template <class T, class Allocator>
668  inline Vector<T, PETScPar, Allocator>& Vector<T, PETScPar, Allocator>
669  ::operator= (const Vector<T, PETScPar, Allocator>& X)
670  {
671  this->Copy(X);
672  return *this;
673  }
674 
675 
677 
680  template <class T, class Allocator>
681  template <class T0>
682  Vector<T, PETScPar, Allocator>&
683  Vector<T, PETScPar, Allocator>::operator= (const T0& x)
684  {
685  this->Fill(x);
686  return *this;
687  }
688 
689 
691 
694  template <class T, class Allocator>
695  template<class T0>
696  inline Vector<T, PETScPar, Allocator>& Vector<T, PETScPar, Allocator>
697  ::operator*= (const T0& alpha)
698  {
699  int ierr;
700  ierr = VecScale(this->petsc_vector_, alpha);
701  CHKERRABORT(this->mpi_communicator_, ierr);
702  return *this;
703  }
704 
705 
707  template <class T, class Allocator>
708  inline void Vector<T, PETScPar, Allocator>::Print() const
709  {
710  int ierr;
711  ierr = VecView(this->petsc_vector_, PETSC_VIEWER_STDOUT_WORLD);
712  }
713 
714 
716 
722  template <class T, class Allocator>
723  inline void Vector<T, PETScPar, Allocator>
724  ::Reallocate(int i, int local_size)
725  {
726  this->Clear();
727  int ierr;
728  this->m_ = i;
729  ierr = VecCreateMPI(this->mpi_communicator_, local_size, i,
730  &this->petsc_vector_);
731  CHKERRABORT(this->mpi_communicator_, ierr);
732  Fill(T(0));
733  this->Flush();
734  this->petsc_vector_deallocated_ = false;
735  }
736 
737 
738 } // namespace Seldon.
739 
740 
741 #define SELDON_FILE_PETSCVECTOR_INLINE_CXX
742 #endif
Seldon::PETScVector::GetPetscVector
Vec & GetPetscVector()
Returns a reference on the inner petsc vector.
Definition: PetscVectorInline.cxx:106
Seldon::PETScVector::GetProcessorRange
void GetProcessorRange(int &i, int &j) const
Returns the range of indices owned by this processor.
Definition: PetscVectorInline.cxx:280
Seldon::Vector< T, PETScPar, Allocator >
Definition: PetscVector.hxx:162
Seldon::PETScVector::GetLocalM
int GetLocalM() const
Returns the number of elements stored.
Definition: PetscVectorInline.cxx:336
Seldon::PETScVector::mpi_communicator_
MPI_Comm mpi_communicator_
The MPI communicator to use.
Definition: PetscVector.hxx:53
Seldon::to_str
std::string to_str(const T &input)
Converts most types to string.
Definition: CommonInline.cxx:137
Seldon::Vector
Definition: SeldonHeader.hxx:207
Seldon::WrongIndex
Definition: Errors.hxx:114
Seldon::PETScVector::Copy
void Copy(const PETScVector< T, Allocator > &X)
Duplicates a vector.
Definition: PetscVectorInline.cxx:296
Seldon::PETScVector::PETScVector
PETScVector()
Default constructor.
Definition: PetscVectorInline.cxx:41
Seldon::PETScVector::Flush
void Flush()
Assembles the PETSc vector.
Definition: PetscVectorInline.cxx:259
Seldon::PETScVector::Append
void Append(const T &x)
Appends an element to the vector.
Definition: PetscVectorInline.cxx:309
Seldon::Vector_Base
Base structure for all vectors.
Definition: SeldonHeader.hxx:201
Seldon::PETScVector
PETSc vector class.
Definition: PetscVector.hxx:41
Seldon::PETScVector::~PETScVector
~PETScVector()
Destructor.
Definition: PetscVectorInline.cxx:95
Seldon::Undefined
Definition: Errors.hxx:62
Seldon::PETScVector::operator()
value_type operator()(int i) const
Access operator.
Definition: PetscVectorInline.cxx:216
Seldon::PETScVector::SetCommunicator
void SetCommunicator(MPI_Comm mpi_communicator)
Sets the MPI communicator.
Definition: PetscVectorInline.cxx:128
Seldon::PETScVector::Clear
void Clear()
Clears the vector.
Definition: PetscVectorInline.cxx:145
Seldon::PETScVector::GetDataSize
int GetDataSize() const
Returns the number of elements stored.
Definition: PetscVectorInline.cxx:325
Seldon::PETScVector::Resize
void Resize(int i)
Changes the length of the vector, and keeps previous values.
Definition: PetscVectorInline.cxx:162
Seldon::Vector< T, PETScSeq, Allocator >
Definition: PetscVector.hxx:105
Seldon::PETScVector::SetData
void SetData(int i, pointer data)
Changes the length of the vector and sets its data array (low level method).
Definition: PetscVectorInline.cxx:184
Seldon::PETScVector::SetBuffer
void SetBuffer(int i, T value, InsertMode insert_mode=INSERT_VALUES)
Inserts or adds values into certain locations of a vector.
Definition: PetscVectorInline.cxx:245
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::PETScVector::petsc_vector_deallocated_
bool petsc_vector_deallocated_
Boolean to indicate if the inner PETSc vector is destroyed or not.
Definition: PetscVector.hxx:55
Seldon::PETScVector::Nullify
void Nullify()
Clears the vector without releasing memory.
Definition: PetscVectorInline.cxx:198