43#include "Ifpack_ConfigDefs.h"
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"
52#include "Ifpack_Preconditioner.h"
53#include "Ifpack_PointRelaxation.h"
55#include "Ifpack_Condest.h"
57static const int IFPACK_JACOBI = 0;
58static const int IFPACK_GS = 1;
59static const int IFPACK_SGS = 2;
64 IsInitialized_(false),
71 ApplyInverseTime_(0.0),
73 ApplyInverseFlops_(0.0),
79 PrecType_(IFPACK_JACOBI),
80 MinDiagonalValue_(0.0),
84 NumGlobalNonzeros_(0),
85 Matrix_(Teuchos::rcp(Matrix_in,false)),
87 ZeroStartingSolution_(true),
91 NumLocalSmoothingIndices_(Matrix_in->NumMyRows()),
92 LocalSmoothingIndices_(0)
103 if (PrecType_ == IFPACK_JACOBI)
105 else if (PrecType_ == IFPACK_GS)
107 else if (PrecType_ == IFPACK_SGS)
108 PT =
"symmetric Gauss-Seidel";
110 PT = List.get(
"relaxation: type", PT);
113 PrecType_ = IFPACK_JACOBI;
114 else if (PT ==
"Gauss-Seidel")
115 PrecType_ = IFPACK_GS;
116 else if (PT ==
"symmetric Gauss-Seidel")
117 PrecType_ = IFPACK_SGS;
122 NumSweeps_ = List.get(
"relaxation: sweeps",NumSweeps_);
123 DampingFactor_ = List.get(
"relaxation: damping factor",
125 MinDiagonalValue_ = List.get(
"relaxation: min diagonal value",
127 ZeroStartingSolution_ = List.get(
"relaxation: zero starting solution",
128 ZeroStartingSolution_);
130 DoBackwardGS_ = List.get(
"relaxation: backward mode",DoBackwardGS_);
132 DoL1Method_ = List.get(
"relaxation: use l1",DoL1Method_);
134 L1Eta_ = List.get(
"relaxation: l1 eta",L1Eta_);
139 NumLocalSmoothingIndices_= List.get(
"relaxation: number of local smoothing indices",NumLocalSmoothingIndices_);
140 LocalSmoothingIndices_ = List.get(
"relaxation: local smoothing indices",LocalSmoothingIndices_);
141 if(PrecType_ == IFPACK_JACOBI && LocalSmoothingIndices_) {
142 NumLocalSmoothingIndices_=NumMyRows_;
143 LocalSmoothingIndices_=0;
144 if(!
Comm().MyPID()) cout<<
"Ifpack_PointRelaxation: WARNING: Reordered/Local Smoothing not implemented for Jacobi. Defaulting to regular Jacobi"<<endl;
155 return(Matrix_->Comm());
161 return(Matrix_->OperatorDomainMap());
167 return(Matrix_->OperatorRangeMap());
173 IsInitialized_ =
false;
175 if (Matrix_ == Teuchos::null)
178 if (Time_ == Teuchos::null)
181 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
184 NumMyRows_ = Matrix_->NumMyRows();
185 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
186 NumGlobalRows_ = Matrix_->NumGlobalRows64();
187 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
189 if (
Comm().NumProc() != 1)
195 InitializeTime_ += Time_->ElapsedTime();
196 IsInitialized_ =
true;
207 Time_->ResetStartTime();
213 if (NumSweeps_ == 0) ierr = 1;
218#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
225 if (Diagonal_ == Teuchos::null)
228 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*Diagonal_));
239 if(DoL1Method_ && IsParallel_) {
241 std::vector<int> Indices(maxLength);
242 std::vector<double> Values(maxLength);
245 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
246 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, maxLength,NumEntries,
247 &Values[0], &Indices[0]));
248 double diagonal_boost=0.0;
249 for (
int k = 0 ; k < NumEntries ; ++k)
250 if(Indices[k] >= NumMyRows_)
251 diagonal_boost+=std::abs(Values[k]/2.0);
252 if ((*Diagonal_)[i] < L1Eta_*diagonal_boost)
253 (*Diagonal_)[i]+=diagonal_boost;
260 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
261 double& diag = (*Diagonal_)[i];
262 if (IFPACK_ABS(diag) < MinDiagonalValue_)
263 diag = MinDiagonalValue_;
267#ifdef IFPACK_FLOPCOUNTERS
268 ComputeFlops_ += NumMyRows_;
274 if ((PrecType_ == IFPACK_JACOBI) || (PrecType_ == IFPACK_GS)) {
275 Diagonal_->Reciprocal(*Diagonal_);
277 ComputeFlops_ += NumMyRows_;
287 if (IsParallel_ && ((PrecType_ == IFPACK_GS) || (PrecType_ == IFPACK_SGS))) {
289 Matrix().RowMatrixRowMap()) );
290 if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
294 ComputeTime_ += Time_->ElapsedTime();
297 IFPACK_CHK_ERR(ierr);
306 double MyMinVal, MyMaxVal;
307 double MinVal, MaxVal;
310 Diagonal_->MinValue(&MyMinVal);
311 Diagonal_->MaxValue(&MyMaxVal);
316 if (!
Comm().MyPID()) {
318 os <<
"================================================================================" << endl;
319 os <<
"Ifpack_PointRelaxation" << endl;
320 os <<
"Sweeps = " << NumSweeps_ << endl;
321 os <<
"damping factor = " << DampingFactor_ << endl;
322 if (PrecType_ == IFPACK_JACOBI)
323 os <<
"Type = Jacobi" << endl;
324 else if (PrecType_ == IFPACK_GS)
325 os <<
"Type = Gauss-Seidel" << endl;
326 else if (PrecType_ == IFPACK_SGS)
327 os <<
"Type = symmetric Gauss-Seidel" << endl;
329 os <<
"Using backward mode (GS only)" << endl;
330 if (ZeroStartingSolution_)
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;
343 os <<
"Initialize() " << std::setw(5) << NumInitialize_
344 <<
" " << std::setw(15) << InitializeTime_
345 <<
" 0.0 0.0" << endl;
346 os <<
"Compute() " << std::setw(5) << NumCompute_
347 <<
" " << std::setw(15) << ComputeTime_
348 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
349 if (ComputeTime_ != 0.0)
350 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
352 os <<
" " << std::setw(15) << 0.0 << endl;
353 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
354 <<
" " << std::setw(15) << ApplyInverseTime_
355 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
356 if (ApplyInverseTime_ != 0.0)
357 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
359 os <<
" " << std::setw(15) << 0.0 << endl;
360 os <<
"================================================================================" << endl;
369Condest(
const Ifpack_CondestType CT,
370 const int MaxIters,
const double Tol,
378 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
384void Ifpack_PointRelaxation::SetLabel()
387 if (PrecType_ == IFPACK_JACOBI)
389 else if (PrecType_ == IFPACK_GS){
392 PT =
"Backward " + PT;
394 else if (PrecType_ == IFPACK_SGS)
397 if(NumLocalSmoothingIndices_!=NumMyRows_) PT=
"Local " + PT;
398 else if(LocalSmoothingIndices_) PT=
"Reordered " + PT;
400 Label_ =
"IFPACK (" + PT +
", sweeps=" + Ifpack_toString(NumSweeps_)
401 +
", damping=" + Ifpack_toString(DampingFactor_) +
")";
422 Time_->ResetStartTime();
426 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
430 Xcopy = Teuchos::rcp( &X,
false );
435 IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
438 IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
441 IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
448 ApplyInverseTime_ += Time_->ElapsedTime();
456int Ifpack_PointRelaxation::
462 if (NumSweeps_ > 0 && ZeroStartingSolution_) {
465 for (
int v = 0; v < NumVectors; v++)
466 IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(RHS(v)), *Diagonal_, 0.0));
471 bool zeroOut =
false;
473 for (
int j = startIter; j < NumSweeps_ ; j++) {
474 IFPACK_CHK_ERR(
Apply(LHS, A_times_LHS));
475 IFPACK_CHK_ERR(A_times_LHS.Update(1.0,RHS,-1.0));
476 for (
int v = 0 ; v < NumVectors ; ++v)
477 IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(A_times_LHS(v)),
490#ifdef IFPACK_FLOPCOUNTERS
491 ApplyInverseFlops_ += NumVectors * (6 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
498int Ifpack_PointRelaxation::
501 if (ZeroStartingSolution_)
510 return(ApplyInverseGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
512 return(ApplyInverseGS_FastCrsMatrix(CrsMatrix, X, Y));
514 return(ApplyInverseGS_CrsMatrix(CrsMatrix, X, Y));
517 return(ApplyInverseGS_RowMatrix(X, Y));
523int Ifpack_PointRelaxation::
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);
543 Diagonal_->ExtractView(&d_ptr);
545 for (
int j = 0; j < NumSweeps_ ; j++) {
549 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
552 if (NumVectors == 1) {
554 double* y0_ptr = y_ptr[0];
555 double* y20_ptr = y2_ptr[0];
556 double* x0_ptr = x_ptr[0];
560 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
561 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
565 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
566 &Values[0], &Indices[0]));
569 for (
int k = 0 ; k < NumEntries ; ++k) {
572 dtemp += Values[k] * y20_ptr[col];
575 y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
580 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
581 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
586 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
587 &Values[0], &Indices[0]));
589 for (
int k = 0 ; k < NumEntries ; ++k) {
591 dtemp += Values[k] * y20_ptr[i];
594 y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
600 for (
int i = 0 ; i < NumMyRows_ ; ++i)
601 y0_ptr[i] = y20_ptr[i];
607 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
611 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
612 &Values[0], &Indices[0]));
614 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
629 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
632 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
633 &Values[0], &Indices[0]));
635 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
652 for (
int m = 0 ; m < NumVectors ; ++m)
653 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
654 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
655 y_ptr[m][i] = y2_ptr[m][i];
660#ifdef IFPACK_FLOPCOUNTERS
661 ApplyInverseFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
668int Ifpack_PointRelaxation::
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);
687 Diagonal_->ExtractView(&d_ptr);
689 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
693 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
697 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
698 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
702 double diag = d_ptr[i];
706 for (
int m = 0 ; m < NumVectors ; ++m) {
710 for (
int k = 0; k < NumEntries; ++k) {
713 dtemp += Values[k] * y2_ptr[m][col];
716 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
722 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
723 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
727 double diag = d_ptr[i];
731 for (
int m = 0 ; m < NumVectors ; ++m) {
734 for (
int k = 0; k < NumEntries; ++k) {
737 dtemp += Values[k] * y2_ptr[m][col];
740 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
747 for (
int m = 0 ; m < NumVectors ; ++m)
748 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
749 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
750 y_ptr[m][i] = y2_ptr[m][i];
754#ifdef IFPACK_FLOPCOUNTERS
755 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
764int Ifpack_PointRelaxation::
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);
785 Diagonal_->ExtractView(&d_ptr);
787 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
791 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
795 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
798 double diag = d_ptr[i];
800 for (
int m = 0 ; m < NumVectors ; ++m) {
804 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
807 dtemp += Values[k] * y2_ptr[m][col];
810 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
816 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
819 double diag = d_ptr[i];
821 for (
int m = 0 ; m < NumVectors ; ++m) {
824 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
827 dtemp += Values[k] * y2_ptr[m][col];
830 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
838 for (
int m = 0 ; m < NumVectors ; ++m)
839 for (
int i = 0 ; i < NumMyRows_ ; ++i)
840 y_ptr[m][i] = y2_ptr[m][i];
843#ifdef IFPACK_FLOPCOUNTERS
844 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
854int Ifpack_PointRelaxation::
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);
875 Diagonal_->ExtractView(&d_ptr);
877 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
881 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
885 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
886 int i=LocalSmoothingIndices_[ii];
889 double diag = d_ptr[i];
891 for (
int m = 0 ; m < NumVectors ; ++m) {
895 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
898 dtemp += Values[k] * y2_ptr[m][col];
901 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
907 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
908 int i=LocalSmoothingIndices_[ii];
911 double diag = d_ptr[i];
913 for (
int m = 0 ; m < NumVectors ; ++m) {
916 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
919 dtemp += Values[k] * y2_ptr[m][col];
922 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
930 for (
int m = 0 ; m < NumVectors ; ++m)
931 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
932 int i=LocalSmoothingIndices_[ii];
933 y_ptr[m][i] = y2_ptr[m][i];
937#ifdef IFPACK_FLOPCOUNTERS
938 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
945int Ifpack_PointRelaxation::
948 if (ZeroStartingSolution_)
957 return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
959 return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
961 return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
964 return(ApplyInverseSGS_RowMatrix(X, Y));
968int Ifpack_PointRelaxation::
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);
987 Diagonal_->ExtractView(&d_ptr);
989 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
993 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
995 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
996 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1000 double diag = d_ptr[i];
1002 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1003 &Values[0], &Indices[0]));
1005 for (
int m = 0 ; m < NumVectors ; ++m) {
1009 for (
int k = 0 ; k < NumEntries ; ++k) {
1012 dtemp += Values[k] * y2_ptr[m][col];
1015 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1018 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1019 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1023 double diag = d_ptr[i];
1025 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1026 &Values[0], &Indices[0]));
1028 for (
int m = 0 ; m < NumVectors ; ++m) {
1031 for (
int k = 0 ; k < NumEntries ; ++k) {
1034 dtemp += Values[k] * y2_ptr[m][col];
1037 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1043 for (
int m = 0 ; m < NumVectors ; ++m)
1044 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1045 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1046 y_ptr[m][i] = y2_ptr[m][i];
1050#ifdef IFPACK_FLOPCOUNTERS
1051 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1057int Ifpack_PointRelaxation::
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);
1076 Diagonal_->ExtractView(&d_ptr);
1078 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1082 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1084 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1085 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1089 double diag = d_ptr[i];
1093 for (
int m = 0 ; m < NumVectors ; ++m) {
1097 for (
int k = 0; k < NumEntries; ++k) {
1100 dtemp += Values[k] * y2_ptr[m][col];
1103 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1106 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1107 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1111 double diag = d_ptr[i];
1115 for (
int m = 0 ; m < NumVectors ; ++m) {
1118 for (
int k = 0; k < NumEntries; ++k) {
1121 dtemp += Values[k] * y2_ptr[m][col];
1124 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1130 for (
int m = 0 ; m < NumVectors ; ++m)
1131 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1132 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1133 y_ptr[m][i] = y2_ptr[m][i];
1137#ifdef IFPACK_FLOPCOUNTERS
1138 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1145int Ifpack_PointRelaxation::
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);
1166 Diagonal_->ExtractView(&d_ptr);
1168 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1172 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1174 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
1177 double diag = d_ptr[i];
1182 for (
int m = 0 ; m < NumVectors ; ++m) {
1186 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1189 dtemp += Values[k] * y2_ptr[m][col];
1192 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1196 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
1199 double diag = d_ptr[i];
1204 for (
int m = 0 ; m < NumVectors ; ++m) {
1207 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1210 dtemp += Values[k] * y2_ptr[m][col];
1213 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1219 for (
int m = 0 ; m < NumVectors ; ++m)
1220 for (
int i = 0 ; i < NumMyRows_ ; ++i)
1221 y_ptr[m][i] = y2_ptr[m][i];
1224#ifdef IFPACK_FLOPCOUNTERS
1225 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1233int Ifpack_PointRelaxation::
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);
1254 Diagonal_->ExtractView(&d_ptr);
1256 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1260 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1262 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1263 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1266 double diag = d_ptr[i];
1271 for (
int m = 0 ; m < NumVectors ; ++m) {
1275 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1278 dtemp += Values[k] * y2_ptr[m][col];
1281 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1285 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1286 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1289 double diag = d_ptr[i];
1294 for (
int m = 0 ; m < NumVectors ; ++m) {
1297 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1300 dtemp += Values[k] * y2_ptr[m][col];
1303 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1309 for (
int m = 0 ; m < NumVectors ; ++m)
1310 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1311 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1312 y_ptr[m][i] = y2_ptr[m][i];
1316#ifdef IFPACK_FLOPCOUNTERS
1317 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
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
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.
virtual int Compute()
Computes the preconditioners.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
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.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
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.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.