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),
470 Matrix_(Teuchos::
rcp(Matrix_in,
false)),
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());
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);
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())
599 Matrix().RowMatrixRowMap()) );
604 ComputeTime_ += Time_.ElapsedTime();
619 if (X.NumVectors() != Y.NumVectors())
622 Time_.ResetStartTime();
626 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
627 if (X.Pointers()[0] == Y.Pointers()[0])
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;
839 Y.ExtractView(&y_ptr);
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;
966 Y.ExtractView(&y_ptr);
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];
1107 template<
typename T>
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;
1165 template<
typename T>
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="
1237 template<
typename T>
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();
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;
1307 #endif // IFPACK_BLOCKPRECONDITIONER_H
Teuchos::RefCountPtr< Epetra_Import > Importer_
int NumSweeps_
Number of preconditioning sweeps.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Containers_[i] contains all the necessary information to solve on each subblock.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool IsComputed_
If true, the preconditioner has been successfully computed.
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
double InitializeFlops_
Contains the number of flops for Initialize().
Teuchos::RefCountPtr< Ifpack_Graph > Graph_
Ifpack_UserPartitioner: A class to define linear partitions.
virtual const char * Label() const
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation &rhs)
operator= (PRIVATE, should not be used).
T & get(ParameterList &l, const std::string &name)
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
static const int IFPACK_SGS
virtual int ApplyInverseJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual int DoJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
Teuchos::ParameterList List_
Parameters list to be used to solve on each subblock.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
int NumLocalBlocks() const
Returns the number local blocks.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
std::string Label_
Label for this object.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully computed.
Ifpack_EquationPartitioner: A class to decompose an Ifpack_Graph so that each block will contain all ...
virtual int SetUseTranspose(bool UseTranspose_in)
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
double ComputeFlops_
Contains the number of flops for Compute().
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
double ComputeTime_
Contains the time for all successful calls to Compute().
int NumInitialize_
Contains the number of successful calls to Initialize().
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual ~Ifpack_BlockRelaxation()
virtual const Epetra_Comm & Comm() const =0
bool ZeroStartingSolution_
If true, starting solution is the zero vector.
virtual int Compute()
Computes the preconditioner.
std::string PartitionerType_
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Ifpack_METISPartitioner: A class to decompose Ifpack_Graph's using METIS.
static const int IFPACK_GS
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
Ifpack_BlockRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_BlockRelaxation constructor with given Epetra_RowMatrix.
Teuchos::RefCountPtr< Epetra_Vector > W_
Weights for overlapping Jacobi only.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual int DoSGS(const Epetra_MultiVector &X, Epetra_MultiVector &Xtmp, Epetra_MultiVector &Y) const
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int DoGaussSeidel(Epetra_MultiVector &X, Epetra_MultiVector &Y) const
#define IFPACK_RETURN(ifpack_err)
virtual double InitializeTime() const
Returns the time spent in Initialize().
std::vector< Teuchos::RefCountPtr< T > > Containers_
virtual int NumProc() const =0
virtual double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
static const int IFPACK_JACOBI
virtual int Initialize()
Initializes the preconditioner.
int NumLocalBlocks_
Number of local blocks.
virtual int ApplyInverseGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Ifpack_LinearPartitioner: A class to define linear partitions.
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
Teuchos::RefCountPtr< Ifpack_Partitioner > Partitioner_
Contains information about non-overlapping partitions.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
int NumCompute_
Contains the number of successful call to Compute().
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 double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
double DampingFactor_
Damping parameter.
#define IFPACK_CHK_ERR(ifpack_err)
Ifpack_GreedyPartitioner: A class to decompose Ifpack_Graph's using a simple greedy algorithm...
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
virtual int ApplyInverseSGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const