IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_SparseContainer.h
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_;
285  int NumVectors_;
287  Teuchos::RefCountPtr<Epetra_Map> Map_;
289  Teuchos::RefCountPtr<Epetra_CrsMatrix> Matrix_;
291  Teuchos::RefCountPtr<Epetra_MultiVector> LHS_;
293  Teuchos::RefCountPtr<Epetra_MultiVector> RHS_;
297  bool IsInitialized_;
299  bool IsComputed_;
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 #if defined(HAVE_TEUCHOSCORE_CXX11)
395  Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy,*Map_,0) );
396 #else
397  Matrix_ = Teuchos::rcp( new Epetra_CrsMatrix(::Copy,*Map_,0) );
398 #endif
399 
400  // create the inverse
401  Inverse_ = Teuchos::rcp( new T(Matrix_.get()) );
402 
403  if (Inverse_ == Teuchos::null)
404  IFPACK_CHK_ERR(-5);
405 
406  IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
407 
408  // Call Inverse_->Initialize() in Compute(). This saves
409  // some time, because I can extract the diagonal blocks faster,
410  // and only once.
411 
412  Label_ = "Ifpack_SparseContainer";
413 
414  IsInitialized_ = true;
415  return(0);
416 
417 }
418 
419 //==============================================================================
420 template<typename T>
421 double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
422 {
423  return(((*LHS_)(Vector))->Values()[i]);
424 }
425 
426 //==============================================================================
427 template<typename T>
428 double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
429 {
430  return(((*RHS_)(Vector))->Values()[i]);
431 }
432 
433 //==============================================================================
434 template<typename T>
436 SetMatrixElement(const int row, const int col, const double value)
437 {
438  if (!IsInitialized())
439  IFPACK_CHK_ERR(-3); // problem not shaped yet
440 
441  if ((row < 0) || (row >= NumRows())) {
442  IFPACK_CHK_ERR(-2); // not in range
443  }
444 
445  if ((col < 0) || (col >= NumRows())) {
446  IFPACK_CHK_ERR(-2); // not in range
447  }
448 
449 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
450  if(Matrix_->RowMatrixRowMap().GlobalIndicesInt()) {
451  int ierr = Matrix_->InsertGlobalValues((int)row,1,(double*)&value,(int*)&col);
452  if (ierr < 0) {
453  ierr = Matrix_->SumIntoGlobalValues((int)row,1,(double*)&value,(int*)&col);
454  if (ierr < 0)
455  IFPACK_CHK_ERR(-1);
456  }
457  }
458  else
459 #endif
460 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
461  if(Matrix_->RowMatrixRowMap().GlobalIndicesLongLong()) {
462  long long col_LL = col;
463  int ierr = Matrix_->InsertGlobalValues(row,1,(double*)&value,&col_LL);
464  if (ierr < 0) {
465  ierr = Matrix_->SumIntoGlobalValues(row,1,(double*)&value,&col_LL);
466  if (ierr < 0)
467  IFPACK_CHK_ERR(-1);
468  }
469  }
470  else
471 #endif
472  throw "Ifpack_SparseContainer<T>::SetMatrixElement: GlobalIndices type unknown";
473 
474  return(0);
475 
476 }
477 
478 //==============================================================================
479 template<typename T>
481 {
482 
483  IsComputed_ = false;
484  if (!IsInitialized()) {
485  IFPACK_CHK_ERR(Initialize());
486  }
487 
488  // extract the submatrices
489  IFPACK_CHK_ERR(Extract(Matrix_in));
490 
491  // initialize the inverse operator
492  IFPACK_CHK_ERR(Inverse_->Initialize());
493 
494  // compute the inverse operator
495  IFPACK_CHK_ERR(Inverse_->Compute());
496 
497  Label_ = "Ifpack_SparseContainer";
498 
499  IsComputed_ = true;
500 
501  return(0);
502 }
503 
504 //==============================================================================
505 template<typename T>
507 {
508  if (IsComputed() == false) {
509  IFPACK_CHK_ERR(-3); // not yet computed
510  }
511 
512  IFPACK_CHK_ERR(Matrix_->Apply(*RHS_, *LHS_));
513 
514  ApplyFlops_ += 2 * Matrix_->NumGlobalNonzeros64();
515  return(0);
516 }
517 
518 //==============================================================================
519 template<typename T>
521 {
522  if (!IsComputed())
523  IFPACK_CHK_ERR(-1);
524 
525  IFPACK_CHK_ERR(Inverse_->ApplyInverse(*RHS_, *LHS_));
526 
527  return(0);
528 }
529 
530 
531 //==============================================================================
532 template<typename T>
534 {
535  IsInitialized_ = false;
536  IsComputed_ = false;
537  return(0);
538 }
539 
540 //==============================================================================
541 template<typename T>
543 {
544  return(GID_[i]);
545 }
546 
547 //==============================================================================
548 template<typename T>
550 SetParameters(Teuchos::ParameterList& List)
551 {
552  List_ = List;
553  return(0);
554 }
555 
556 //==============================================================================
557 // FIXME: optimize performances of this guy...
558 template<typename T>
560 {
561 
562  for (int j = 0 ; j < NumRows_ ; ++j) {
563  // be sure that the user has set all the ID's
564  if (ID(j) == -1)
565  IFPACK_CHK_ERR(-1);
566  // be sure that all are local indices
567  if (ID(j) > Matrix_in.NumMyRows())
568  IFPACK_CHK_ERR(-1);
569  }
570 
571  int Length = Matrix_in.MaxNumEntries();
572  std::vector<double> Values;
573  Values.resize(Length);
574  std::vector<int> Indices;
575  Indices.resize(Length);
576 
577  for (int j = 0 ; j < NumRows_ ; ++j) {
578 
579  int LRID = ID(j);
580 
581  int NumEntries;
582 
583  int ierr =
584  Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries,
585  &Values[0], &Indices[0]);
586  IFPACK_CHK_ERR(ierr);
587 
588  for (int k = 0 ; k < NumEntries ; ++k) {
589 
590  int LCID = Indices[k];
591 
592  // skip off-processor elements
593  if (LCID >= Matrix_in.NumMyRows())
594  continue;
595 
596  // for local column IDs, look for each ID in the list
597  // of columns hosted by this object
598  // FIXME: use STL
599  int jj = -1;
600  for (int kk = 0 ; kk < NumRows_ ; ++kk)
601  if (ID(kk) == LCID)
602  jj = kk;
603 
604  if (jj != -1)
605  SetMatrixElement(j,jj,Values[k]);
606 
607  }
608  }
609 
610  IFPACK_CHK_ERR(Matrix_->FillComplete());
611 
612  return(0);
613 }
614 
615 //==============================================================================
616 template<typename T>
617 std::ostream& Ifpack_SparseContainer<T>::Print(std::ostream & os) const
618 {
619  using std::endl;
620 
621  os << "================================================================================" << endl;
622  os << "Ifpack_SparseContainer" << endl;
623  os << "Number of rows = " << NumRows() << endl;
624  os << "Number of vectors = " << NumVectors() << endl;
625  os << "IsInitialized() = " << IsInitialized() << endl;
626  os << "IsComputed() = " << IsComputed() << endl;
627  os << "Flops in Initialize() = " << InitializeFlops() << endl;
628  os << "Flops in Compute() = " << ComputeFlops() << endl;
629  os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl;
630  os << "================================================================================" << endl;
631  os << endl;
632 
633  return(os);
634 }
635 #endif // IFPACK_SPARSECONTAINER_H
virtual double ComputeFlops() const
Returns the flops in Compute().
virtual double InitializeFlops() const
Returns the flops in Compute().
Teuchos::RCP< const Epetra_MultiVector > LHS() const
Returns a pointer to the internally stored solution multi-vector.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
Ifpack_SparseContainer & operator=(const Ifpack_SparseContainer< T > &rhs)
Operator =.
Teuchos::RCP< const Epetra_MultiVector > RHS() const
Returns a pointer to the internally stored rhs multi-vector.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all necessary parameters.
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
Teuchos::RCP< const Epetra_CrsMatrix > Matrix() const
Returns a pointer to the internally stored 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.
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
virtual const char * Label() const
Returns the label of this container.
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.
virtual int MaxNumEntries() const =0
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
Ifpack_SparseContainer: a class for storing and solving linear systems using sparse matrices...
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.
Ifpack_SparseContainer(const int NumRows, const int NumVectors=1)
Constructor.
virtual int NumMyRows() const =0
virtual double ApplyFlops() const
Returns the flops in Apply().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.
virtual ~Ifpack_SparseContainer()
Destructor.
Teuchos::RCP< const T > Inverse() const
Returns a pointer to the internally stored inverse operator.
const Epetra_IntSerialDenseVector & ID() const
Returns a pointer to the internally stored ID&#39;s.
virtual int Destroy()
Destroys all data.
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
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
virtual int Apply()
Apply the matrix to RHS, result is stored in LHS.