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)
181 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
189 if (
Comm().NumProc() != 1)
207 Time_->ResetStartTime();
218#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
241 std::vector<int> Indices(maxLength);
242 std::vector<double> Values(maxLength);
247 &Values[0], &Indices[0]));
248 double diagonal_boost=0.0;
249 for (
int k = 0 ; k < NumEntries ; ++k)
251 diagonal_boost+=std::abs(Values[k]/2.0);
253 (*Diagonal_)[i]+=diagonal_boost;
261 double& diag = (*Diagonal_)[i];
267#ifdef IFPACK_FLOPCOUNTERS
289 Matrix().RowMatrixRowMap()) );
306 double MyMinVal, MyMaxVal;
307 double MinVal, MaxVal;
316 if (!
Comm().MyPID()) {
318 os <<
"================================================================================" << endl;
319 os <<
"Ifpack_PointRelaxation" << endl;
323 os <<
"Type = Jacobi" << endl;
325 os <<
"Type = Gauss-Seidel" << endl;
327 os <<
"Type = symmetric Gauss-Seidel" << endl;
329 os <<
"Using backward mode (GS only)" << endl;
331 os <<
"Using zero starting solution" << endl;
333 os <<
"Using input starting solution" << endl;
334 os <<
"Condition number estimate = " <<
Condest() << endl;
335 os <<
"Global number of rows = " <<
Matrix_->NumGlobalRows64() << endl;
337 os <<
"Minimum value on stored diagonal = " << MinVal << endl;
338 os <<
"Maximum value on stored diagonal = " << MaxVal << endl;
341 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
342 os <<
"----- ------- -------------- ------------ --------" << endl;
345 <<
" 0.0 0.0" << endl;
352 os <<
" " << std::setw(15) << 0.0 << endl;
359 os <<
" " << std::setw(15) << 0.0 << endl;
360 os <<
"================================================================================" << endl;
370 const int MaxIters,
const double Tol,
392 PT =
"Backward " + PT;
422 Time_->ResetStartTime();
426 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
430 Xcopy = Teuchos::rcp( &X,
false );
471 bool zeroOut =
false;
473 for (
int j = startIter; j <
NumSweeps_ ; j++) {
490#ifdef IFPACK_FLOPCOUNTERS
529 std::vector<int> Indices(Length);
530 std::vector<double> Values(Length);
532 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
536 Y2 = Teuchos::rcp( &Y,
false );
539 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
542 Y2->ExtractView(&y2_ptr);
554 double* y0_ptr = y_ptr[0];
555 double* y20_ptr = y2_ptr[0];
556 double* x0_ptr = x_ptr[0];
566 &Values[0], &Indices[0]));
569 for (
int k = 0 ; k < NumEntries ; ++k) {
572 dtemp += Values[k] * y20_ptr[col];
587 &Values[0], &Indices[0]));
589 for (
int k = 0 ; k < NumEntries ; ++k) {
591 dtemp += Values[k] * y20_ptr[i];
601 y0_ptr[i] = y20_ptr[i];
612 &Values[0], &Indices[0]));
617 for (
int k = 0 ; k < NumEntries ; ++k) {
620 dtemp += Values[k] * y2_ptr[m][col];
623 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
633 &Values[0], &Indices[0]));
638 for (
int k = 0 ; k < NumEntries ; ++k) {
641 dtemp += Values[k] * y2_ptr[m][col];
644 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
655 y_ptr[m][i] = y2_ptr[m][i];
660#ifdef IFPACK_FLOPCOUNTERS
676 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
681 Y2 = Teuchos::rcp( &Y,
false );
683 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
686 Y2->ExtractView(&y2_ptr);
689 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
702 double diag = d_ptr[i];
710 for (
int k = 0; k < NumEntries; ++k) {
713 dtemp += Values[k] * y2_ptr[m][col];
727 double diag = d_ptr[i];
734 for (
int k = 0; k < NumEntries; ++k) {
737 dtemp += Values[k] * y2_ptr[m][col];
750 y_ptr[m][i] = y2_ptr[m][i];
754#ifdef IFPACK_FLOPCOUNTERS
774 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
779 Y2 = Teuchos::rcp( &Y,
false );
781 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
784 Y2->ExtractView(&y2_ptr);
787 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
798 double diag = d_ptr[i];
804 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
807 dtemp += Values[k] * y2_ptr[m][col];
819 double diag = d_ptr[i];
824 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
827 dtemp += Values[k] * y2_ptr[m][col];
840 y_ptr[m][i] = y2_ptr[m][i];
843#ifdef IFPACK_FLOPCOUNTERS
864 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
869 Y2 = Teuchos::rcp( &Y,
false );
871 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
874 Y2->ExtractView(&y2_ptr);
877 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
889 double diag = d_ptr[i];
895 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
898 dtemp += Values[k] * y2_ptr[m][col];
911 double diag = d_ptr[i];
916 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
919 dtemp += Values[k] * y2_ptr[m][col];
933 y_ptr[m][i] = y2_ptr[m][i];
937#ifdef IFPACK_FLOPCOUNTERS
973 std::vector<int> Indices(Length);
974 std::vector<double> Values(Length);
976 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
981 Y2 = Teuchos::rcp( &Y,
false );
983 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
986 Y2->ExtractView(&y2_ptr);
989 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1000 double diag = d_ptr[i];
1003 &Values[0], &Indices[0]));
1009 for (
int k = 0 ; k < NumEntries ; ++k) {
1012 dtemp += Values[k] * y2_ptr[m][col];
1023 double diag = d_ptr[i];
1026 &Values[0], &Indices[0]));
1031 for (
int k = 0 ; k < NumEntries ; ++k) {
1034 dtemp += Values[k] * y2_ptr[m][col];
1046 y_ptr[m][i] = y2_ptr[m][i];
1050#ifdef IFPACK_FLOPCOUNTERS
1065 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1070 Y2 = Teuchos::rcp( &Y,
false );
1072 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1075 Y2->ExtractView(&y2_ptr);
1078 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1089 double diag = d_ptr[i];
1097 for (
int k = 0; k < NumEntries; ++k) {
1100 dtemp += Values[k] * y2_ptr[m][col];
1111 double diag = d_ptr[i];
1118 for (
int k = 0; k < NumEntries; ++k) {
1121 dtemp += Values[k] * y2_ptr[m][col];
1133 y_ptr[m][i] = y2_ptr[m][i];
1137#ifdef IFPACK_FLOPCOUNTERS
1155 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1160 Y2 = Teuchos::rcp( &Y,
false );
1162 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1165 Y2->ExtractView(&y2_ptr);
1168 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1177 double diag = d_ptr[i];
1186 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1189 dtemp += Values[k] * y2_ptr[m][col];
1196 for (
int i =
NumMyRows_ - 1 ; i > -1 ; --i) {
1199 double diag = d_ptr[i];
1207 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1210 dtemp += Values[k] * y2_ptr[m][col];
1221 y_ptr[m][i] = y2_ptr[m][i];
1224#ifdef IFPACK_FLOPCOUNTERS
1243 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1248 Y2 = Teuchos::rcp( &Y,
false );
1250 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1253 Y2->ExtractView(&y2_ptr);
1256 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1266 double diag = d_ptr[i];
1275 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1278 dtemp += Values[k] * y2_ptr[m][col];
1289 double diag = d_ptr[i];
1297 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1300 dtemp += Values[k] * y2_ptr[m][col];
1312 y_ptr[m][i] = y2_ptr[m][i];
1316#ifdef IFPACK_FLOPCOUNTERS
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.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
static const int IFPACK_JACOBI
static const int IFPACK_GS
static const int IFPACK_SGS
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
bool StorageOptimized() const
const Epetra_BlockMap & Map() const
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 int MaxNumEntries() const=0
const Epetra_BlockMap & RowMap() const
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual int ApplyInverseJacobi(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the Jacobi preconditioner to X, returns the result in Y.
bool IsParallel_
If true, more than 1 processor is currently used.
virtual int ApplyInverseSGS_RowMatrix(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int Compute()
Computes the preconditioners.
bool DoL1Method_
Do L1 Jacobi/GS/SGS.
Teuchos::RefCountPtr< Epetra_Import > Importer_
Importer for parallel GS and SGS.
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual int ApplyInverseGS_CrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int NumMyNonzeros_
Number of local nonzeros.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int ApplyInverseSGS_CrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
long long NumGlobalRows_
Number of global rows.
bool DoBackwardGS_
Backward-Mode Gauss Seidel.
bool IsComputed_
If true, the preconditioner has been computed successfully.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
double DampingFactor_
Damping factor.
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 InitializeTime_
Contains the time for all successful calls to Initialize().
int NumLocalSmoothingIndices_
Number of (local) unknowns for local smoothing.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int NumSweeps_
Number of application of the preconditioner (should be greater than 0).
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
virtual int ApplyInverseGS_FastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
virtual int ApplyInverseGS_RowMatrix(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
int * LocalSmoothingIndices_
List of (local) unknowns for local smoothing (if any)
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
virtual void SetLabel()
Sets the label.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
virtual int ApplyInverseSGS_LocalFastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Teuchos::RefCountPtr< Epetra_Vector > Diagonal_
Contains the diagonal elements of Matrix.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
double ComputeFlops_
Contains the number of flops for Compute().
double Condest_
Contains the estimated condition number.
long long NumGlobalNonzeros_
Number of global nonzeros.
virtual int ApplyInverseSGS_FastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ApplyInverseGS_LocalFastCrsMatrix(const Epetra_CrsMatrix *A, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int ApplyInverseGS(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the Gauss-Seidel preconditioner to X, returns the result in Y.
std::string Label_
Contains the label of this object.
int NumMyRows_
Number of local rows.
int NumCompute_
Contains the number of successful call to Compute().
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool IsInitialized_
If true, the preconditioner has been computed successfully.
double L1Eta_
Eta parameter for modified L1 method.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.