Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
92template<typename T>
94
95public:
96
98
99 Ifpack_SparseContainer(const int NumRows, const int NumVectors = 1);
100
103
105 virtual ~Ifpack_SparseContainer();
107
109
113
115
116 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
238 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
277private:
278
280 virtual int Extract(const Epetra_RowMatrix& Matrix_in);
281
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_;
308
309};
310
311//==============================================================================
312template<typename T>
314Ifpack_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//==============================================================================
331template<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//==============================================================================
360template<typename T>
362{
363 Destroy();
364}
365
366//==============================================================================
367template<typename T>
369{
370 if (IsInitialized() == false)
371 return(0);
372 else
373 return(NumRows_);
374}
375
376//==============================================================================
377template<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//==============================================================================
420template<typename T>
421double& Ifpack_SparseContainer<T>::LHS(const int i, const int Vector)
422{
423 return(((*LHS_)(Vector))->Values()[i]);
424}
425
426//==============================================================================
427template<typename T>
428double& Ifpack_SparseContainer<T>::RHS(const int i, const int Vector)
429{
430 return(((*RHS_)(Vector))->Values()[i]);
431}
432
433//==============================================================================
434template<typename T>
436SetMatrixElement(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//==============================================================================
479template<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//==============================================================================
505template<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//==============================================================================
519template<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//==============================================================================
532template<typename T>
534{
535 IsInitialized_ = false;
536 IsComputed_ = false;
537 return(0);
538}
539
540//==============================================================================
541template<typename T>
543{
544 return(GID_[i]);
545}
546
547//==============================================================================
548template<typename T>
550SetParameters(Teuchos::ParameterList& List)
551{
552 List_ = List;
553 return(0);
554}
555
556//==============================================================================
557// FIXME: optimize performances of this guy...
558template<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//==============================================================================
616template<typename T>
617std::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
#define IFPACK_CHK_ERR(ifpack_err)
virtual int NumMyRows() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Ifpack_Container: a pure virtual class for creating and solving local linear problems.
Ifpack_SparseContainer: a class for storing and solving linear systems using sparse matrices.
Teuchos::RefCountPtr< Epetra_CrsMatrix > Matrix_
Pointer to the local matrix.
Teuchos::RefCountPtr< T > Inverse_
Pointer to an Ifpack_Preconditioner object whose ApplyInverse() defined the action of the inverse of ...
virtual double InitializeFlops() const
Returns the flops in Compute().
std::string Label_
Label for this object.
virtual double ApplyFlops() const
Returns the flops in Apply().
virtual int Apply()
Apply the matrix to RHS, result is stored in LHS.
Teuchos::RefCountPtr< Epetra_MultiVector > RHS_
right-hand side for local problems.
virtual double ComputeFlops() const
Returns the flops in Compute().
Teuchos::RefCountPtr< Epetra_Comm > SerialComm_
Serial communicator (containing only MPI_COMM_SELF if MPI is used).
virtual int SetNumVectors(const int NumVectors_in)
Sets the number of vectors for LHS/RHS.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
int NumVectors_
Number of vectors in the local linear system.
virtual int SetMatrixElement(const int row, const int col, const double value)
Set the matrix element (row,col) to value.
virtual bool IsComputed() const
Returns true is the container has been successfully computed.
virtual int NumRows() const
Returns the number of rows of the matrix and LHS/RHS.
const Epetra_IntSerialDenseVector & ID() const
Returns a pointer to the internally stored ID's.
virtual int Initialize()
Initializes the container, by completing all the operations based on matrix structure.
Teuchos::ParameterList List_
virtual double & LHS(const int i, const int Vector=0)
Returns the i-th component of the vector Vector of LHS.
Teuchos::RCP< const Epetra_MultiVector > LHS() const
Returns a pointer to the internally stored solution multi-vector.
virtual int Destroy()
Destroys all data.
virtual bool IsInitialized() const
Returns true is the container has been successfully initialized.
bool IsComputed_
If true, the container has been successfully computed.
virtual ~Ifpack_SparseContainer()
Destructor.
virtual double ApplyInverseFlops() const
Returns the flops in ApplyInverse().
virtual int Compute(const Epetra_RowMatrix &Matrix_in)
Finalizes the linear system matrix and prepares for the application of the inverse.
int NumRows_
Number of rows in the local matrix.
Ifpack_SparseContainer & operator=(const Ifpack_SparseContainer< T > &rhs)
Operator =.
virtual int NumVectors() const
Returns the number of vectors in LHS/RHS.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all necessary parameters.
Teuchos::RCP< const Epetra_Map > Map() const
Returns a pointer to the internally stored map.
virtual int ApplyInverse()
Apply the inverse of the matrix to RHS, result is stored in LHS.
Teuchos::RCP< const T > Inverse() const
Returns a pointer to the internally stored inverse operator.
Teuchos::RCP< const Epetra_MultiVector > RHS() const
Returns a pointer to the internally stored rhs multi-vector.
virtual int Extract(const Epetra_RowMatrix &Matrix_in)
Extract the submatrices identified by the ID set int ID().
Ifpack_SparseContainer(const int NumRows, const int NumVectors=1)
Constructor.
Teuchos::RefCountPtr< Epetra_Map > Map_
Linear map on which the local matrix is based.
Teuchos::RefCountPtr< Epetra_MultiVector > LHS_
Solution vector.
Epetra_IntSerialDenseVector GID_
Contains the subrows/subcols of A that will be inserted in Matrix_.
Teuchos::RCP< const Epetra_CrsMatrix > Matrix() const
Returns a pointer to the internally stored matrix.
virtual double & RHS(const int i, const int Vector=0)
Returns the i-th component of the vector Vector of RHS.
virtual int & ID(const int i)
Returns the ID associated to local row i.
virtual const char * Label() const
Returns the label of this container.
bool IsInitialized_
If true, the container has been successfully initialized.
#define false
const int NumVectors
Definition: performance.cpp:71