21 #ifndef SELDON_FILE_MATRIX_PETSCMATRIX_INLINE_CXX
23 #include "PetscMatrix.hxx"
34 template <
class T,
class Prop,
class Storage,
class Allocator>
38 mpi_communicator_ = MPI_COMM_WORLD;
39 MatCreate(PETSC_COMM_WORLD, &petsc_matrix_);
40 petsc_matrix_deallocated_ =
false;
49 template <
class T,
class Prop,
class Storage,
class Allocator>
54 mpi_communicator_ = MPI_COMM_WORLD;
55 MatCreate(mpi_communicator_, &petsc_matrix_);
56 petsc_matrix_deallocated_ =
false;
61 template <
class T,
class Prop,
class Storage,
class Allocator>
73 template <
class T,
class Prop,
class Storage,
class Allocator>
77 mpi_communicator_ = mpi_communicator;
85 template <
class T,
class Prop,
class Storage,
class Allocator>
89 return mpi_communicator_;
94 template <
class T,
class Prop,
class Storage,
class Allocator>
107 template <
class T,
class Prop,
class Storage,
class Allocator>
111 if (petsc_matrix_deallocated_)
114 ierr = MatDestroy(&petsc_matrix_);
115 CHKERRABORT(mpi_communicator_, ierr);
116 petsc_matrix_deallocated_ =
true;
126 template <
class T,
class Prop,
class Storage,
class Allocator>
130 throw Undefined(
"PetscMatrix<T, Prop, Storage, Allocator>"
139 template <
class T,
class Prop,
class Storage,
class Allocator>
143 return petsc_matrix_;
151 template <
class T,
class Prop,
class Storage,
class Allocator>
155 return petsc_matrix_;
160 template <
class T,
class Prop,
class Storage,
class Allocator>
161 inline size_t PetscMatrix<T, Prop, Storage, Allocator>
162 ::GetMemorySize()
const
176 template <
class T,
class Prop,
class Storage,
class Allocator>
180 throw Undefined(
"void PetscMatrix<T, Prop, Storage, Allocator>"
181 "::Resize(int i, int j)");
201 template <
class T,
class Prop,
class Storage,
class Allocator>
205 throw Undefined(
"void PetscMatrix<T, Prop, Storage, Allocator>"
206 "::SetData(int i, int j, pointer data)");
217 template <
class T,
class Prop,
class Storage,
class Allocator>
218 inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
221 throw Undefined(
"PetscMatrix::Val(int i, int j) const");
232 template <
class T,
class Prop,
class Storage,
class Allocator>
233 inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
236 throw Undefined(
"PetscMatrix::Val(int i, int j)");
246 template <
class T,
class Prop,
class Storage,
class Allocator>
247 inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
250 throw Undefined(
"PetscMatrix::operator[] (int i)");
260 template <
class T,
class Prop,
class Storage,
class Allocator>
261 inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
264 throw Undefined(
"PetscMatrix::operator[] (int i) const");
268 template <
class T,
class Prop,
class Storage,
class Allocator>
271 throw Undefined(
"PetscMatrix::Set(int i, int j, T value)");
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)
290 ierr = MatSetValues(petsc_matrix_, 1, &i, 1,
291 &j, &value, insert_mode);
292 CHKERRABORT(mpi_communicator_, ierr);
297 template <
class T,
class Prop,
class Storage,
class Allocator>
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);
317 template <
class T,
class Prop,
class Storage,
class Allocator>
322 ierr = MatGetOwnershipRange(petsc_matrix_, &i, &j);
323 CHKERRABORT(mpi_communicator_, ierr);
341 template <
class T,
class Prop,
class Allocator>
353 template <
class T,
class Prop,
class Allocator>
357 MatSetType(this->petsc_matrix_, MATSEQDENSE);
358 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
366 template <
class T,
class Prop,
class Allocator>
377 template <
class T,
class Prop,
class Allocator>
381 this->petsc_matrix_deallocated_ =
true;
393 template <
class T,
class Prop,
class Allocator>
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;
414 template <
class T,
class Prop,
class Allocator>
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, ")
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) +
".");
429 PetscInt idxm[1] = {i};
430 PetscInt idxn[1] = {j};
432 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
444 template <
class T,
class Prop,
class Allocator>
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, ")
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, ")
458 +
"], but is equal to " +
to_str(j) +
".");
460 PetscInt idxm[1] = {i};
461 PetscInt idxn[1] = {j};
463 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
474 template <
class T,
class Prop,
class Allocator>
480 this->petsc_matrix_deallocated_ =
false;
490 template <
class T,
class Prop,
class Allocator>
514 template <
class T,
class Prop,
class Allocator>
526 template <
class T,
class Prop,
class Allocator>
530 MatSetType(this->petsc_matrix_, MATMPIDENSE);
531 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
539 template <
class T,
class Prop,
class Allocator>
550 template <
class T,
class Prop,
class Allocator>
554 this->petsc_matrix_deallocated_ =
true;
566 template <
class T,
class Prop,
class Allocator>
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);
577 this->petsc_matrix_deallocated_ =
false;
589 template <
class T,
class Prop,
class Allocator>
593 #ifdef SELDON_CHECK_BOUNDS
595 this->GetProcessorRowRange(start, end);
596 if (i < start || i >= end)
597 throw WrongRow(
"PetscMatrix<PETScMPIDense>::operator()",
598 string(
"Index should be in [")
601 +
"], but is equal to "
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) +
".");
608 PetscInt idxm[1] = {i};
609 PetscInt idxn[1] = {j};
611 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
623 template <
class T,
class Prop,
class Allocator>
627 #ifdef SELDON_CHECK_BOUNDS
629 this->GetProcessorRowRange(start, end);
630 if (i < start || i >= end)
631 throw WrongRow(
"PetscMatrix<PETScMPIDense>::operator()",
632 string(
"Index should be in [")
634 +
", " +
to_str(end - 1) +
"], but is equal to "
636 if (j < 0 || j >= this->n_)
637 throw WrongCol(
"PetscMatrix<PETScMPIDense>::operator()",
638 string(
"Index should be in [0, ")
640 +
"], but is equal to " +
to_str(j) +
".");
642 PetscInt idxm[1] = {i};
643 PetscInt idxn[1] = {j};
645 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
656 template <
class T,
class Prop,
class Allocator>
662 this->petsc_matrix_deallocated_ =
false;
672 template <
class T,
class Prop,
class Allocator>
696 template <
class T,
class Prop,
class Allocator>
708 template <
class T,
class Prop,
class Allocator>
712 MatSetType(this->petsc_matrix_, MATMPIAIJ);
713 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
721 template <
class T,
class Prop,
class Allocator>
732 template <
class T,
class Prop,
class Allocator>
736 this->petsc_matrix_deallocated_ =
true;
748 template <
class T,
class Prop,
class Allocator>
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);
759 this->petsc_matrix_deallocated_ =
false;
771 template <
class T,
class Prop,
class Allocator>
775 #ifdef SELDON_CHECK_BOUNDS
777 this->GetProcessorRowRange(start, end);
778 if (i < start || i >= end)
779 throw WrongRow(
"PetscMatrix<PETScMPIAIJ>::operator()",
780 string(
"Index should be in [")
783 +
"], but is equal to "
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) +
".");
790 PetscInt idxm[1] = {i};
791 PetscInt idxn[1] = {j};
793 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
805 template <
class T,
class Prop,
class Allocator>
809 #ifdef SELDON_CHECK_BOUNDS
811 this->GetProcessorRowRange(start, end);
812 if (i < start || i >= end)
813 throw WrongRow(
"PetscMatrix<PETScMPIAIJ>::operator()",
814 string(
"Index should be in [")
816 +
", " +
to_str(end - 1) +
"], but is equal to "
818 if (j < 0 || j >= this->n_)
819 throw WrongCol(
"PetscMatrix<PETScMPIAIJ>::operator()",
820 string(
"Index should be in [0, ")
822 +
"], but is equal to " +
to_str(j) +
".");
824 PetscInt idxm[1] = {i};
825 PetscInt idxn[1] = {j};
827 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
838 template <
class T,
class Prop,
class Allocator>
844 this->petsc_matrix_deallocated_ =
false;
854 template <
class T,
class Prop,
class Allocator>
866 #define SELDON_FILE_MATRIX_PETSCMATRIX_INLINE_CXX