44#ifndef IFPACK_BLOCKPRECONDITIONER_H
45#define IFPACK_BLOCKPRECONDITIONER_H
59#include "Teuchos_ParameterList.hpp"
60#include "Teuchos_RefCountPtr.hpp"
62#include "Epetra_Map.h"
63#include "Epetra_RowMatrix.h"
64#include "Epetra_MultiVector.h"
65#include "Epetra_Vector.h"
66#include "Epetra_Time.h"
67#include "Epetra_Import.h"
199 virtual const char*
Label()
const;
268 std::ostream&
Print(std::ostream& os)
const;
309#ifdef IFPACK_FLOPCOUNTERS
317 for (
unsigned int i = 0 ; i <
Containers_.size() ; ++i)
327#ifdef IFPACK_FLOPCOUNTERS
332 for (
unsigned int i = 0 ; i <
Containers_.size() ; ++i)
342#ifdef IFPACK_FLOPCOUNTERS
347 for (
unsigned int i = 0 ; i <
Containers_.size() ; ++i) {
428 Teuchos::RefCountPtr< const Epetra_RowMatrix >
Matrix_;
440 Teuchos::RefCountPtr<Ifpack_Graph>
Graph_;
442 Teuchos::RefCountPtr<Epetra_Vector>
W_;
456 IsInitialized_(
false),
461 InitializeTime_(0.0),
463 ApplyInverseTime_(0.0),
464 InitializeFlops_(0.0),
466 ApplyInverseFlops_(0.0),
471 Diagonal_( Matrix_in->Map()),
472 PartitionerType_(
"greedy"),
474 ZeroStartingSolution_(
true),
493 return(Label_.c_str());
501 int ierr = Matrix().Apply(X,Y);
510 return(Matrix().Comm());
518 return(Matrix().OperatorDomainMap());
526 return(Matrix().OperatorRangeMap());
534 if (Partitioner_ == Teuchos::null)
537 NumLocalBlocks_ = Partitioner_->NumLocalParts();
539 Containers_.resize(NumLocalBlocks());
542 Matrix_->ExtractDiagonalCopy(Diagonal_);
544 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
546 int rows = Partitioner_->NumRowsInPart(i);
551 Containers_[i] = Teuchos::rcp(
new T(rows) );
558 for (
int j = 0 ; j < rows ; ++j) {
559 int LRID = (*Partitioner_)(i,j);
560 Containers_[i]->ID(j) = LRID;
568#ifdef SINGLETON_DEBUG
571 for (
int i = 0 ; i < NumLocalBlocks() ; ++i)
572 issing += (
int) ( Partitioner_->NumRowsInPart(i) == 1);
573 printf(
" %d of %d containers are singleton \n",issing,NumLocalBlocks());
584 if (!IsInitialized())
587 Time_.ResetStartTime();
591 if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
598 Importer_ = Teuchos::rcp(
new Epetra_Import(Matrix().RowMatrixColMap(),
599 Matrix().RowMatrixRowMap()) );
604 ComputeTime_ += Time_.ElapsedTime();
622 Time_.ResetStartTime();
626 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
630 Xcopy = Teuchos::rcp( &X,
false );
644 ApplyInverseTime_ += Time_.ElapsedTime();
660 if (ZeroStartingSolution_)
664 if (NumSweeps_ == 1 && ZeroStartingSolution_) {
665 int ierr = DoJacobi(X,Y);
671 for (
int j = 0; j < NumSweeps_ ; j++) {
673 ApplyInverseFlops_ += X.
NumVectors() * 2 * Matrix_->NumGlobalNonzeros64();
675 ApplyInverseFlops_ += X.
NumVectors() * 2 * Matrix_->NumGlobalRows64();
691 if (OverlapLevel_ == 0) {
693 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
695 int rows = Partitioner_->NumRowsInPart(i);
704 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
705 LID = Containers_[i]->ID(j);
707 Containers_[i]->RHS(j,k) = X[k][LID];
718 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
719 LID = Containers_[i]->ID(j);
721 Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
727 int LRID = (*Partitioner_)(i,0);
728 double b = X[0][LRID];
729 double a = Diagonal_[LRID];
730 Y[0][LRID] += DampingFactor_* b/a;
734#ifdef IFPACK_FLOPCOUNTERS
735 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalRows();
742 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
744 int rows = Partitioner_->NumRowsInPart(i);
753 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
754 LID = Containers_[i]->ID(j);
756 Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
765 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
766 LID = Containers_[i]->ID(j);
768 Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
773 int LRID = (*Partitioner_)(i,0);
774 double w = (*W_)[LRID];
775 double b = w * X[0][LRID];
776 double a = Diagonal_[LRID];
778 Y[0][LRID] += DampingFactor_ * w * b / a;
784#ifdef IFPACK_FLOPCOUNTERS
785 ApplyInverseFlops_ +=
NumVectors * 4 * Matrix_->NumGlobalRows();
799 if (ZeroStartingSolution_)
803 for (
int j = 0; j < NumSweeps_ ; j++) {
805 if (j != NumSweeps_ - 1)
821 int Length = Matrix().MaxNumEntries();
822 std::vector<int> Indices(Length);
823 std::vector<double> Values(Length);
825 int NumMyRows = Matrix().NumMyRows();
831 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
835 Y2 = Teuchos::rcp( &Y,
false );
840 Y2->ExtractView(&y2_ptr);
846 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
847 int rows = Partitioner_->NumRowsInPart(i);
850 if (rows!=1 && Containers_[i]->NumRows() == 0)
857 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
858 LID = (*Partitioner_)(i,j);
862 &Values[0], &Indices[0]));
864 for (
int k = 0 ; k < NumEntries ; ++k) {
865 int col = Indices[k];
868 X[kk][LID] -= Values[k] * y2_ptr[kk][col];
876 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
877 LID = Containers_[i]->ID(j);
879 Containers_[i]->RHS(j,k) = X[k][LID];
884#ifdef IFPACK_FLOPCOUNTERS
885 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
888 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
889 LID = Containers_[i]->ID(j);
891 double temp = DampingFactor_ * Containers_[i]->LHS(j,k);
892 y2_ptr[k][LID] += temp;
897 int LRID = (*Partitioner_)(i,0);
898 double b = X[0][LRID];
899 double a = Diagonal_[LRID];
900 y2_ptr[0][LRID]+= DampingFactor_* b/a;
906#ifdef IFPACK_FLOPCOUNTERS
907 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalNonzeros();
908 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalRows();
915 for (
int i = 0 ; i < NumMyRows ; ++i)
916 y_ptr[m][i] = y2_ptr[m][i];
927 if (ZeroStartingSolution_)
931 for (
int j = 0; j < NumSweeps_ ; j++) {
933 if (j != NumSweeps_ - 1)
946 int NumMyRows = Matrix().NumMyRows();
949 int Length = Matrix().MaxNumEntries();
950 std::vector<int> Indices;
951 std::vector<double> Values;
952 Indices.resize(Length);
953 Values.resize(Length);
958 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
962 Y2 = Teuchos::rcp( &Y,
false );
967 Y2->ExtractView(&y2_ptr);
973 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
974 int rows = Partitioner_->NumRowsInPart(i);
976 if (rows !=1 && Containers_[i]->NumRows() == 0)
983 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
984 LID = (*Partitioner_)(i,j);
987 &Values[0], &Indices[0]));
989 for (
int k = 0 ; k < NumEntries ; ++k) {
990 int col = Indices[k];
993 Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
1001 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1002 LID = Containers_[i]->ID(j);
1004 Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1009#ifdef IFPACK_FLOPCOUNTERS
1010 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1013 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1014 LID = Containers_[i]->ID(j);
1016 y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1021 int LRID = (*Partitioner_)(i,0);
1022 double b = Xcopy[0][LRID];
1023 double a = Diagonal_[LRID];
1024 y2_ptr[0][LRID]+= DampingFactor_* b/a;
1028#ifdef IFPACK_FLOPCOUNTERS
1030 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1031 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalRows();
1036 for (
int i = NumLocalBlocks() - 1; i >=0 ; --i) {
1037 int rows = Partitioner_->NumRowsInPart(i);
1038 if (rows != 1 &&Containers_[i]->NumRows() == 0)
1045 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1046 LID = (*Partitioner_)(i,j);
1050 &Values[0], &Indices[0]));
1052 for (
int k = 0 ; k < NumEntries ; ++k) {
1053 int col = Indices[k];
1056 Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
1063 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1064 LID = Containers_[i]->ID(j);
1066 Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1071#ifdef IFPACK_FLOPCOUNTERS
1072 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1075 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1076 LID = Containers_[i]->ID(j);
1078 y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1083 int LRID = (*Partitioner_)(i,0);
1084 double b = Xcopy[0][LRID];
1085 double a = Diagonal_[LRID];
1086 y2_ptr[0][LRID]+= DampingFactor_* b/a;
1090#ifdef IFPACK_FLOPCOUNTERS
1092 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1093 ApplyInverseFlops_ +=
NumVectors * 2 * Matrix_->NumGlobalRows();
1100 for (
int i = 0 ; i < NumMyRows ; ++i)
1101 y_ptr[m][i] = y2_ptr[m][i];
1116 PT =
"Gauss-Seidel";
1118 PT =
"symmetric Gauss-Seidel";
1120 if (!Comm().MyPID()) {
1122 os <<
"================================================================================" << endl;
1123 os <<
"Ifpack_BlockRelaxation, " << PT << endl;
1124 os <<
"Sweeps = " << NumSweeps_ << endl;
1125 os <<
"Damping factor = " << DampingFactor_;
1126 if (ZeroStartingSolution_)
1127 os <<
", using zero starting solution" << endl;
1129 os <<
", using input starting solution" << endl;
1130 os <<
"Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
1132 os <<
"Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1134 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1135 os <<
"----- ------- -------------- ------------ --------" << endl;
1136 os <<
"Initialize() " << std::setw(5) << NumInitialize()
1137 <<
" " << std::setw(15) << InitializeTime()
1138 <<
" " << std::setw(15) << 1.0e-6 * InitializeFlops();
1139 if (InitializeTime() != 0.0)
1140 os <<
" " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
1142 os <<
" " << std::setw(15) << 0.0 << endl;
1143 os <<
"Compute() " << std::setw(5) << NumCompute()
1144 <<
" " << std::setw(15) << ComputeTime()
1145 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops();
1146 if (ComputeTime() != 0.0)
1147 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
1149 os <<
" " << std::setw(15) << 0.0 << endl;
1150 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse()
1151 <<
" " << std::setw(15) << ApplyInverseTime()
1152 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
1153 if (ApplyInverseTime() != 0.0)
1154 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
1156 os <<
" " << std::setw(15) << 0.0 << endl;
1157 os <<
"================================================================================" << endl;
1175 PT =
"Gauss-Seidel";
1177 PT =
"symmetric Gauss-Seidel";
1179 PT = List.get(
"relaxation: type", PT);
1181 if (PT ==
"Jacobi") {
1184 else if (PT ==
"Gauss-Seidel") {
1187 else if (PT ==
"symmetric Gauss-Seidel") {
1190 cerr <<
"Option `relaxation: type' has an incorrect value ("
1191 << PT <<
")'" << endl;
1192 cerr <<
"(file " << __FILE__ <<
", line " << __LINE__ <<
")" << endl;
1196 NumSweeps_ = List.get(
"relaxation: sweeps", NumSweeps_);
1197 DampingFactor_ = List.get(
"relaxation: damping factor",
1199 ZeroStartingSolution_ = List.get(
"relaxation: zero starting solution",
1200 ZeroStartingSolution_);
1201 PartitionerType_ = List.get(
"partitioner: type",
1203 NumLocalBlocks_ = List.get(
"partitioner: local parts",
1206 OverlapLevel_ = List.get(
"partitioner: overlap",
1212 if (NumLocalBlocks_ < 0)
1213 NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
1228 Label_ =
"IFPACK (" + PT2 +
", sweeps="
1240 IsInitialized_ =
false;
1241 Time_.ResetStartTime();
1246 if (PartitionerType_ ==
"linear")
1248 else if (PartitionerType_ ==
"greedy")
1250 else if (PartitionerType_ ==
"metis")
1252 else if (PartitionerType_ ==
"equation")
1254 else if (PartitionerType_ ==
"user")
1256 else if (PartitionerType_ ==
"line")
1268 NumLocalBlocks_ = Partitioner_->NumLocalParts();
1271 W_ = Teuchos::rcp(
new Epetra_Vector(Matrix().RowMatrixRowMap()) );
1274 for (
int i = 0 ; i < NumLocalBlocks() ; ++i) {
1276 for (
int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1277 int LID = (*Partitioner_)(i,j);
1281 W_->Reciprocal(*W_);
1284 if (PartitionerType_ ==
"line") {
1293 Label_ =
"IFPACK (" + PT2 +
", auto-line, sweeps="
1299 InitializeTime_ += Time_.ElapsedTime();
1300 IsInitialized_ =
true;
static const int IFPACK_JACOBI
static const int IFPACK_GS
static const int IFPACK_SGS
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
@ Ifpack_Cheap
cheap estimate
#define IFPACK_RETURN(ifpack_err)
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int NumProc() const=0
int ExtractView(double **A, int *MyLDA) const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
double ** Pointers() const
int PutScalar(double ScalarConstant)
virtual const Epetra_Comm & Comm() const=0
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual int ApplyInverseGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int DoSGS(const Epetra_MultiVector &X, Epetra_MultiVector &Xtmp, Epetra_MultiVector &Y) const
Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation &rhs)
operator= (PRIVATE, should not be used).
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
double InitializeFlops_
Contains the number of flops for Initialize().
Ifpack_BlockRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_BlockRelaxation constructor with given Epetra_RowMatrix.
bool IsComputed_
If true, the preconditioner has been successfully computed.
virtual int NumCompute() const
Returns the number of calls to Compute().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Ifpack_BlockRelaxation(const Ifpack_BlockRelaxation &rhs)
Copy constructor (PRIVATE, should not be used).
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int SetUseTranspose(bool UseTranspose_in)
double DampingFactor_
Damping parameter.
Teuchos::RefCountPtr< Epetra_Import > Importer_
Teuchos::RefCountPtr< Ifpack_Graph > Graph_
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual double ComputeTime() const
Returns the time spent in Compute().
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual int DoGaussSeidel(Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
virtual int Initialize()
Initializes the preconditioner.
virtual double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
int NumCompute_
Contains the number of successful call to Compute().
int NumLocalBlocks_
Number of local blocks.
bool ZeroStartingSolution_
If true, starting solution is the zero vector.
std::vector< Teuchos::RefCountPtr< T > > Containers_
virtual double Condest(const Ifpack_CondestType=Ifpack_Cheap, const int=1550, const double=1e-9, Epetra_RowMatrix *=0)
Computes the condition number estimate, returns its value.
virtual ~Ifpack_BlockRelaxation()
int NumLocalBlocks() const
Returns the number local blocks.
std::string Label_
Label for this object.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully computed.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Teuchos::RefCountPtr< Ifpack_Partitioner > Partitioner_
Contains information about non-overlapping partitions.
virtual const char * Label() const
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
double ComputeFlops_
Contains the number of flops for Compute().
int NumSweeps_
Number of preconditioning sweeps.
Teuchos::RefCountPtr< Epetra_Vector > W_
Weights for overlapping Jacobi only.
int NumInitialize_
Contains the number of successful calls to Initialize().
Teuchos::ParameterList List_
Parameters list to be used to solve on each subblock.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int ApplyInverseSGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Containers_[i] contains all the necessary information to solve on each subblock.
virtual int ApplyInverseJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int DoJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int Compute()
Computes the preconditioner.
std::string PartitionerType_
Ifpack_EquationPartitioner: A class to decompose an Ifpack_Graph so that each block will contain all ...
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
Ifpack_GreedyPartitioner: A class to decompose Ifpack_Graph's using a simple greedy algorithm.
Ifpack_LinearPartitioner: A class to define linear partitions.
Ifpack_METISPartitioner: A class to decompose Ifpack_Graph's using METIS.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
Ifpack_UserPartitioner: A class to define linear partitions.