00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef SELDON_FILE_MATRIX_PETSCMATRIX_CXX
00022
00023 #include "PetscMatrix.hxx"
00024
00025
00026 namespace Seldon
00027 {
00028
00029
00031
00034 template <class T, class Prop, class Storage, class Allocator>
00035 inline PetscMatrix<T, Prop, Storage, Allocator>::PetscMatrix():
00036 Matrix_Base<T, Allocator>()
00037 {
00038 mpi_communicator_ = MPI_COMM_WORLD;
00039 MatCreate(PETSC_COMM_WORLD, &petsc_matrix_);
00040 petsc_matrix_deallocated_ = false;
00041 }
00042
00043
00045
00049 template <class T, class Prop, class Storage, class Allocator>
00050 inline PetscMatrix<T, Prop, Storage, Allocator>::PetscMatrix(int i, int j):
00051 Matrix_Base<T, Allocator>(i, j)
00052 {
00053 int ierr;
00054 mpi_communicator_ = MPI_COMM_WORLD;
00055 MatCreate(mpi_communicator_, &petsc_matrix_);
00056 petsc_matrix_deallocated_ = false;
00057 }
00058
00059
00061 template <class T, class Prop, class Storage, class Allocator>
00062 inline PetscMatrix<T, Prop, Storage, Allocator>
00063 ::PetscMatrix(Mat& A): Matrix_Base<T, Allocator>()
00064 {
00065 Copy(A);
00066 }
00067
00068
00070
00073 template <class T, class Prop, class Storage, class Allocator>
00074 inline void PetscMatrix<T, Prop, Storage, Allocator>
00075 ::SetCommunicator(MPI_Comm mpi_communicator)
00076 {
00077 mpi_communicator_ = mpi_communicator;
00078 }
00079
00080
00082
00085 template <class T, class Prop, class Storage, class Allocator>
00086 inline MPI_Comm PetscMatrix<T, Prop, Storage, Allocator>
00087 ::GetCommunicator() const
00088 {
00089 return mpi_communicator_;
00090 }
00091
00092
00094 template <class T, class Prop, class Storage, class Allocator>
00095 inline PetscMatrix<T, Prop, Storage, Allocator>
00096 ::~PetscMatrix()
00097 {
00098 Clear();
00099 }
00100
00101
00103
00107 template <class T, class Prop, class Storage, class Allocator>
00108 inline void PetscMatrix<T, Prop, Storage, Allocator>
00109 ::Clear()
00110 {
00111 if (petsc_matrix_deallocated_)
00112 return;
00113 int ierr;
00114 ierr = MatDestroy(&petsc_matrix_);
00115 CHKERRABORT(mpi_communicator_, ierr);
00116 petsc_matrix_deallocated_ = true;
00117 }
00118
00119
00121
00126 template <class T, class Prop, class Storage, class Allocator>
00127 inline void PetscMatrix<T, Prop, Storage, Allocator>
00128 ::Nullify()
00129 {
00130 throw Undefined("PetscMatrix<T, Prop, Storage, Allocator>"
00131 "::Nullify()");
00132 }
00133
00134
00136
00139 template <class T, class Prop, class Storage, class Allocator>
00140 inline Mat& PetscMatrix<T, Prop, Storage, Allocator>
00141 ::GetPetscMatrix()
00142 {
00143 return petsc_matrix_;
00144 }
00145
00146
00148
00151 template <class T, class Prop, class Storage, class Allocator>
00152 inline const Mat& PetscMatrix<T, Prop, Storage, Allocator>
00153 ::GetPetscMatrix() const
00154 {
00155 return petsc_matrix_;
00156 }
00157
00158
00160
00167 template <class T, class Prop, class Storage, class Allocator>
00168 inline void PetscMatrix<T, Prop, Storage, Allocator>
00169 ::Resize(int i, int j)
00170 {
00171 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00172 "::Resize(int i, int j)");
00173 }
00174
00175
00178
00192 template <class T, class Prop, class Storage, class Allocator>
00193 inline void PetscMatrix<T, Prop, Storage, Allocator>
00194 ::SetData(int i, int j, pointer data)
00195 {
00196 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00197 "::SetData(int i, int j, pointer data)");
00198 }
00199
00200
00202
00208 template <class T, class Prop, class Storage, class Allocator>
00209 inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
00210 PetscMatrix<T, Prop, Storage, Allocator>::Val(int i, int j) const
00211 {
00212 throw Undefined("PetscMatrix::Val(int i, int j) const");
00213 }
00214
00215
00217
00223 template <class T, class Prop, class Storage, class Allocator>
00224 inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
00225 PetscMatrix<T, Prop, Storage, Allocator>::Val(int i, int j)
00226 {
00227 throw Undefined("PetscMatrix::Val(int i, int j)");
00228 }
00229
00230
00232
00237 template <class T, class Prop, class Storage, class Allocator>
00238 inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
00239 PetscMatrix<T, Prop, Storage, Allocator>::operator[] (int i)
00240 {
00241 throw Undefined("PetscMatrix::operator[] (int i)");
00242 }
00243
00244
00246
00251 template <class T, class Prop, class Storage, class Allocator>
00252 inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
00253 PetscMatrix<T, Prop, Storage, Allocator>::operator[] (int i) const
00254 {
00255 throw Undefined("PetscMatrix::operator[] (int i) const");
00256 }
00257
00258
00259 template <class T, class Prop, class Storage, class Allocator>
00260 void PetscMatrix<T, Prop, Storage, Allocator>::Set(int i, int j, T value)
00261 {
00262 throw Undefined("PetscMatrix::Set(int i, int j, T value)");
00263 }
00264
00265
00267
00276 template <class T, class Prop, class Storage, class Allocator>
00277 void PetscMatrix<T, Prop, Storage, Allocator>
00278 ::SetBuffer(int i, int j, T value, InsertMode insert_mode = INSERT_VALUES)
00279 {
00280 int ierr;
00281 ierr = MatSetValues(petsc_matrix_, 1, &i, 1,
00282 &j, &value, insert_mode);
00283 CHKERRABORT(mpi_communicator_, ierr);
00284 }
00285
00286
00288 template <class T, class Prop, class Storage, class Allocator>
00289 void PetscMatrix<T, Prop, Storage, Allocator>::Flush() const
00290 {
00291 int ierr;
00292 ierr = MatAssemblyBegin(petsc_matrix_, MAT_FINAL_ASSEMBLY);
00293 CHKERRABORT(mpi_communicator_, ierr);
00294 ierr = MatAssemblyEnd(petsc_matrix_, MAT_FINAL_ASSEMBLY);
00295 CHKERRABORT(mpi_communicator_, ierr);
00296 }
00297
00298
00300
00308 template <class T, class Prop, class Storage, class Allocator>
00309 void PetscMatrix<T, Prop, Storage, Allocator>
00310 ::GetProcessorRowRange(int& i, int& j) const
00311 {
00312 int ierr;
00313 ierr = MatGetOwnershipRange(petsc_matrix_, &i, &j);
00314 CHKERRABORT(mpi_communicator_, ierr);
00315 }
00316
00317
00319
00324 template <class T, class Prop, class Storage, class Allocator>
00325 inline void PetscMatrix<T, Prop, Storage, Allocator>
00326 ::Copy(const Mat& A)
00327 {
00328 Clear();
00329 int ierr;
00330 ierr = MatDuplicate(A, MAT_COPY_VALUES, &petsc_matrix_);
00331 CHKERRABORT(mpi_communicator_, ierr);
00332 MatGetSize(A, &this->m_, &this->n_);
00333 mpi_communicator_ = MPI_COMM_WORLD;
00334 petsc_matrix_deallocated_ = false;
00335 }
00336
00337
00339
00343 template <class T, class Prop, class Storage, class Allocator>
00344 void PetscMatrix<T, Prop, Storage, Allocator>::Zero()
00345 {
00346 MatScale(petsc_matrix_, T(0));
00347 }
00348
00349
00351 template <class T, class Prop, class Storage, class Allocator>
00352 void PetscMatrix<T, Prop, Storage, Allocator>::SetIdentity()
00353 {
00354 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00355 "::SetIdentity()");
00356 }
00357
00358
00360
00364 template <class T, class Prop, class Storage, class Allocator>
00365 void PetscMatrix<T, Prop, Storage, Allocator>::Fill()
00366 {
00367 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00368 "::Fill()");
00369 }
00370
00371
00373
00376 template <class T, class Prop, class Storage, class Allocator>
00377 template <class T0>
00378 void PetscMatrix<T, Prop, Storage, Allocator>::Fill(const T0& x)
00379 {
00380 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00381 "::Fill(const T0& x)");
00382 }
00383
00384
00386
00389 template <class T, class Prop, class Storage, class Allocator>
00390 void PetscMatrix<T, Prop, Storage, Allocator>::FillRand()
00391 {
00392 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00393 "::FillRand()");
00394 }
00395
00396
00398
00409 template <class T, class Prop, class Storage, class Allocator>
00410 void PetscMatrix<T, Prop, Storage, Allocator>
00411 ::Print(int a, int b, int m, int n) const
00412 {
00413 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00414 "::Print(int a, int b, int m, int n) const");
00415 }
00416
00417
00419
00427 template <class T, class Prop, class Storage, class Allocator>
00428 void PetscMatrix<T, Prop, Storage, Allocator>::Print(int l) const
00429 {
00430 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00431 "::Print(int l) const");
00432 }
00433
00434
00435
00436
00437
00438
00439
00441
00448 template <class T, class Prop, class Storage, class Allocator>
00449 void PetscMatrix<T, Prop, Storage, Allocator>
00450 ::Write(string FileName, bool with_size) const
00451 {
00452 ofstream FileStream;
00453 FileStream.open(FileName.c_str());
00454
00455 #ifdef SELDON_CHECK_IO
00456
00457 if (!FileStream.is_open())
00458 throw IOError("PetscMatrix::Write(string FileName)",
00459 string("Unable to open file \"") + FileName + "\".");
00460 #endif
00461
00462 this->Write(FileStream, with_size);
00463
00464 FileStream.close();
00465 }
00466
00467
00469
00476 template <class T, class Prop, class Storage, class Allocator>
00477 void PetscMatrix<T, Prop, Storage, Allocator>
00478 ::Write(ostream& FileStream, bool with_size) const
00479 {
00480 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00481 "::Write(ostream& FileStream, bool with_size) const");
00482 }
00483
00484
00486
00493 template <class T, class Prop, class Storage, class Allocator>
00494 void PetscMatrix<T, Prop, Storage, Allocator>
00495 ::WriteText(string FileName) const
00496 {
00497 ofstream FileStream;
00498 FileStream.precision(cout.precision());
00499 FileStream.flags(cout.flags());
00500 FileStream.open(FileName.c_str());
00501
00502 #ifdef SELDON_CHECK_IO
00503
00504 if (!FileStream.is_open())
00505 throw IOError("PetscMatrix::WriteText(string FileName)",
00506 string("Unable to open file \"") + FileName + "\".");
00507 #endif
00508
00509 this->WriteText(FileStream);
00510
00511 FileStream.close();
00512 }
00513
00514
00516
00523 template <class T, class Prop, class Storage, class Allocator>
00524 void PetscMatrix<T, Prop, Storage, Allocator>
00525 ::WriteText(ostream& FileStream) const
00526 {
00527 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00528 "::WriteText(ostream& FileStream) const");
00529 }
00530
00531
00533
00540 template <class T, class Prop, class Storage, class Allocator>
00541 void PetscMatrix<T, Prop, Storage, Allocator>::Read(string FileName)
00542 {
00543 ifstream FileStream;
00544 FileStream.open(FileName.c_str());
00545
00546 #ifdef SELDON_CHECK_IO
00547
00548 if (!FileStream.is_open())
00549 throw IOError("PetscMatrix::Read(string FileName)",
00550 string("Unable to open file \"") + FileName + "\".");
00551 #endif
00552
00553 this->Read(FileStream);
00554
00555 FileStream.close();
00556 }
00557
00558
00560
00567 template <class T, class Prop, class Storage, class Allocator>
00568 void PetscMatrix<T, Prop, Storage, Allocator>
00569 ::Read(istream& FileStream)
00570 {
00571 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00572 "::Read(istream& FileStream) const");
00573 }
00574
00575
00577
00581 template <class T, class Prop, class Storage, class Allocator>
00582 void PetscMatrix<T, Prop, Storage, Allocator>::ReadText(string FileName)
00583 {
00584 ifstream FileStream;
00585 FileStream.open(FileName.c_str());
00586
00587 #ifdef SELDON_CHECK_IO
00588
00589 if (!FileStream.is_open())
00590 throw IOError("Matrix_Pointers::ReadText(string FileName)",
00591 string("Unable to open file \"") + FileName + "\".");
00592 #endif
00593
00594 this->ReadText(FileStream);
00595
00596 FileStream.close();
00597 }
00598
00599
00601
00605 template <class T, class Prop, class Storage, class Allocator>
00606 void PetscMatrix<T, Prop, Storage, Allocator>
00607 ::ReadText(istream& FileStream)
00608 {
00609 throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
00610 "::ReadText(istream& FileStream)");
00611 }
00612
00613
00615
00617
00618
00619
00620
00621
00622
00623
00625
00628 template <class T, class Prop, class Allocator>
00629 inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix():
00630 PetscMatrix<T, Prop, RowMajor, Allocator>()
00631 {
00632 }
00633
00634
00636
00640 template <class T, class Prop, class Allocator>
00641 inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix(int i, int j):
00642 PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
00643 {
00644 MatSetType(this->petsc_matrix_, MATSEQDENSE);
00645 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
00646 }
00647
00648
00650
00653 template <class T, class Prop, class Allocator>
00654 inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix(Mat& A):
00655 PetscMatrix<T, Prop, RowMajor, Allocator>(A)
00656 {
00657 }
00658
00659
00661
00664 template <class T, class Prop, class Allocator>
00665 inline Matrix<T, Prop, PETScSeqDense, Allocator>
00666 ::Matrix(const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
00667 {
00668 this->petsc_matrix_deallocated_ = true;
00669 Copy(A);
00670 }
00671
00672
00674
00680 template <class T, class Prop, class Allocator>
00681 inline void Matrix<T, Prop, PETScSeqDense, Allocator>
00682 ::Reallocate(int i, int j)
00683 {
00684 this->Clear();
00685 int ierr;
00686 MatCreate(MPI_COMM_SELF, &this->petsc_matrix_);
00687 MatSetSizes(this->petsc_matrix_, i, j, i, j);
00688 MatSetType(this->petsc_matrix_, MATSEQDENSE);
00689 this->petsc_matrix_deallocated_ = false;
00690 this->Flush();
00691 }
00692
00693
00695
00701 template <class T, class Prop, class Allocator>
00702 inline typename Matrix<T, Prop, PETScSeqDense, Allocator>::value_type
00703 Matrix<T, Prop, PETScSeqDense, Allocator>::operator() (int i, int j)
00704 {
00705 #ifdef SELDON_CHECK_BOUNDS
00706 if (i < 0 || i >= this->m_)
00707 throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
00708 string("Index should be in [0, ")
00709 + to_str(this->m_-1)
00710 + "], but is equal to " + to_str(i) + ".");
00711 if (j < 0 || j >= this->n_)
00712 throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
00713 string("Index should be in [0, ") + to_str(this->n_-1)
00714 + "], but is equal to " + to_str(j) + ".");
00715 #endif
00716 PetscInt idxm[1] = {i};
00717 PetscInt idxn[1] = {j};
00718 PetscScalar v[1];
00719 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
00720 return v[0];
00721 }
00722
00723
00725
00731 template <class T, class Prop, class Allocator>
00732 inline typename Matrix<T, Prop, PETScSeqDense, Allocator>::value_type
00733 Matrix<T, Prop, PETScSeqDense, Allocator>::operator() (int i, int j) const
00734 {
00735 #ifdef SELDON_CHECK_BOUNDS
00736 if (i < 0 || i >= this->m_)
00737 throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
00738 string("Index should be in [0, ")
00739 + to_str(this->m_-1)
00740 + "], but is equal to " + to_str(i) + ".");
00741 if (j < 0 || j >= this->n_)
00742 throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
00743 string("Index should be in [0, ")
00744 + to_str(this->n_-1)
00745 + "], but is equal to " + to_str(j) + ".");
00746 #endif
00747 PetscInt idxm[1] = {i};
00748 PetscInt idxn[1] = {j};
00749 PetscScalar v[1];
00750 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
00751 return v[0];
00752 }
00753
00754
00756
00761 template <class T, class Prop, class Allocator>
00762 inline void Matrix<T, Prop, PETScSeqDense, Allocator>
00763 ::Copy(const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
00764 {
00765 this->mpi_communicator_ = A.mpi_communicator_;
00766 PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
00767 this->petsc_matrix_deallocated_ = false;
00768 }
00769
00770
00772
00777 template <class T, class Prop, class Allocator>
00778 inline Matrix<T, Prop, PETScSeqDense, Allocator>&
00779 Matrix<T, Prop, PETScSeqDense, Allocator>
00780 ::operator= (const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
00781 {
00782 this->Copy(A);
00783 return *this;
00784 }
00785
00786
00788
00793 template <class T, class Prop, class Allocator>
00794 void Matrix<T, Prop, PETScSeqDense, Allocator>::Print() const
00795 {
00796 int ierr;
00797 ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_SELF);
00798 CHKERRABORT(this->mpi_communicator_, ierr);
00799 }
00800
00801
00803
00805
00806
00807
00808
00809
00810
00811
00813
00816 template <class T, class Prop, class Allocator>
00817 inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix():
00818 PetscMatrix<T, Prop, RowMajor, Allocator>()
00819 {
00820 }
00821
00822
00824
00828 template <class T, class Prop, class Allocator>
00829 inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix(int i, int j):
00830 PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
00831 {
00832 MatSetType(this->petsc_matrix_, MATMPIDENSE);
00833 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
00834 }
00835
00836
00838
00841 template <class T, class Prop, class Allocator>
00842 inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix(Mat& A):
00843 PetscMatrix<T, Prop, RowMajor, Allocator>(A)
00844 {
00845 }
00846
00847
00849
00852 template <class T, class Prop, class Allocator>
00853 inline Matrix<T, Prop, PETScMPIDense, Allocator>
00854 ::Matrix(const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
00855 {
00856 this->petsc_matrix_deallocated_ = true;
00857 Copy(A);
00858 }
00859
00860
00862
00868 template <class T, class Prop, class Allocator>
00869 inline void Matrix<T, Prop, PETScMPIDense, Allocator>
00870 ::Reallocate(int i, int j, int local_i, int local_j)
00871 {
00872 this->Clear();
00873 int ierr;
00874 MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
00875 MatSetType(this->petsc_matrix_, MATMPIDENSE);
00876 MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
00877 this->m_ = i;
00878 this->n_ = j;
00879 this->petsc_matrix_deallocated_ = false;
00880 this->Flush();
00881 }
00882
00883
00885
00891 template <class T, class Prop, class Allocator>
00892 inline typename Matrix<T, Prop, PETScMPIDense, Allocator>::value_type
00893 Matrix<T, Prop, PETScMPIDense, Allocator>::operator() (int i, int j)
00894 {
00895 #ifdef SELDON_CHECK_BOUNDS
00896 int start, end;
00897 this->GetProcessorRowRange(start, end);
00898 if (i < start || i >= end)
00899 throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
00900 string("Index should be in [")
00901 + to_str(start)
00902 + ", " + to_str(end - 1)
00903 + "], but is equal to "
00904 + to_str(i) + ".");
00905 if (j < 0 || j >= this->n_)
00906 throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
00907 string("Index should be in [0, ") + to_str(this->n_-1)
00908 + "], but is equal to " + to_str(j) + ".");
00909 #endif
00910 PetscInt idxm[1] = {i};
00911 PetscInt idxn[1] = {j};
00912 PetscScalar v[1];
00913 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
00914 return v[0];
00915 }
00916
00917
00919
00925 template <class T, class Prop, class Allocator>
00926 inline typename Matrix<T, Prop, PETScMPIDense, Allocator>::value_type
00927 Matrix<T, Prop, PETScMPIDense, Allocator>::operator() (int i, int j) const
00928 {
00929 #ifdef SELDON_CHECK_BOUNDS
00930 int start, end;
00931 this->GetProcessorRowRange(start, end);
00932 if (i < start || i >= end)
00933 throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
00934 string("Index should be in [")
00935 + to_str(start)
00936 + ", " + to_str(end - 1) + "], but is equal to "
00937 + to_str(i) + ".");
00938 if (j < 0 || j >= this->n_)
00939 throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
00940 string("Index should be in [0, ")
00941 + to_str(this->n_-1)
00942 + "], but is equal to " + to_str(j) + ".");
00943 #endif
00944 PetscInt idxm[1] = {i};
00945 PetscInt idxn[1] = {j};
00946 PetscScalar v[1];
00947 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
00948 return v[0];
00949 }
00950
00951
00953
00958 template <class T, class Prop, class Allocator>
00959 inline void Matrix<T, Prop, PETScMPIDense, Allocator>
00960 ::Copy(const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
00961 {
00962 this->mpi_communicator_ = A.mpi_communicator_;
00963 PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
00964 this->petsc_matrix_deallocated_ = false;
00965 }
00966
00967
00969
00974 template <class T, class Prop, class Allocator>
00975 inline Matrix<T, Prop, PETScMPIDense, Allocator>&
00976 Matrix<T, Prop, PETScMPIDense, Allocator>
00977 ::operator= (const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
00978 {
00979 this->Copy(A);
00980 return *this;
00981 }
00982
00983
00985
00990 template <class T, class Prop, class Allocator>
00991 void Matrix<T, Prop, PETScMPIDense, Allocator>::Print() const
00992 {
00993 int ierr;
00994 ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_WORLD);
00995 CHKERRABORT(this->mpi_communicator_, ierr);
00996 }
00997
00998
01000
01002
01003
01004
01005
01006
01007
01008
01010
01013 template <class T, class Prop, class Allocator>
01014 inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix():
01015 PetscMatrix<T, Prop, RowMajor, Allocator>()
01016 {
01017 }
01018
01019
01021
01025 template <class T, class Prop, class Allocator>
01026 inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix(int i, int j):
01027 PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
01028 {
01029 MatSetType(this->petsc_matrix_, MATMPIAIJ);
01030 MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
01031 }
01032
01033
01035
01038 template <class T, class Prop, class Allocator>
01039 inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix(Mat& A):
01040 PetscMatrix<T, Prop, RowMajor, Allocator>(A)
01041 {
01042 }
01043
01044
01046
01049 template <class T, class Prop, class Allocator>
01050 inline Matrix<T, Prop, PETScMPIAIJ, Allocator>
01051 ::Matrix(const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
01052 {
01053 this->petsc_matrix_deallocated_ = true;
01054 Copy(A);
01055 }
01056
01057
01059
01065 template <class T, class Prop, class Allocator>
01066 inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
01067 ::Reallocate(int i, int j, int local_i, int local_j)
01068 {
01069 this->Clear();
01070 int ierr;
01071 MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
01072 MatSetType(this->petsc_matrix_, MATMPIAIJ);
01073 MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
01074 this->m_ = i;
01075 this->n_ = j;
01076 this->petsc_matrix_deallocated_ = false;
01077 this->Flush();
01078 }
01079
01080
01082
01088 template <class T, class Prop, class Allocator>
01089 inline typename Matrix<T, Prop, PETScMPIAIJ, Allocator>::value_type
01090 Matrix<T, Prop, PETScMPIAIJ, Allocator>::operator() (int i, int j)
01091 {
01092 #ifdef SELDON_CHECK_BOUNDS
01093 int start, end;
01094 this->GetProcessorRowRange(start, end);
01095 if (i < start || i >= end)
01096 throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
01097 string("Index should be in [")
01098 + to_str(start)
01099 + ", " + to_str(end - 1)
01100 + "], but is equal to "
01101 + to_str(i) + ".");
01102 if (j < 0 || j >= this->n_)
01103 throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
01104 string("Index should be in [0, ") + to_str(this->n_-1)
01105 + "], but is equal to " + to_str(j) + ".");
01106 #endif
01107 PetscInt idxm[1] = {i};
01108 PetscInt idxn[1] = {j};
01109 PetscScalar v[1];
01110 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
01111 return v[0];
01112 }
01113
01114
01116
01122 template <class T, class Prop, class Allocator>
01123 inline typename Matrix<T, Prop, PETScMPIAIJ, Allocator>::value_type
01124 Matrix<T, Prop, PETScMPIAIJ, Allocator>::operator() (int i, int j) const
01125 {
01126 #ifdef SELDON_CHECK_BOUNDS
01127 int start, end;
01128 this->GetProcessorRowRange(start, end);
01129 if (i < start || i >= end)
01130 throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
01131 string("Index should be in [")
01132 + to_str(start)
01133 + ", " + to_str(end - 1) + "], but is equal to "
01134 + to_str(i) + ".");
01135 if (j < 0 || j >= this->n_)
01136 throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
01137 string("Index should be in [0, ")
01138 + to_str(this->n_-1)
01139 + "], but is equal to " + to_str(j) + ".");
01140 #endif
01141 PetscInt idxm[1] = {i};
01142 PetscInt idxn[1] = {j};
01143 PetscScalar v[1];
01144 MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
01145 return v[0];
01146 }
01147
01148
01150
01155 template <class T, class Prop, class Allocator>
01156 inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
01157 ::Copy(const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
01158 {
01159 this->mpi_communicator_ = A.mpi_communicator_;
01160 PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
01161 this->petsc_matrix_deallocated_ = false;
01162 }
01163
01164
01166
01171 template <class T, class Prop, class Allocator>
01172 inline Matrix<T, Prop, PETScMPIAIJ, Allocator>&
01173 Matrix<T, Prop, PETScMPIAIJ, Allocator>
01174 ::operator= (const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
01175 {
01176 this->Copy(A);
01177 return *this;
01178 }
01179
01180
01182
01187 template <class T, class Prop, class Allocator>
01188 template <class T0, class Allocator0>
01189 inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
01190 ::Copy(const Matrix<T0, General, RowSparse, Allocator0>& A)
01191 {
01192 this->Clear();
01193
01194 int ierr;
01195 ierr = MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
01196 CHKERRABORT(this->mpi_communicator_, ierr);
01197
01198 int ma = A.GetM();
01199 int na = A.GetN();
01200 int nnz = A.GetDataSize();
01201 double *value = A.GetData();
01202 int *column = A.GetInd();
01203 int *ptr = A.GetPtr();
01204
01205 this->m_ = ma;
01206 this->n_ = na;
01207
01208 ierr = MatSetType(this->petsc_matrix_, MATMPIAIJ);
01209 CHKERRABORT(this->mpi_communicator_, ierr);
01210 ierr = MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE,
01211 this->m_, this->n_);
01212 CHKERRABORT(this->mpi_communicator_, ierr);
01213 int loc_start, loc_end;
01214 ierr = MatGetOwnershipRange(this->petsc_matrix_,
01215 &loc_start, &loc_end);
01216 CHKERRABORT(this->mpi_communicator_, ierr);
01217 for (int i = loc_start; i < loc_end; i++)
01218 for (int j = ptr[i]; j < ptr[i + 1]; j++)
01219 ierr = MatSetValues(this->petsc_matrix_, 1, &i, 1, &column[j],
01220 &value[j], INSERT_VALUES);
01221 ierr = MatAssemblyBegin(this->petsc_matrix_,MAT_FINAL_ASSEMBLY);
01222 CHKERRABORT(this->mpi_communicator_, ierr);
01223 ierr = MatAssemblyEnd(this->petsc_matrix_,MAT_FINAL_ASSEMBLY);
01224 CHKERRABORT(this->mpi_communicator_, ierr);
01225 }
01226
01227
01229
01232 template <class T, class Prop, class Allocator>
01233 int Matrix<T, Prop, PETScMPIAIJ, Allocator>::GetLocalM() const
01234 {
01235 int ierr;
01236 int m;
01237 ierr = MatGetLocalSize(this->petsc_matrix_, &m, PETSC_NULL);
01238 CHKERRABORT(this->mpi_communicator_, ierr);
01239 return m;
01240 }
01241
01242
01244
01247 template <class T, class Prop, class Allocator>
01248 int Matrix<T, Prop, PETScMPIAIJ, Allocator>::GetLocalN() const
01249 {
01250 int ierr;
01251 int n;
01252 ierr = MatGetLocalSize(this->petsc_matrix_, PETSC_NULL, &n);
01253 CHKERRABORT(this->mpi_communicator_, ierr);
01254 return n;
01255 }
01256
01257
01259
01264 template <class T, class Prop, class Allocator>
01265 void Matrix<T, Prop, PETScMPIAIJ, Allocator>::Print() const
01266 {
01267 int ierr;
01268 ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_WORLD);
01269 CHKERRABORT(this->mpi_communicator_, ierr);
01270 }
01271
01272
01273 }
01274
01275
01276 #define SELDON_FILE_MATRIX_PETSCMATRIX_CXX
01277 #endif