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();
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 286 Matrix().RowMatrixRowMap()) );
303 double MyMinVal, MyMaxVal;
304 double MinVal, MaxVal;
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,
389 PT =
"Backward " + PT;
419 Time_->ResetStartTime();
423 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
468 bool zeroOut =
false;
470 for (
int j = startIter; j <
NumSweeps_ ; j++) {
487 #ifdef IFPACK_FLOPCOUNTERS 526 std::vector<int> Indices(Length);
527 std::vector<double> Values(Length);
529 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
536 double** y_ptr, ** y2_ptr, ** x_ptr, *d_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) {
588 dtemp += Values[k] * y20_ptr[i];
598 y0_ptr[i] = y20_ptr[i];
609 &Values[0], &Indices[0]));
614 for (
int k = 0 ; k < NumEntries ; ++k) {
617 dtemp += Values[k] * y2_ptr[m][col];
620 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
630 &Values[0], &Indices[0]));
635 for (
int k = 0 ; k < NumEntries ; ++k) {
638 dtemp += Values[k] * y2_ptr[m][col];
641 y2_ptr[m][i] +=
DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
652 y_ptr[m][i] = y2_ptr[m][i];
657 #ifdef IFPACK_FLOPCOUNTERS 673 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
680 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
683 Y2->ExtractView(&y2_ptr);
686 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
699 double diag = d_ptr[i];
707 for (
int k = 0; k < NumEntries; ++k) {
710 dtemp += Values[k] * y2_ptr[m][col];
724 double diag = d_ptr[i];
731 for (
int k = 0; k < NumEntries; ++k) {
734 dtemp += Values[k] * y2_ptr[m][col];
747 y_ptr[m][i] = y2_ptr[m][i];
751 #ifdef IFPACK_FLOPCOUNTERS 767 IFPACK_CHK_ERR(
A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
771 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
778 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
781 Y2->ExtractView(&y2_ptr);
784 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
795 double diag = d_ptr[i];
801 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
804 dtemp += Values[k] * y2_ptr[m][col];
816 double diag = d_ptr[i];
821 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
824 dtemp += Values[k] * y2_ptr[m][col];
837 y_ptr[m][i] = y2_ptr[m][i];
840 #ifdef IFPACK_FLOPCOUNTERS 857 IFPACK_CHK_ERR(
A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
861 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
868 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
871 Y2->ExtractView(&y2_ptr);
874 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
886 double diag = d_ptr[i];
892 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
895 dtemp += Values[k] * y2_ptr[m][col];
908 double diag = d_ptr[i];
913 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
916 dtemp += Values[k] * y2_ptr[m][col];
930 y_ptr[m][i] = y2_ptr[m][i];
934 #ifdef IFPACK_FLOPCOUNTERS 970 std::vector<int> Indices(Length);
971 std::vector<double> Values(Length);
973 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
980 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
983 Y2->ExtractView(&y2_ptr);
986 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
997 double diag = d_ptr[i];
1000 &Values[0], &Indices[0]));
1006 for (
int k = 0 ; k < NumEntries ; ++k) {
1009 dtemp += Values[k] * y2_ptr[m][col];
1020 double diag = d_ptr[i];
1023 &Values[0], &Indices[0]));
1028 for (
int k = 0 ; k < NumEntries ; ++k) {
1031 dtemp += Values[k] * y2_ptr[m][col];
1043 y_ptr[m][i] = y2_ptr[m][i];
1047 #ifdef IFPACK_FLOPCOUNTERS 1062 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1069 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1072 Y2->ExtractView(&y2_ptr);
1075 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1086 double diag = d_ptr[i];
1094 for (
int k = 0; k < NumEntries; ++k) {
1097 dtemp += Values[k] * y2_ptr[m][col];
1108 double diag = d_ptr[i];
1115 for (
int k = 0; k < NumEntries; ++k) {
1118 dtemp += Values[k] * y2_ptr[m][col];
1130 y_ptr[m][i] = y2_ptr[m][i];
1134 #ifdef IFPACK_FLOPCOUNTERS 1148 IFPACK_CHK_ERR(
A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1152 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1159 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1162 Y2->ExtractView(&y2_ptr);
1165 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1174 double diag = d_ptr[i];
1183 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1186 dtemp += Values[k] * y2_ptr[m][col];
1193 for (
int i =
NumMyRows_ - 1 ; i > -1 ; --i) {
1196 double diag = d_ptr[i];
1204 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1207 dtemp += Values[k] * y2_ptr[m][col];
1218 y_ptr[m][i] = y2_ptr[m][i];
1221 #ifdef IFPACK_FLOPCOUNTERS 1236 IFPACK_CHK_ERR(
A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1240 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1247 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1250 Y2->ExtractView(&y2_ptr);
1253 for (
int iter = 0 ; iter <
NumSweeps_ ; ++iter) {
1263 double diag = d_ptr[i];
1272 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1275 dtemp += Values[k] * y2_ptr[m][col];
1286 double diag = d_ptr[i];
1294 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1297 dtemp += Values[k] * y2_ptr[m][col];
1309 y_ptr[m][i] = y2_ptr[m][i];
1313 #ifdef IFPACK_FLOPCOUNTERS bool StorageOptimized() const
int NumCompute_
Contains the number of successful call to Compute().
const Epetra_BlockMap & Map() const
std::string Label_
Contains the label of this object.
const Epetra_BlockMap & RowMap() const
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
T & get(ParameterList &l, const std::string &name)
int PutScalar(double ScalarConstant)
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 MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
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.
virtual int MaxNumEntries() const=0
static const int IFPACK_JACOBI
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
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 ExtractView(double **A, int *MyLDA) const
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.
double ** Pointers() const
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