00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef IFPACK_SPARSECONTAINER_H
00031 #define IFPACK_SPARSECONTAINER_H
00032
00033 #include "Ifpack_Container.h"
00034 #include "Epetra_IntSerialDenseVector.h"
00035 #include "Epetra_MultiVector.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_RowMatrix.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Epetra_LinearProblem.h"
00041 #include "Epetra_IntSerialDenseVector.h"
00042 #include "Teuchos_ParameterList.hpp"
00043 #include "Teuchos_RefCountPtr.hpp"
00044 #ifdef HAVE_MPI
00045 #include "Epetra_MpiComm.h"
00046 #else
00047 #include "Epetra_SerialComm.h"
00048 #endif
00049
00079 template<typename T>
00080 class Ifpack_SparseContainer : public Ifpack_Container {
00081
00082 public:
00083
00085
00086 Ifpack_SparseContainer(const int NumRows, const int NumVectors = 1);
00087
00089 Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs);
00090
00092 virtual ~Ifpack_SparseContainer();
00094
00096
00098 Ifpack_SparseContainer& operator=(const Ifpack_SparseContainer<T>& rhs);
00100
00102
00103 virtual int NumRows() const;
00104
00106 virtual int NumVectors() const
00107 {
00108 return(NumVectors_);
00109 }
00110
00112 virtual int SetNumVectors(const int NumVectors_in)
00113 {
00114 if (NumVectors_ == NumVectors_in)
00115 return(0);
00116 IFPACK_CHK_ERR(-99);
00117 }
00118
00120 virtual double& LHS(const int i, const int Vector = 0);
00121
00123 virtual double& RHS(const int i, const int Vector = 0);
00124
00126
00135 virtual int& ID(const int i);
00136
00138 virtual int SetMatrixElement(const int row, const int col,
00139 const double value);
00140
00141
00143 virtual bool IsInitialized() const
00144 {
00145 return(IsInitialized_);
00146 }
00147
00149 virtual bool IsComputed() const
00150 {
00151 return(IsComputed_);
00152 }
00153
00155 virtual int SetParameters(Teuchos::ParameterList& List);
00156
00158 virtual const char* Label() const
00159 {
00160 return(Label_.c_str());
00161 }
00162
00164 const Epetra_Map* Map() const
00165 {
00166 return(Map_);
00167 }
00168
00170 const Epetra_MultiVector* LHS() const
00171 {
00172 return(LHS_);
00173 }
00174
00176 const Epetra_MultiVector* RHS() const
00177 {
00178 return(RHS_);
00179 }
00180
00182 const Epetra_CrsMatrix* Matrix() const
00183 {
00184 return(Matrix_);
00185 }
00186
00188 const Epetra_IntSerialDenseVector* ID() const
00189 {
00190 return(GID_);
00191 }
00192
00194 const T* Inverse() const
00195 {
00196 return(Inverse_);
00197 }
00199
00201
00208 virtual int Initialize();
00210 virtual int Compute(const Epetra_RowMatrix& Matrix_in);
00212 virtual int Apply();
00213
00215 virtual int ApplyInverse();
00216
00218
00220
00221 virtual int Destroy();
00223
00225 virtual double InitializeFlops() const
00226 {
00227 if (Inverse_ == Teuchos::null)
00228 return (0.0);
00229 else
00230 return(Inverse_->InitializeFlops());
00231 }
00232
00234 virtual double ComputeFlops() const
00235 {
00236 if (Inverse_ == Teuchos::null)
00237 return (0.0);
00238 else
00239 return(Inverse_->ComputeFlops());
00240 }
00241
00243 virtual double ApplyFlops() const
00244 {
00245 return(ApplyFlops_);
00246 }
00247
00249 virtual double ApplyInverseFlops() const
00250 {
00251 if (Inverse_ == Teuchos::null)
00252 return (0.0);
00253 else
00254 return(Inverse_->ApplyInverseFlops());
00255 }
00256
00258 virtual ostream& Print(std::ostream& os) const;
00259
00260 private:
00261
00263 virtual int Extract(const Epetra_RowMatrix& Matrix_in);
00264
00266 int NumRows_;
00268 int NumVectors_;
00270 Teuchos::RefCountPtr<Epetra_Map> Map_;
00272 Teuchos::RefCountPtr<Epetra_CrsMatrix> Matrix_;
00274 Teuchos::RefCountPtr<Epetra_MultiVector> LHS_;
00276 Teuchos::RefCountPtr<Epetra_MultiVector> RHS_;
00278 Epetra_IntSerialDenseVector GID_;
00280 bool IsInitialized_;
00282 bool IsComputed_;
00284 Teuchos::RefCountPtr<Epetra_Comm> SerialComm_;
00286 Teuchos::RefCountPtr<T> Inverse_;
00288 string Label_;
00289 Teuchos::ParameterList List_;
00290 double ApplyFlops_;
00291
00292 };
00293
00294
00295 template<typename T>
00296 Ifpack_SparseContainer<T>::
00297 Ifpack_SparseContainer(const int NumRows_in, const int NumVectors_in) :
00298 NumRows_(NumRows_in),
00299 NumVectors_(NumVectors_in),
00300 IsInitialized_(false),
00301 IsComputed_(false),
00302 ApplyFlops_(0.0)
00303 {
00304
00305 #ifdef HAVE_MPI
00306 SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00307 #else
00308 SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00309 #endif
00310
00311 }
00312
00313
00314 template<typename T>
00315 Ifpack_SparseContainer<T>::
00316 Ifpack_SparseContainer(const Ifpack_SparseContainer<T>& rhs) :
00317 NumRows_(rhs.NumRows()),
00318 NumVectors_(rhs.NumVectors()),
00319 IsInitialized_(rhs.IsInitialized()),
00320 IsComputed_(rhs.IsComputed())
00321 {
00322
00323 #ifdef HAVE_MPI
00324 SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00325 #else
00326 SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00327 #endif
00328
00329 if (rhs.Map())
00330 Map_ = Teuchos::rcp( new Epetra_Map(*rhs.Map()) );
00331
00332 if (rhs.Matrix())
00333 Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(*rhs.Matrix()) );
00334
00335 if (rhs.LHS())
00336 LHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.LHS()) );
00337
00338 if (rhs.RHS())
00339 RHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.RHS()) );
00340
00341 }
00342
00343 template<typename T>
00344 Ifpack_SparseContainer<T>::~Ifpack_SparseContainer()
00345 {
00346 Destroy();
00347 }
00348
00349
00350 template<typename T>
00351 int Ifpack_SparseContainer<T>::NumRows() const
00352 {
00353 if (IsInitialized() == false)
00354 return(0);
00355 else
00356 return(NumRows_);
00357 }
00358
00359
00360 template<typename T>
00361 int Ifpack_SparseContainer<T>::Initialize()
00362 {
00363
00364 if (IsInitialized_ == true)
00365 Destroy();
00366
00367 IsInitialized_ = false;
00368
00369 Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
00370
00371 LHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00372 RHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
00373 GID_.Reshape(NumRows_,1);
00374
00375 Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy,*Map_,0) );
00376
00377
00378 Inverse_ = Teuchos::rcp( new T(Matrix_.get()) );
00379
00380 if (Inverse_ == Teuchos::null)
00381 IFPACK_CHK_ERR(-5);
00382
00383 IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00384
00385
00386
00387
00388
00389 Label_ = "Ifpack_SparseContainer";
00390
00391 IsInitialized_ = true;
00392 return(0);
00393
00394 }
00395
00396
00397 template<typename T>
00398 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
00399 {
00400 return(((*LHS_)(Vector))->Values()[i]);
00401 }
00402
00403
00404 template<typename T>
00405 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
00406 {
00407 return(((*RHS_)(Vector))->Values()[i]);
00408 }
00409
00410
00411 template<typename T>
00412 int Ifpack_SparseContainer<T>::
00413 SetMatrixElement(const int row, const int col, const double value)
00414 {
00415 if (!IsInitialized())
00416 IFPACK_CHK_ERR(-3);
00417
00418 if ((row < 0) || (row >= NumRows())) {
00419 IFPACK_CHK_ERR(-2);
00420 }
00421
00422 if ((col < 0) || (col >= NumRows())) {
00423 IFPACK_CHK_ERR(-2);
00424 }
00425
00426 int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
00427 if (ierr < 0) {
00428 ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
00429 if (ierr < 0)
00430 IFPACK_CHK_ERR(-1);
00431 }
00432
00433 return(0);
00434
00435 }
00436
00437
00438 template<typename T>
00439 int Ifpack_SparseContainer<T>::Compute(const Epetra_RowMatrix& Matrix_in)
00440 {
00441
00442 IsComputed_ = false;
00443 if (!IsInitialized()) {
00444 IFPACK_CHK_ERR(Initialize());
00445 }
00446
00447
00448 IFPACK_CHK_ERR(Extract(Matrix_in));
00449
00450
00451 IFPACK_CHK_ERR(Inverse_->Initialize());
00452
00453
00454 IFPACK_CHK_ERR(Inverse_->Compute());
00455
00456 Label_ = "Ifpack_SparseContainer";
00457
00458 IsComputed_ = true;
00459
00460 return(0);
00461 }
00462
00463
00464 template<typename T>
00465 int Ifpack_SparseContainer<T>::Apply()
00466 {
00467 if (IsComputed() == false) {
00468 IFPACK_CHK_ERR(-3);
00469 }
00470
00471 IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
00472
00473 ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros();
00474 return(0);
00475 }
00476
00477
00478 template<typename T>
00479 int Ifpack_SparseContainer<T>::ApplyInverse()
00480 {
00481 if (!IsComputed())
00482 IFPACK_CHK_ERR(-1);
00483
00484 IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
00485
00486 return(0);
00487 }
00488
00489
00490
00491 template<typename T>
00492 int Ifpack_SparseContainer<T>::Destroy()
00493 {
00494 IsInitialized_ = false;
00495 IsComputed_ = false;
00496 return(0);
00497 }
00498
00499
00500 template<typename T>
00501 int& Ifpack_SparseContainer<T>::ID(const int i)
00502 {
00503 return(GID_[i]);
00504 }
00505
00506
00507 template<typename T>
00508 int Ifpack_SparseContainer<T>::
00509 SetParameters(Teuchos::ParameterList& List)
00510 {
00511 List_ = List;
00512 return(0);
00513 }
00514
00515
00516
00517 template<typename T>
00518 int Ifpack_SparseContainer<T>::Extract(const Epetra_RowMatrix& Matrix_in)
00519 {
00520
00521 for (int j = 0 ; j < NumRows_ ; ++j) {
00522
00523 if (ID(j) == -1)
00524 IFPACK_CHK_ERR(-1);
00525
00526 if (ID(j) > Matrix_in.NumMyRows())
00527 IFPACK_CHK_ERR(-1);
00528 }
00529
00530 int Length = Matrix_in.MaxNumEntries();
00531 std::vector<double> Values;
00532 Values.resize(Length);
00533 std::vector<int> Indices;
00534 Indices.resize(Length);
00535
00536 for (int j = 0 ; j < NumRows_ ; ++j) {
00537
00538 int LRID = ID(j);
00539
00540 int NumEntries;
00541
00542 int ierr =
00543 Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries,
00544 &Values[0], &Indices[0]);
00545 IFPACK_CHK_ERR(ierr);
00546
00547 for (int k = 0 ; k < NumEntries ; ++k) {
00548
00549 int LCID = Indices[k];
00550
00551
00552 if (LCID >= Matrix_in.NumMyRows())
00553 continue;
00554
00555
00556
00557
00558 int jj = -1;
00559 for (int kk = 0 ; kk < NumRows_ ; ++kk)
00560 if (ID(kk) == LCID)
00561 jj = kk;
00562
00563 if (jj != -1)
00564 SetMatrixElement(j,jj,Values[k]);
00565
00566 }
00567 }
00568
00569 IFPACK_CHK_ERR(Matrix_->FillComplete());
00570
00571 return(0);
00572 }
00573
00574
00575 template<typename T>
00576 ostream& Ifpack_SparseContainer<T>::Print(ostream & os) const
00577 {
00578 os << "================================================================================" << endl;
00579 os << "Ifpack_SparseContainer" << endl;
00580 os << "Number of rows = " << NumRows() << endl;
00581 os << "Number of vectors = " << NumVectors() << endl;
00582 os << "IsInitialized() = " << IsInitialized() << endl;
00583 os << "IsComputed() = " << IsComputed() << endl;
00584 os << "Flops in Initialize() = " << InitializeFlops() << endl;
00585 os << "Flops in Compute() = " << ComputeFlops() << endl;
00586 os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl;
00587 os << "================================================================================" << endl;
00588 os << endl;
00589
00590 return(os);
00591 }
00592 #endif // IFPACK_SPARSECONTAINER_H