PetscMatrixInline.cxx
1 // Copyright (C) 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_MATRIX_PETSCMATRIX_INLINE_CXX
22 
23 #include "PetscMatrix.hxx"
24 
25 
26 namespace Seldon
27 {
28 
29 
31 
34  template <class T, class Prop, class Storage, class Allocator>
36  Matrix_Base<T, Allocator>()
37  {
38  mpi_communicator_ = MPI_COMM_WORLD;
39  MatCreate(PETSC_COMM_WORLD, &petsc_matrix_);
40  petsc_matrix_deallocated_ = false;
41  }
42 
43 
45 
49  template <class T, class Prop, class Storage, class Allocator>
51  Matrix_Base<T, Allocator>(i, j)
52  {
53  int ierr;
54  mpi_communicator_ = MPI_COMM_WORLD;
55  MatCreate(mpi_communicator_, &petsc_matrix_);
56  petsc_matrix_deallocated_ = false;
57  }
58 
59 
61  template <class T, class Prop, class Storage, class Allocator>
63  ::PetscMatrix(Mat& A): Matrix_Base<T, Allocator>()
64  {
65  Copy(A);
66  }
67 
68 
70 
73  template <class T, class Prop, class Storage, class Allocator>
75  ::SetCommunicator(MPI_Comm mpi_communicator)
76  {
77  mpi_communicator_ = mpi_communicator;
78  }
79 
80 
82 
85  template <class T, class Prop, class Storage, class Allocator>
88  {
89  return mpi_communicator_;
90  }
91 
92 
94  template <class T, class Prop, class Storage, class Allocator>
97  {
98  Clear();
99  }
100 
101 
103 
107  template <class T, class Prop, class Storage, class Allocator>
110  {
111  if (petsc_matrix_deallocated_)
112  return;
113  int ierr;
114  ierr = MatDestroy(&petsc_matrix_);
115  CHKERRABORT(mpi_communicator_, ierr);
116  petsc_matrix_deallocated_ = true;
117  }
118 
119 
121 
126  template <class T, class Prop, class Storage, class Allocator>
129  {
130  throw Undefined("PetscMatrix<T, Prop, Storage, Allocator>"
131  "::Nullify()");
132  }
133 
134 
136 
139  template <class T, class Prop, class Storage, class Allocator>
142  {
143  return petsc_matrix_;
144  }
145 
146 
148 
151  template <class T, class Prop, class Storage, class Allocator>
153  ::GetPetscMatrix() const
154  {
155  return petsc_matrix_;
156  }
157 
158 
160  template <class T, class Prop, class Storage, class Allocator>
161  inline size_t PetscMatrix<T, Prop, Storage, Allocator>
162  ::GetMemorySize() const
163  {
164  return 0;
165  }
166 
167 
169 
176  template <class T, class Prop, class Storage, class Allocator>
178  ::Resize(int i, int j)
179  {
180  throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
181  "::Resize(int i, int j)");
182  }
183 
184 
187 
201  template <class T, class Prop, class Storage, class Allocator>
203  ::SetData(int i, int j, pointer data)
204  {
205  throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
206  "::SetData(int i, int j, pointer data)");
207  }
208 
209 
211 
217  template <class T, class Prop, class Storage, class Allocator>
218  inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
220  {
221  throw Undefined("PetscMatrix::Val(int i, int j) const");
222  }
223 
224 
226 
232  template <class T, class Prop, class Storage, class Allocator>
233  inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
235  {
236  throw Undefined("PetscMatrix::Val(int i, int j)");
237  }
238 
239 
241 
246  template <class T, class Prop, class Storage, class Allocator>
247  inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
249  {
250  throw Undefined("PetscMatrix::operator[] (int i)");
251  }
252 
253 
255 
260  template <class T, class Prop, class Storage, class Allocator>
261  inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
263  {
264  throw Undefined("PetscMatrix::operator[] (int i) const");
265  }
266 
267 
268  template <class T, class Prop, class Storage, class Allocator>
269  inline void PetscMatrix<T, Prop, Storage, Allocator>::Set(int i, int j, T value)
270  {
271  throw Undefined("PetscMatrix::Set(int i, int j, T value)");
272  }
273 
274 
276 
285  template <class T, class Prop, class Storage, class Allocator>
286  inline void PetscMatrix<T, Prop, Storage, Allocator>
287  ::SetBuffer(int i, int j, T value, InsertMode insert_mode = INSERT_VALUES)
288  {
289  int ierr;
290  ierr = MatSetValues(petsc_matrix_, 1, &i, 1,
291  &j, &value, insert_mode);
292  CHKERRABORT(mpi_communicator_, ierr);
293  }
294 
295 
297  template <class T, class Prop, class Storage, class Allocator>
299  {
300  int ierr;
301  ierr = MatAssemblyBegin(petsc_matrix_, MAT_FINAL_ASSEMBLY);
302  CHKERRABORT(mpi_communicator_, ierr);
303  ierr = MatAssemblyEnd(petsc_matrix_, MAT_FINAL_ASSEMBLY);
304  CHKERRABORT(mpi_communicator_, ierr);
305  }
306 
307 
309 
317  template <class T, class Prop, class Storage, class Allocator>
319  ::GetProcessorRowRange(int& i, int& j) const
320  {
321  int ierr;
322  ierr = MatGetOwnershipRange(petsc_matrix_, &i, &j);
323  CHKERRABORT(mpi_communicator_, ierr);
324  }
325 
326 
328  // Matrix<PETScSeqDense> //
330 
331 
332  /****************
333  * CONSTRUCTORS *
334  ****************/
335 
336 
338 
341  template <class T, class Prop, class Allocator>
343  PetscMatrix<T, Prop, RowMajor, Allocator>()
344  {
345  }
346 
347 
349 
353  template <class T, class Prop, class Allocator>
355  PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
356  {
357  MatSetType(this->petsc_matrix_, MATSEQDENSE);
358  MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
359  }
360 
361 
363 
366  template <class T, class Prop, class Allocator>
368  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
369  {
370  }
371 
372 
374 
377  template <class T, class Prop, class Allocator>
380  {
381  this->petsc_matrix_deallocated_ = true;
382  Copy(A);
383  }
384 
385 
387 
393  template <class T, class Prop, class Allocator>
395  ::Reallocate(int i, int j)
396  {
397  this->Clear();
398  int ierr;
399  MatCreate(MPI_COMM_SELF, &this->petsc_matrix_);
400  MatSetSizes(this->petsc_matrix_, i, j, i, j);
401  MatSetType(this->petsc_matrix_, MATSEQDENSE);
402  this->petsc_matrix_deallocated_ = false;
403  this->Flush();
404  }
405 
406 
408 
414  template <class T, class Prop, class Allocator>
417  {
418 #ifdef SELDON_CHECK_BOUNDS
419  if (i < 0 || i >= this->m_)
420  throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
421  string("Index should be in [0, ")
422  + to_str(this->m_-1)
423  + "], but is equal to " + to_str(i) + ".");
424  if (j < 0 || j >= this->n_)
425  throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
426  string("Index should be in [0, ") + to_str(this->n_-1)
427  + "], but is equal to " + to_str(j) + ".");
428 #endif
429  PetscInt idxm[1] = {i};
430  PetscInt idxn[1] = {j};
431  PetscScalar v[1];
432  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
433  return v[0];
434  }
435 
436 
438 
444  template <class T, class Prop, class Allocator>
447  {
448 #ifdef SELDON_CHECK_BOUNDS
449  if (i < 0 || i >= this->m_)
450  throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
451  string("Index should be in [0, ")
452  + to_str(this->m_-1)
453  + "], but is equal to " + to_str(i) + ".");
454  if (j < 0 || j >= this->n_)
455  throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
456  string("Index should be in [0, ")
457  + to_str(this->n_-1)
458  + "], but is equal to " + to_str(j) + ".");
459 #endif
460  PetscInt idxm[1] = {i};
461  PetscInt idxn[1] = {j};
462  PetscScalar v[1];
463  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
464  return v[0];
465  }
466 
467 
469 
474  template <class T, class Prop, class Allocator>
477  {
478  this->mpi_communicator_ = A.mpi_communicator_;
480  this->petsc_matrix_deallocated_ = false;
481  }
482 
483 
485 
490  template <class T, class Prop, class Allocator>
494  {
495  this->Copy(A);
496  return *this;
497  }
498 
499 
501  // Matrix<PETScMPIDense> //
503 
504 
505  /****************
506  * CONSTRUCTORS *
507  ****************/
508 
509 
511 
514  template <class T, class Prop, class Allocator>
516  PetscMatrix<T, Prop, RowMajor, Allocator>()
517  {
518  }
519 
520 
522 
526  template <class T, class Prop, class Allocator>
528  PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
529  {
530  MatSetType(this->petsc_matrix_, MATMPIDENSE);
531  MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
532  }
533 
534 
536 
539  template <class T, class Prop, class Allocator>
541  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
542  {
543  }
544 
545 
547 
550  template <class T, class Prop, class Allocator>
553  {
554  this->petsc_matrix_deallocated_ = true;
555  Copy(A);
556  }
557 
558 
560 
566  template <class T, class Prop, class Allocator>
568  ::Reallocate(int i, int j, int local_i, int local_j)
569  {
570  this->Clear();
571  int ierr;
572  MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
573  MatSetType(this->petsc_matrix_, MATMPIDENSE);
574  MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
575  this->m_ = i;
576  this->n_ = j;
577  this->petsc_matrix_deallocated_ = false;
578  this->Flush();
579  }
580 
581 
583 
589  template <class T, class Prop, class Allocator>
592  {
593 #ifdef SELDON_CHECK_BOUNDS
594  int start, end;
595  this->GetProcessorRowRange(start, end);
596  if (i < start || i >= end)
597  throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
598  string("Index should be in [")
599  + to_str(start)
600  + ", " + to_str(end - 1)
601  + "], but is equal to "
602  + to_str(i) + ".");
603  if (j < 0 || j >= this->n_)
604  throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
605  string("Index should be in [0, ") + to_str(this->n_-1)
606  + "], but is equal to " + to_str(j) + ".");
607 #endif
608  PetscInt idxm[1] = {i};
609  PetscInt idxn[1] = {j};
610  PetscScalar v[1];
611  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
612  return v[0];
613  }
614 
615 
617 
623  template <class T, class Prop, class Allocator>
626  {
627 #ifdef SELDON_CHECK_BOUNDS
628  int start, end;
629  this->GetProcessorRowRange(start, end);
630  if (i < start || i >= end)
631  throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
632  string("Index should be in [")
633  + to_str(start)
634  + ", " + to_str(end - 1) + "], but is equal to "
635  + to_str(i) + ".");
636  if (j < 0 || j >= this->n_)
637  throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
638  string("Index should be in [0, ")
639  + to_str(this->n_-1)
640  + "], but is equal to " + to_str(j) + ".");
641 #endif
642  PetscInt idxm[1] = {i};
643  PetscInt idxn[1] = {j};
644  PetscScalar v[1];
645  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
646  return v[0];
647  }
648 
649 
651 
656  template <class T, class Prop, class Allocator>
659  {
660  this->mpi_communicator_ = A.mpi_communicator_;
662  this->petsc_matrix_deallocated_ = false;
663  }
664 
665 
667 
672  template <class T, class Prop, class Allocator>
676  {
677  this->Copy(A);
678  return *this;
679  }
680 
681 
683  // Matrix<PETScMPIAIJ> //
685 
686 
687  /****************
688  * CONSTRUCTORS *
689  ****************/
690 
691 
693 
696  template <class T, class Prop, class Allocator>
698  PetscMatrix<T, Prop, RowMajor, Allocator>()
699  {
700  }
701 
702 
704 
708  template <class T, class Prop, class Allocator>
710  PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
711  {
712  MatSetType(this->petsc_matrix_, MATMPIAIJ);
713  MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
714  }
715 
716 
718 
721  template <class T, class Prop, class Allocator>
723  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
724  {
725  }
726 
727 
729 
732  template <class T, class Prop, class Allocator>
735  {
736  this->petsc_matrix_deallocated_ = true;
737  Copy(A);
738  }
739 
740 
742 
748  template <class T, class Prop, class Allocator>
750  ::Reallocate(int i, int j, int local_i, int local_j)
751  {
752  this->Clear();
753  int ierr;
754  MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
755  MatSetType(this->petsc_matrix_, MATMPIAIJ);
756  MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
757  this->m_ = i;
758  this->n_ = j;
759  this->petsc_matrix_deallocated_ = false;
760  this->Flush();
761  }
762 
763 
765 
771  template <class T, class Prop, class Allocator>
774  {
775 #ifdef SELDON_CHECK_BOUNDS
776  int start, end;
777  this->GetProcessorRowRange(start, end);
778  if (i < start || i >= end)
779  throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
780  string("Index should be in [")
781  + to_str(start)
782  + ", " + to_str(end - 1)
783  + "], but is equal to "
784  + to_str(i) + ".");
785  if (j < 0 || j >= this->n_)
786  throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
787  string("Index should be in [0, ") + to_str(this->n_-1)
788  + "], but is equal to " + to_str(j) + ".");
789 #endif
790  PetscInt idxm[1] = {i};
791  PetscInt idxn[1] = {j};
792  PetscScalar v[1];
793  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
794  return v[0];
795  }
796 
797 
799 
805  template <class T, class Prop, class Allocator>
808  {
809 #ifdef SELDON_CHECK_BOUNDS
810  int start, end;
811  this->GetProcessorRowRange(start, end);
812  if (i < start || i >= end)
813  throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
814  string("Index should be in [")
815  + to_str(start)
816  + ", " + to_str(end - 1) + "], but is equal to "
817  + to_str(i) + ".");
818  if (j < 0 || j >= this->n_)
819  throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
820  string("Index should be in [0, ")
821  + to_str(this->n_-1)
822  + "], but is equal to " + to_str(j) + ".");
823 #endif
824  PetscInt idxm[1] = {i};
825  PetscInt idxn[1] = {j};
826  PetscScalar v[1];
827  MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
828  return v[0];
829  }
830 
831 
833 
838  template <class T, class Prop, class Allocator>
841  {
842  this->mpi_communicator_ = A.mpi_communicator_;
844  this->petsc_matrix_deallocated_ = false;
845  }
846 
847 
849 
854  template <class T, class Prop, class Allocator>
858  {
859  this->Copy(A);
860  return *this;
861  }
862 
863 }
864 
865 
866 #define SELDON_FILE_MATRIX_PETSCMATRIX_INLINE_CXX
867 #endif
Seldon::RowMajor
Definition: Storage.hxx:45
Seldon::Matrix< T, Prop, PETScMPIDense, Allocator >
Definition: PetscMatrix.hxx:154
Seldon::Matrix_Base
Base class for all matrices.
Definition: Matrix_Base.hxx:143
Seldon::to_str
std::string to_str(const T &input)
Converts most types to string.
Definition: CommonInline.cxx:137
Seldon::WrongCol
Definition: Errors.hxx:138
Seldon::Matrix< T, Prop, PETScSeqDense, Allocator >
Definition: PetscMatrix.hxx:121
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::PetscMatrix::GetPetscMatrix
Mat & GetPetscMatrix()
Returns a reference on the inner petsc matrix.
Definition: PetscMatrixInline.cxx:141
Seldon::WrongRow
Definition: Errors.hxx:126
Seldon::PetscMatrix
Matrix class based on PETSc matrix.
Definition: PetscMatrix.hxx:40
Seldon::Matrix< T, Prop, PETScMPIAIJ, Allocator >
Definition: PetscMatrix.hxx:188
Seldon::Undefined
Definition: Errors.hxx:62
Seldon::PetscMatrix::mpi_communicator_
MPI_Comm mpi_communicator_
The MPI communicator to use.
Definition: PetscMatrix.hxx:58
Seldon::PetscMatrix::PetscMatrix
PetscMatrix()
Default constructor.
Definition: PetscMatrixInline.cxx:35
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::PetscMatrix::Copy
void Copy(const Mat &A)
Duplicates a matrix.
Definition: PetscMatrix.cxx:37