Ifpack Package Browser (Single Doxygen Collection)  Development
Ifpack_SparseContainer.h
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK_SPARSECONTAINER_H
44 #define IFPACK_SPARSECONTAINER_H
45 
46 #include "Ifpack_Container.h"
47 #include "Epetra_IntSerialDenseVector.h"
48 #include "Epetra_MultiVector.h"
49 #include "Epetra_Vector.h"
50 #include "Epetra_Map.h"
51 #include "Epetra_RowMatrix.h"
52 #include "Epetra_CrsMatrix.h"
53 #include "Epetra_LinearProblem.h"
54 #include "Epetra_IntSerialDenseVector.h"
55 #include "Teuchos_ParameterList.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #ifdef HAVE_MPI
58 #include "Epetra_MpiComm.h"
59 #else
60 #include "Epetra_SerialComm.h"
61 #endif
62 
92 template<typename T>
94 
95 public:
96 
98  Ifpack_SparseContainer(const int NumRows, const int NumVectors = 1);
100 
103 
105  virtual ~Ifpack_SparseContainer();
107 
109 
113 
115  virtual int NumRows() const;
117 
119  virtual int NumVectors() const
120  {
121  return(NumVectors_);
122  }
123 
125  virtual int SetNumVectors(const int NumVectors_in)
126  {
127  if (NumVectors_ != NumVectors_in)
128  {
129  NumVectors_=NumVectors_in;
130  LHS_=Teuchos::rcp(new Epetra_MultiVector(*Map_,NumVectors_));
131  RHS_=Teuchos::rcp(new Epetra_MultiVector(*Map_,NumVectors_));
132  }
133  return(0);
134  }
135 
137  virtual double& LHS(const int i, const int Vector = 0);
138 
140  virtual double& RHS(const int i, const int Vector = 0);
141 
143 
152  virtual int& ID(const int i);
153 
155  virtual int SetMatrixElement(const int row, const int col,
156  const double value);
157 
158 
160  virtual bool IsInitialized() const
161  {
162  return(IsInitialized_);
163  }
164 
166  virtual bool IsComputed() const
167  {
168  return(IsComputed_);
169  }
170 
172  virtual int SetParameters(Teuchos::ParameterList& List);
173 
175  virtual const char* Label() const
176  {
177  return(Label_.c_str());
178  }
179 
181  Teuchos::RCP<const Epetra_Map> Map() const
182  {
183  return(Map_);
184  }
185 
187  Teuchos::RCP<const Epetra_MultiVector> LHS() const
188  {
189  return(LHS_);
190  }
191 
193  Teuchos::RCP<const Epetra_MultiVector> RHS() const
194  {
195  return(RHS_);
196  }
197 
199  Teuchos::RCP<const Epetra_CrsMatrix> Matrix() const
200  {
201  return(Matrix_);
202  }
203 
206  {
207  return(GID_);
208  }
209 
211  Teuchos::RCP<const T> Inverse() const
212  {
213  return(Inverse_);
214  }
216 
218 
225  virtual int Initialize();
227  virtual int Compute(const Epetra_RowMatrix& Matrix_in);
229  virtual int Apply();
230 
232  virtual int ApplyInverse();
233 
235 
237  virtual int Destroy();
240 
242  virtual double InitializeFlops() const
243  {
244  if (Inverse_ == Teuchos::null)
245  return (0.0);
246  else
247  return(Inverse_->InitializeFlops());
248  }
249 
251  virtual double ComputeFlops() const
252  {
253  if (Inverse_ == Teuchos::null)
254  return (0.0);
255  else
256  return(Inverse_->ComputeFlops());
257  }
258 
260  virtual double ApplyFlops() const
261  {
262  return(ApplyFlops_);
263  }
264 
266  virtual double ApplyInverseFlops() const
267  {
268  if (Inverse_ == Teuchos::null)
269  return (0.0);
270  else
271  return(Inverse_->ApplyInverseFlops());
272  }
273 
275  virtual std::ostream& Print(std::ostream& os) const;
276 
277 private:
278 
280  virtual int Extract(const Epetra_RowMatrix& Matrix_in);
281 
283  int NumRows_;
287  Teuchos::RefCountPtr<Epetra_Map> Map_;
289  Teuchos::RefCountPtr<Epetra_CrsMatrix> Matrix_;
291  Teuchos::RefCountPtr<Epetra_MultiVector> LHS_;
293  Teuchos::RefCountPtr<Epetra_MultiVector> RHS_;
301  Teuchos::RefCountPtr<Epetra_Comm> SerialComm_;
303  Teuchos::RefCountPtr<T> Inverse_;
305  std::string Label_;
306  Teuchos::ParameterList List_;
307  double ApplyFlops_;
308 
309 };
310 
311 //==============================================================================
312 template<typename T>
314 Ifpack_SparseContainer(const int NumRows_in, const int NumVectors_in) :
315  NumRows_(NumRows_in),
316  NumVectors_(NumVectors_in),
317  IsInitialized_(false),
318  IsComputed_(false),
319  ApplyFlops_(0.0)
320 {
321 
322 #ifdef HAVE_MPI
323  SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
324 #else
325  SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
326 #endif
327 
328 }
329 
330 //==============================================================================
331 template<typename T>
334  NumRows_(rhs.NumRows()),
335  NumVectors_(rhs.NumVectors()),
336  IsInitialized_(rhs.IsInitialized()),
337  IsComputed_(rhs.IsComputed())
338 {
339 
340 #ifdef HAVE_MPI
341  SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
342 #else
343  SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
344 #endif
345 
346  if (!rhs.Map().is_null())
347  Map_ = Teuchos::rcp( new Epetra_Map(*rhs.Map()) );
348 
349  if (!rhs.Matrix().is_null())
350  Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(*rhs.Matrix()) );
351 
352  if (!rhs.LHS().is_null())
353  LHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.LHS()) );
354 
355  if (!rhs.RHS().is_null())
356  RHS_ = Teuchos::rcp( new Epetra_MultiVector(*rhs.RHS()) );
357 
358 }
359 //==============================================================================
360 template<typename T>
362 {
363  Destroy();
364 }
365 
366 //==============================================================================
367 template<typename T>
369 {
370  if (IsInitialized() == false)
371  return(0);
372  else
373  return(NumRows_);
374 }
375 
376 //==============================================================================
377 template<typename T>
379 {
380 
381  if (IsInitialized_ == true)
382  Destroy();
383 
384  IsInitialized_ = false;
385 
386 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
387  Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
388 #endif
389 
390  LHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
391  RHS_ = Teuchos::rcp( new Epetra_MultiVector(*Map_,NumVectors_) );
392  GID_.Reshape(NumRows_,1);
393 
394  Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy,*Map_,0) );
395 
396  // create the inverse
397  Inverse_ = Teuchos::rcp( new T(Matrix_.get()) );
398 
399  if (Inverse_ == Teuchos::null)
400  IFPACK_CHK_ERR(-5);
401 
402  IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
403 
404  // Call Inverse_->Initialize() in Compute(). This saves
405  // some time, because I can extract the diagonal blocks faster,
406  // and only once.
407 
408  Label_ = "Ifpack_SparseContainer";
409 
410  IsInitialized_ = true;
411  return(0);
412 
413 }
414 
415 //==============================================================================
416 template<typename T>
417 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
418 {
419  return(((*LHS_)(Vector))->Values()[i]);
420 }
421 
422 //==============================================================================
423 template<typename T>
424 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
425 {
426  return(((*RHS_)(Vector))->Values()[i]);
427 }
428 
429 //==============================================================================
430 template<typename T>
432 SetMatrixElement(const int row, const int col, const double value)
433 {
434  if (!IsInitialized())
435  IFPACK_CHK_ERR(-3); // problem not shaped yet
436 
437  if ((row < 0) || (row >= NumRows())) {
438  IFPACK_CHK_ERR(-2); // not in range
439  }
440 
441  if ((col < 0) || (col >= NumRows())) {
442  IFPACK_CHK_ERR(-2); // not in range
443  }
444 
445 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
446  if(Matrix_->RowMatrixRowMap().GlobalIndicesInt()) {
447  int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
448  if (ierr < 0) {
449  ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
450  if (ierr < 0)
451  IFPACK_CHK_ERR(-1);
452  }
453  }
454  else
455 #endif
456 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
457  if(Matrix_->RowMatrixRowMap().GlobalIndicesLongLong()) {
458  long long col_LL = col;
459  int ierr = Matrix_->InsertGlobalValues(row,1,(double*)&value,&col_LL);
460  if (ierr < 0) {
461  ierr = Matrix_->SumIntoGlobalValues(row,1,(double*)&value,&col_LL);
462  if (ierr < 0)
463  IFPACK_CHK_ERR(-1);
464  }
465  }
466  else
467 #endif
468  throw "Ifpack_SparseContainer<T>::SetMatrixElement: GlobalIndices type unknown";
469 
470  return(0);
471 
472 }
473 
474 //==============================================================================
475 template<typename T>
477 {
478 
479  IsComputed_ = false;
480  if (!IsInitialized()) {
481  IFPACK_CHK_ERR(Initialize());
482  }
483 
484  // extract the submatrices
485  IFPACK_CHK_ERR(Extract(Matrix_in));
486 
487  // initialize the inverse operator
488  IFPACK_CHK_ERR(Inverse_->Initialize());
489 
490  // compute the inverse operator
491  IFPACK_CHK_ERR(Inverse_->Compute());
492 
493  Label_ = "Ifpack_SparseContainer";
494 
495  IsComputed_ = true;
496 
497  return(0);
498 }
499 
500 //==============================================================================
501 template<typename T>
503 {
504  if (IsComputed() == false) {
505  IFPACK_CHK_ERR(-3); // not yet computed
506  }
507 
508  IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
509 
510  ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros64();
511  return(0);
512 }
513 
514 //==============================================================================
515 template<typename T>
517 {
518  if (!IsComputed())
519  IFPACK_CHK_ERR(-1);
520 
521  IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
522 
523  return(0);
524 }
525 
526 
527 //==============================================================================
528 template<typename T>
530 {
531  IsInitialized_ = false;
532  IsComputed_ = false;
533  return(0);
534 }
535 
536 //==============================================================================
537 template<typename T>
539 {
540  return(GID_[i]);
541 }
542 
543 //==============================================================================
544 template<typename T>
546 SetParameters(Teuchos::ParameterList& List)
547 {
548  List_ = List;
549  return(0);
550 }
551 
552 //==============================================================================
553 // FIXME: optimize performances of this guy...
554 template<typename T>
556 {
557 
558  for (int j = 0 ; j < NumRows_ ; ++j) {
559  // be sure that the user has set all the ID's
560  if (ID(j) == -1)
561  IFPACK_CHK_ERR(-1);
562  // be sure that all are local indices
563  if (ID(j) > Matrix_in.NumMyRows())
564  IFPACK_CHK_ERR(-1);
565  }
566 
567  int Length = Matrix_in.MaxNumEntries();
568  std::vector<double> Values;
569  Values.resize(Length);
570  std::vector<int> Indices;
571  Indices.resize(Length);
572 
573  for (int j = 0 ; j < NumRows_ ; ++j) {
574 
575  int LRID = ID(j);
576 
577  int NumEntries;
578 
579  int ierr =
580  Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries,
581  &Values[0], &Indices[0]);
582  IFPACK_CHK_ERR(ierr);
583 
584  for (int k = 0 ; k < NumEntries ; ++k) {
585 
586  int LCID = Indices[k];
587 
588  // skip off-processor elements
589  if (LCID >= Matrix_in.NumMyRows())
590  continue;
591 
592  // for local column IDs, look for each ID in the list
593  // of columns hosted by this object
594  // FIXME: use STL
595  int jj = -1;
596  for (int kk = 0 ; kk < NumRows_ ; ++kk)
597  if (ID(kk) == LCID)
598  jj = kk;
599 
600  if (jj != -1)
601  SetMatrixElement(j,jj,Values[k]);
602 
603  }
604  }
605 
606  IFPACK_CHK_ERR(Matrix_->FillComplete());
607 
608  return(0);
609 }
610 
611 //==============================================================================
612 template<typename T>
613 std::ostream& Ifpack_SparseContainer<T>::Print(std::ostream & os) const
614 {
615  using std::endl;
616 
617  os << "================================================================================" << endl;
618  os << "Ifpack_SparseContainer" << endl;
619  os << "Number of rows = " << NumRows() << endl;
620  os << "Number of vectors = " << NumVectors() << endl;
621  os << "IsInitialized() = " << IsInitialized() << endl;
622  os << "IsComputed() = " << IsComputed() << endl;
623  os << "Flops in Initialize() = " << InitializeFlops() << endl;
624  os << "Flops in Compute() = " << ComputeFlops() << endl;
625  os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl;
626  os << "================================================================================" << endl;
627  os << endl;
628 
629  return(os);
630 }
631 #endif // IFPACK_SPARSECONTAINER_H
int NumRows_
Number of rows in the local matrix.
Ifpack_SparseContainer & operator=(const Ifpack_SparseContainer< T > &rhs)
Operator =.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all necessary parameters.
virtual double InitializeFlops() const
Returns the flops in Compute().
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
virtual double ApplyFlops() const
Returns the flops in Apply().
Epetra_IntSerialDenseVector GID_
Contains the subrows/subcols of A that will be inserted in Matrix_.
virtual double & LHS(const int i, const int Vector=0)
Returns the i-th component of the vector Vector of LHS.
virtual int Initialize()
Initializes the container, by completing all the operations based on matrix structure.
Teuchos::RefCountPtr< T > Inverse_
Pointer to an Ifpack_Preconditioner object whose ApplyInverse() defined the action of the inverse of ...
bool IsInitialized_
If true, the container has been successfully initialized.
const int NumVectors
Definition: performance.cpp:71
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
Teuchos::RefCountPtr< Epetra_Map > Map_
Linear map on which the local matrix is based.
Teuchos::RCP< const Epetra_CrsMatrix > Matrix() const
Returns a pointer to the internally stored matrix.
Teuchos::RefCountPtr< Epetra_CrsMatrix > Matrix_
Pointer to the local matrix.
virtual double & RHS(const int i, const int Vector=0)
Returns the i-th component of the vector Vector of RHS.
Teuchos::RCP< const Epetra_Map > Map() const
Returns a pointer to the internally stored map.
Teuchos::RefCountPtr< Epetra_Comm > SerialComm_
Serial communicator (containing only MPI_COMM_SELF if MPI is used).
Teuchos::ParameterList List_
bool IsComputed_
If true, the container has been successfully computed.
virtual int MaxNumEntries() const=0
std::string Label_
Label for this object.
Ifpack_SparseContainer: a class for storing and solving linear systems using sparse matrices...
Teuchos::RCP< const Epetra_MultiVector > RHS() const
Returns a pointer to the internally stored rhs multi-vector.
virtual int SetNumVectors(const int NumVectors_in)
Sets the number of vectors for LHS/RHS.
virtual int ApplyInverse()
Apply the inverse of the matrix to RHS, result is stored in LHS.
virtual const char * Label() const
Returns the label of this container.
Ifpack_SparseContainer(const int NumRows, const int NumVectors=1)
Constructor.
virtual int NumMyRows() const=0
const Epetra_IntSerialDenseVector & ID() const
Returns a pointer to the internally stored ID&#39;s.
#define false
Teuchos::RCP< const Epetra_MultiVector > LHS() const
Returns a pointer to the internally stored solution multi-vector.
virtual ~Ifpack_SparseContainer()
Destructor.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
int NumVectors_
Number of vectors in the local linear system.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
Teuchos::RefCountPtr< Epetra_MultiVector > RHS_
right-hand side for local problems.
virtual int Destroy()
Destroys all data.
Teuchos::RefCountPtr< Epetra_MultiVector > LHS_
Solution vector.
Ifpack_Container: a pure virtual class for creating and solving local linear problems.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Teuchos::RCP< const T > Inverse() const
Returns a pointer to the internally stored inverse operator.
virtual double ComputeFlops() const
Returns the flops in Compute().
#define IFPACK_CHK_ERR(ifpack_err)
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
virtual int Extract(const Epetra_RowMatrix &Matrix_in)
Extract the submatrices identified by the ID set int ID().
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
virtual int Apply()
Apply the matrix to RHS, result is stored in LHS.