46 #include "Epetra_Operator.h" 47 #include "Epetra_CrsMatrix.h" 48 #include "Epetra_VbrMatrix.h" 49 #include "Epetra_Comm.h" 50 #include "Epetra_Map.h" 51 #include "Epetra_MultiVector.h" 64 IsInitialized_(
false),
71 ApplyInverseTime_(0.0),
73 ApplyInverseFlops_(0.0),
80 MinDiagonalValue_(0.0),
84 NumGlobalNonzeros_(0),
87 ZeroStartingSolution_(
true),
91 NumLocalSmoothingIndices_(Matrix_in->NumMyRows()),
92 LocalSmoothingIndices_(0)
108 PT =
"symmetric Gauss-Seidel";
110 PT = List.get(
"relaxation: type", PT);
114 else if (PT ==
"Gauss-Seidel")
116 else if (PT ==
"symmetric Gauss-Seidel")
144 if(!
Comm().MyPID()) cout<<
"Ifpack_PointRelaxation: WARNING: Reordered/Local Smoothing not implemented for Jacobi. Defaulting to regular Jacobi"<<endl;
161 return(
Matrix_->OperatorDomainMap());
167 return(
Matrix_->OperatorRangeMap());
178 if (
Time_ == Teuchos::null)
179 Time_ = Teuchos::rcp(
new Epetra_Time(
Comm()) );
181 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
189 if (
Comm().NumProc() != 1)
207 Time_->ResetStartTime();
218 const Epetra_VbrMatrix * VbrMat =
dynamic_cast<const Epetra_VbrMatrix*
>(&
Matrix());
219 if(VbrMat)
Diagonal_ = Teuchos::rcp(
new Epetra_Vector(VbrMat->RowMap()) );
220 else Diagonal_ = Teuchos::rcp(
new Epetra_Vector(
Matrix().RowMatrixRowMap()) );
237 int maxLength =
Matrix().MaxNumEntries();
238 std::vector<int> Indices(maxLength);
239 std::vector<double> Values(maxLength);
244 &Values[0], &Indices[0]));
245 double diagonal_boost=0.0;
246 for (
int k = 0 ; k < NumEntries ; ++k)
248 diagonal_boost+=std::abs(Values[k]/2.0);
250 (*Diagonal_)[i]+=diagonal_boost;
258 double& diag = (*Diagonal_)[i];
264 #ifdef IFPACK_FLOPCOUNTERS 285 Importer_ = Teuchos::rcp(
new Epetra_Import(
Matrix().RowMatrixColMap(),
286 Matrix().RowMatrixRowMap()) );
303 double MyMinVal, MyMaxVal;
304 double MinVal, MaxVal;
309 Comm().MinAll(&MyMinVal,&MinVal,1);
310 Comm().MinAll(&MyMaxVal,&MaxVal,1);
313 if (!
Comm().MyPID()) {
315 os <<
"================================================================================" << endl;
316 os <<
"Ifpack_PointRelaxation" << endl;
320 os <<
"Type = Jacobi" << endl;
322 os <<
"Type = Gauss-Seidel" << endl;
324 os <<
"Type = symmetric Gauss-Seidel" << endl;
326 os <<
"Using backward mode (GS only)" << endl;
328 os <<
"Using zero starting solution" << endl;
330 os <<
"Using input starting solution" << endl;
331 os <<
"Condition number estimate = " <<
Condest() << endl;
332 os <<
"Global number of rows = " <<
Matrix_->NumGlobalRows64() << endl;
334 os <<
"Minimum value on stored diagonal = " << MinVal << endl;
335 os <<
"Maximum value on stored diagonal = " << MaxVal << endl;
338 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
339 os <<
"----- ------- -------------- ------------ --------" << endl;
342 <<
" 0.0 0.0" << endl;
349 os <<
" " << std::setw(15) << 0.0 << endl;
356 os <<
" " << std::setw(15) << 0.0 << endl;
357 os <<
"================================================================================" << endl;
367 const int MaxIters,
const double Tol,
368 Epetra_RowMatrix* Matrix_in)
389 PT =
"Backward " + PT;
416 if (X.NumVectors() != Y.NumVectors())
419 Time_->ResetStartTime();
423 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
424 if (X.Pointers()[0] == Y.Pointers()[0])
425 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
427 Xcopy = Teuchos::rcp( &X,
false );
468 bool zeroOut =
false;
469 Epetra_MultiVector A_times_LHS(LHS.Map(),
NumVectors, zeroOut);
470 for (
int j = startIter; j <
NumSweeps_ ; j++) {
487 #ifdef IFPACK_FLOPCOUNTERS 501 const Epetra_CrsMatrix* CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*
Matrix_);
508 else if (CrsMatrix->StorageOptimized())
525 int Length =
Matrix().MaxNumEntries();
526 std::vector<int> Indices(Length);
527 std::vector<double> Values(Length);
529 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
533 Y2 = Teuchos::rcp( &Y,
false );
536 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
537 X.ExtractView(&x_ptr);
538 Y.ExtractView(&y_ptr);
539 Y2->ExtractView(&y2_ptr);
551 double* y0_ptr = y_ptr[0];
552 double* y20_ptr = y2_ptr[0];
553 double* x0_ptr = x_ptr[0];
563 &Values[0], &Indices[0]));
566 for (
int k = 0 ; k < NumEntries ; ++k) {
569 dtemp += Values[k] * y20_ptr[col];
584 &Values[0], &Indices[0]));
586 for (
int k = 0 ; k < NumEntries ; ++k) {
589 dtemp += Values[k] * y20_ptr[i];
599 y0_ptr[i] = y20_ptr[i];
610 &Values[0], &Indices[0]));
615 for (
int k = 0 ; k < NumEntries ; ++k) {
618 dtemp += Values[k] * y2_ptr[m][col];
621 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
631 &Values[0], &Indices[0]));
636 for (
int k = 0 ; k < NumEntries ; ++k) {
639 dtemp += Values[k] * y2_ptr[m][col];
642 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
653 y_ptr[m][i] = y2_ptr[m][i];
658 #ifdef IFPACK_FLOPCOUNTERS 674 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
679 Y2 = Teuchos::rcp( &Y,
false );
681 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
682 X.ExtractView(&x_ptr);
683 Y.ExtractView(&y_ptr);
684 Y2->ExtractView(&y2_ptr);
687 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
700 double diag = d_ptr[i];
702 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
708 for (
int k = 0; k < NumEntries; ++k) {
711 dtemp += Values[k] * y2_ptr[m][col];
725 double diag = d_ptr[i];
727 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
732 for (
int k = 0; k < NumEntries; ++k) {
735 dtemp += Values[k] * y2_ptr[m][col];
748 y_ptr[m][i] = y2_ptr[m][i];
752 #ifdef IFPACK_FLOPCOUNTERS 768 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
772 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
777 Y2 = Teuchos::rcp( &Y,
false );
779 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
780 X.ExtractView(&x_ptr);
781 Y.ExtractView(&y_ptr);
782 Y2->ExtractView(&y2_ptr);
785 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
796 double diag = d_ptr[i];
802 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
805 dtemp += Values[k] * y2_ptr[m][col];
817 double diag = d_ptr[i];
822 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
825 dtemp += Values[k] * y2_ptr[m][col];
838 y_ptr[m][i] = y2_ptr[m][i];
841 #ifdef IFPACK_FLOPCOUNTERS 858 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
862 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
867 Y2 = Teuchos::rcp( &Y,
false );
869 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
870 X.ExtractView(&x_ptr);
871 Y.ExtractView(&y_ptr);
872 Y2->ExtractView(&y2_ptr);
875 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
887 double diag = d_ptr[i];
893 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
896 dtemp += Values[k] * y2_ptr[m][col];
909 double diag = d_ptr[i];
914 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
917 dtemp += Values[k] * y2_ptr[m][col];
931 y_ptr[m][i] = y2_ptr[m][i];
935 #ifdef IFPACK_FLOPCOUNTERS 949 const Epetra_CrsMatrix* CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*
Matrix_);
956 else if (CrsMatrix->StorageOptimized())
970 int Length =
Matrix().MaxNumEntries();
971 std::vector<int> Indices(Length);
972 std::vector<double> Values(Length);
974 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
979 Y2 = Teuchos::rcp( &Y,
false );
981 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
982 X.ExtractView(&x_ptr);
983 Y.ExtractView(&y_ptr);
984 Y2->ExtractView(&y2_ptr);
987 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
998 double diag = d_ptr[i];
1001 &Values[0], &Indices[0]));
1007 for (
int k = 0 ; k < NumEntries ; ++k) {
1010 dtemp += Values[k] * y2_ptr[m][col];
1021 double diag = d_ptr[i];
1024 &Values[0], &Indices[0]));
1029 for (
int k = 0 ; k < NumEntries ; ++k) {
1032 dtemp += Values[k] * y2_ptr[m][col];
1044 y_ptr[m][i] = y2_ptr[m][i];
1048 #ifdef IFPACK_FLOPCOUNTERS 1063 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1068 Y2 = Teuchos::rcp( &Y,
false );
1070 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1071 X.ExtractView(&x_ptr);
1072 Y.ExtractView(&y_ptr);
1073 Y2->ExtractView(&y2_ptr);
1076 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1087 double diag = d_ptr[i];
1089 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1095 for (
int k = 0; k < NumEntries; ++k) {
1098 dtemp += Values[k] * y2_ptr[m][col];
1109 double diag = d_ptr[i];
1111 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1116 for (
int k = 0; k < NumEntries; ++k) {
1119 dtemp += Values[k] * y2_ptr[m][col];
1131 y_ptr[m][i] = y2_ptr[m][i];
1135 #ifdef IFPACK_FLOPCOUNTERS 1149 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1153 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1158 Y2 = Teuchos::rcp( &Y,
false );
1160 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1161 X.ExtractView(&x_ptr);
1162 Y.ExtractView(&y_ptr);
1163 Y2->ExtractView(&y2_ptr);
1166 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1175 double diag = d_ptr[i];
1184 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1187 dtemp += Values[k] * y2_ptr[m][col];
1194 for (
int i =
NumMyRows_ - 1 ; i > -1 ; --i) {
1197 double diag = d_ptr[i];
1205 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1208 dtemp += Values[k] * y2_ptr[m][col];
1219 y_ptr[m][i] = y2_ptr[m][i];
1222 #ifdef IFPACK_FLOPCOUNTERS 1237 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1241 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1246 Y2 = Teuchos::rcp( &Y,
false );
1248 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1249 X.ExtractView(&x_ptr);
1250 Y.ExtractView(&y_ptr);
1251 Y2->ExtractView(&y2_ptr);
1254 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1264 double diag = d_ptr[i];
1273 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1276 dtemp += Values[k] * y2_ptr[m][col];
1287 double diag = d_ptr[i];
1295 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1298 dtemp += Values[k] * y2_ptr[m][col];
1310 y_ptr[m][i] = y2_ptr[m][i];
1314 #ifdef IFPACK_FLOPCOUNTERS int NumCompute_
Contains the number of successful call to Compute().
std::string Label_
Contains the label of this object.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
bool IsInitialized_
If true, the preconditioner has been computed successfully.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int ApplyInverseGS_LocalFastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
static const int IFPACK_GS
virtual int ApplyInverseSGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the symmetric Gauss-Seidel preconditioner to X, returns the result in Y.
double ComputeFlops_
Contains the number of flops for Compute().
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual int ApplyInverseSGS_CrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
virtual int ApplyInverseGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the Gauss-Seidel preconditioner to X, returns the result in Y.
virtual int ApplyInverseSGS_RowMatrix(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
static const int IFPACK_SGS
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
long long NumGlobalNonzeros_
Number of global nonzeros.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
Teuchos::RefCountPtr< Epetra_Vector > Diagonal_
Contains the diagonal elements of Matrix.
int NumLocalSmoothingIndices_
Number of (local) unknowns for local smoothing.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
static const int IFPACK_JACOBI
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
long long NumGlobalRows_
Number of global rows.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
double ComputeTime_
Contains the time for all successful calls to Compute().
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int ApplyInverseJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the Jacobi preconditioner to X, returns the result in Y.
double L1Eta_
Eta parameter for modified L1 method.
virtual int Compute()
Computes the preconditioners.
int NumInitialize_
Contains the number of successful calls to Initialize().
bool IsComputed_
If true, the preconditioner has been computed successfully.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
int NumMyRows_
Number of local rows.
bool IsParallel_
If true, more than 1 processor is currently used.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
double DampingFactor_
Damping factor.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual int ApplyInverseSGS_FastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
bool DoBackwardGS_
Backward-Mode Gauss Seidel.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
int NumSweeps_
Number of application of the preconditioner (should be greater than 0).
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
int * LocalSmoothingIndices_
List of (local) unknowns for local smoothing (if any)
virtual int ApplyInverseGS_FastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual void SetLabel()
Sets the label.
#define IFPACK_CHK_ERR(ifpack_err)
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double Condest_
Contains the estimated condition number.
int NumMyNonzeros_
Number of local nonzeros.
virtual int ApplyInverseGS_CrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::RefCountPtr< Epetra_Import > Importer_
Importer for parallel GS and SGS.
virtual int ApplyInverseGS_RowMatrix(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
bool DoL1Method_
Do L1 Jacobi/GS/SGS.
virtual int ApplyInverseSGS_LocalFastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const