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" 57 static const int IFPACK_JACOBI = 0;
58 static const int IFPACK_GS = 1;
59 static 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;
222 if (Diagonal_ == Teuchos::null)
225 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*Diagonal_));
236 if(DoL1Method_ && IsParallel_) {
238 std::vector<int> Indices(maxLength);
239 std::vector<double> Values(maxLength);
242 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
243 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, maxLength,NumEntries,
244 &Values[0], &Indices[0]));
245 double diagonal_boost=0.0;
246 for (
int k = 0 ; k < NumEntries ; ++k)
247 if(Indices[k] > NumMyRows_)
248 diagonal_boost+=std::abs(Values[k]/2.0);
249 if ((*Diagonal_)[i] < L1Eta_*diagonal_boost)
250 (*Diagonal_)[i]+=diagonal_boost;
257 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
258 double& diag = (*Diagonal_)[i];
259 if (IFPACK_ABS(diag) < MinDiagonalValue_)
260 diag = MinDiagonalValue_;
264 #ifdef IFPACK_FLOPCOUNTERS 265 ComputeFlops_ += NumMyRows_;
271 if ((PrecType_ == IFPACK_JACOBI) || (PrecType_ == IFPACK_GS)) {
272 Diagonal_->Reciprocal(*Diagonal_);
274 ComputeFlops_ += NumMyRows_;
284 if (IsParallel_ && ((PrecType_ == IFPACK_GS) || (PrecType_ == IFPACK_SGS))) {
286 Matrix().RowMatrixRowMap()) );
287 if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
291 ComputeTime_ += Time_->ElapsedTime();
294 IFPACK_CHK_ERR(ierr);
303 double MyMinVal, MyMaxVal;
304 double MinVal, MaxVal;
307 Diagonal_->MinValue(&MyMinVal);
308 Diagonal_->MaxValue(&MyMaxVal);
313 if (!
Comm().MyPID()) {
315 os <<
"================================================================================" << endl;
316 os <<
"Ifpack_PointRelaxation" << endl;
317 os <<
"Sweeps = " << NumSweeps_ << endl;
318 os <<
"damping factor = " << DampingFactor_ << endl;
319 if (PrecType_ == IFPACK_JACOBI)
320 os <<
"Type = Jacobi" << endl;
321 else if (PrecType_ == IFPACK_GS)
322 os <<
"Type = Gauss-Seidel" << endl;
323 else if (PrecType_ == IFPACK_SGS)
324 os <<
"Type = symmetric Gauss-Seidel" << endl;
326 os <<
"Using backward mode (GS only)" << endl;
327 if (ZeroStartingSolution_)
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;
340 os <<
"Initialize() " << std::setw(5) << NumInitialize_
341 <<
" " << std::setw(15) << InitializeTime_
342 <<
" 0.0 0.0" << endl;
343 os <<
"Compute() " << std::setw(5) << NumCompute_
344 <<
" " << std::setw(15) << ComputeTime_
345 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
346 if (ComputeTime_ != 0.0)
347 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
349 os <<
" " << std::setw(15) << 0.0 << endl;
350 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
351 <<
" " << std::setw(15) << ApplyInverseTime_
352 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
353 if (ApplyInverseTime_ != 0.0)
354 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
356 os <<
" " << std::setw(15) << 0.0 << endl;
357 os <<
"================================================================================" << endl;
367 const int MaxIters,
const double Tol,
375 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
381 void Ifpack_PointRelaxation::SetLabel()
384 if (PrecType_ == IFPACK_JACOBI)
386 else if (PrecType_ == IFPACK_GS){
389 PT =
"Backward " + PT;
391 else if (PrecType_ == IFPACK_SGS)
394 if(NumLocalSmoothingIndices_!=NumMyRows_) PT=
"Local " + PT;
395 else if(LocalSmoothingIndices_) PT=
"Reordered " + PT;
419 Time_->ResetStartTime();
423 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
427 Xcopy = Teuchos::rcp( &X,
false );
432 IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
435 IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
438 IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
445 ApplyInverseTime_ += Time_->ElapsedTime();
453 int Ifpack_PointRelaxation::
459 if (NumSweeps_ > 0 && ZeroStartingSolution_) {
462 for (
int v = 0; v < NumVectors; v++)
463 IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(RHS(v)), *Diagonal_, 0.0));
468 bool zeroOut =
false;
470 for (
int j = startIter; j < NumSweeps_ ; j++) {
471 IFPACK_CHK_ERR(
Apply(LHS, A_times_LHS));
472 IFPACK_CHK_ERR(A_times_LHS.Update(1.0,RHS,-1.0));
473 for (
int v = 0 ; v < NumVectors ; ++v)
474 IFPACK_CHK_ERR(LHS(v)->
Multiply(DampingFactor_, *(A_times_LHS(v)),
487 #ifdef IFPACK_FLOPCOUNTERS 488 ApplyInverseFlops_ += NumVectors * (6 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
495 int Ifpack_PointRelaxation::
498 if (ZeroStartingSolution_)
507 return(ApplyInverseGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
509 return(ApplyInverseGS_FastCrsMatrix(CrsMatrix, X, Y));
511 return(ApplyInverseGS_CrsMatrix(CrsMatrix, X, Y));
514 return(ApplyInverseGS_RowMatrix(X, Y));
520 int Ifpack_PointRelaxation::
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;
539 Y2->ExtractView(&y2_ptr);
540 Diagonal_->ExtractView(&d_ptr);
542 for (
int j = 0; j < NumSweeps_ ; j++) {
546 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
549 if (NumVectors == 1) {
551 double* y0_ptr = y_ptr[0];
552 double* y20_ptr = y2_ptr[0];
553 double* x0_ptr = x_ptr[0];
557 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
558 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
562 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
563 &Values[0], &Indices[0]));
566 for (
int k = 0 ; k < NumEntries ; ++k) {
569 dtemp += Values[k] * y20_ptr[col];
572 y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
577 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
578 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
583 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
584 &Values[0], &Indices[0]));
586 for (
int k = 0 ; k < NumEntries ; ++k) {
589 dtemp += Values[k] * y20_ptr[i];
592 y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
598 for (
int i = 0 ; i < NumMyRows_ ; ++i)
599 y0_ptr[i] = y20_ptr[i];
605 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
609 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
610 &Values[0], &Indices[0]));
612 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
627 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
630 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
631 &Values[0], &Indices[0]));
633 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
650 for (
int m = 0 ; m < NumVectors ; ++m)
651 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
652 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
653 y_ptr[m][i] = y2_ptr[m][i];
658 #ifdef IFPACK_FLOPCOUNTERS 659 ApplyInverseFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
666 int Ifpack_PointRelaxation::
674 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
679 Y2 = Teuchos::rcp( &Y,
false );
681 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
684 Y2->ExtractView(&y2_ptr);
685 Diagonal_->ExtractView(&d_ptr);
687 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
691 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
695 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
696 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
700 double diag = d_ptr[i];
704 for (
int m = 0 ; m < NumVectors ; ++m) {
708 for (
int k = 0; k < NumEntries; ++k) {
711 dtemp += Values[k] * y2_ptr[m][col];
714 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
720 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
721 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
725 double diag = d_ptr[i];
729 for (
int m = 0 ; m < NumVectors ; ++m) {
732 for (
int k = 0; k < NumEntries; ++k) {
735 dtemp += Values[k] * y2_ptr[m][col];
738 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
745 for (
int m = 0 ; m < NumVectors ; ++m)
746 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
747 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
748 y_ptr[m][i] = y2_ptr[m][i];
752 #ifdef IFPACK_FLOPCOUNTERS 753 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
762 int Ifpack_PointRelaxation::
772 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
777 Y2 = Teuchos::rcp( &Y,
false );
779 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
782 Y2->ExtractView(&y2_ptr);
783 Diagonal_->ExtractView(&d_ptr);
785 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
789 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
793 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
796 double diag = d_ptr[i];
798 for (
int m = 0 ; m < NumVectors ; ++m) {
802 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
805 dtemp += Values[k] * y2_ptr[m][col];
808 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
814 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
817 double diag = d_ptr[i];
819 for (
int m = 0 ; m < NumVectors ; ++m) {
822 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
825 dtemp += Values[k] * y2_ptr[m][col];
828 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
836 for (
int m = 0 ; m < NumVectors ; ++m)
837 for (
int i = 0 ; i < NumMyRows_ ; ++i)
838 y_ptr[m][i] = y2_ptr[m][i];
841 #ifdef IFPACK_FLOPCOUNTERS 842 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
852 int Ifpack_PointRelaxation::
862 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
867 Y2 = Teuchos::rcp( &Y,
false );
869 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
872 Y2->ExtractView(&y2_ptr);
873 Diagonal_->ExtractView(&d_ptr);
875 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
879 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
883 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
884 int i=LocalSmoothingIndices_[ii];
887 double diag = d_ptr[i];
889 for (
int m = 0 ; m < NumVectors ; ++m) {
893 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
896 dtemp += Values[k] * y2_ptr[m][col];
899 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
905 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
906 int i=LocalSmoothingIndices_[ii];
909 double diag = d_ptr[i];
911 for (
int m = 0 ; m < NumVectors ; ++m) {
914 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
917 dtemp += Values[k] * y2_ptr[m][col];
920 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
928 for (
int m = 0 ; m < NumVectors ; ++m)
929 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
930 int i=LocalSmoothingIndices_[ii];
931 y_ptr[m][i] = y2_ptr[m][i];
935 #ifdef IFPACK_FLOPCOUNTERS 936 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
943 int Ifpack_PointRelaxation::
946 if (ZeroStartingSolution_)
955 return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
957 return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
959 return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
962 return(ApplyInverseSGS_RowMatrix(X, Y));
966 int Ifpack_PointRelaxation::
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;
984 Y2->ExtractView(&y2_ptr);
985 Diagonal_->ExtractView(&d_ptr);
987 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
991 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
993 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
994 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
998 double diag = d_ptr[i];
1000 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1001 &Values[0], &Indices[0]));
1003 for (
int m = 0 ; m < NumVectors ; ++m) {
1007 for (
int k = 0 ; k < NumEntries ; ++k) {
1010 dtemp += Values[k] * y2_ptr[m][col];
1013 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1016 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1017 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1021 double diag = d_ptr[i];
1023 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1024 &Values[0], &Indices[0]));
1026 for (
int m = 0 ; m < NumVectors ; ++m) {
1029 for (
int k = 0 ; k < NumEntries ; ++k) {
1032 dtemp += Values[k] * y2_ptr[m][col];
1035 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1041 for (
int m = 0 ; m < NumVectors ; ++m)
1042 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1043 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1044 y_ptr[m][i] = y2_ptr[m][i];
1048 #ifdef IFPACK_FLOPCOUNTERS 1049 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1055 int Ifpack_PointRelaxation::
1063 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1068 Y2 = Teuchos::rcp( &Y,
false );
1070 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1073 Y2->ExtractView(&y2_ptr);
1074 Diagonal_->ExtractView(&d_ptr);
1076 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1080 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1082 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1083 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1087 double diag = d_ptr[i];
1091 for (
int m = 0 ; m < NumVectors ; ++m) {
1095 for (
int k = 0; k < NumEntries; ++k) {
1098 dtemp += Values[k] * y2_ptr[m][col];
1101 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1104 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1105 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1109 double diag = d_ptr[i];
1113 for (
int m = 0 ; m < NumVectors ; ++m) {
1116 for (
int k = 0; k < NumEntries; ++k) {
1119 dtemp += Values[k] * y2_ptr[m][col];
1122 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1128 for (
int m = 0 ; m < NumVectors ; ++m)
1129 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1130 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1131 y_ptr[m][i] = y2_ptr[m][i];
1135 #ifdef IFPACK_FLOPCOUNTERS 1136 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1143 int Ifpack_PointRelaxation::
1153 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1158 Y2 = Teuchos::rcp( &Y,
false );
1160 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1163 Y2->ExtractView(&y2_ptr);
1164 Diagonal_->ExtractView(&d_ptr);
1166 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1170 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1172 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
1175 double diag = d_ptr[i];
1180 for (
int m = 0 ; m < NumVectors ; ++m) {
1184 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1187 dtemp += Values[k] * y2_ptr[m][col];
1190 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1194 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
1197 double diag = d_ptr[i];
1202 for (
int m = 0 ; m < NumVectors ; ++m) {
1205 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1208 dtemp += Values[k] * y2_ptr[m][col];
1211 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1217 for (
int m = 0 ; m < NumVectors ; ++m)
1218 for (
int i = 0 ; i < NumMyRows_ ; ++i)
1219 y_ptr[m][i] = y2_ptr[m][i];
1222 #ifdef IFPACK_FLOPCOUNTERS 1223 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1231 int Ifpack_PointRelaxation::
1241 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1246 Y2 = Teuchos::rcp( &Y,
false );
1248 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1251 Y2->ExtractView(&y2_ptr);
1252 Diagonal_->ExtractView(&d_ptr);
1254 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1258 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1260 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1261 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1264 double diag = d_ptr[i];
1269 for (
int m = 0 ; m < NumVectors ; ++m) {
1273 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1276 dtemp += Values[k] * y2_ptr[m][col];
1279 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1283 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1284 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1287 double diag = d_ptr[i];
1292 for (
int m = 0 ; m < NumVectors ; ++m) {
1295 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1298 dtemp += Values[k] * y2_ptr[m][col];
1301 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1307 for (
int m = 0 ; m < NumVectors ; ++m)
1308 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1309 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1310 y_ptr[m][i] = y2_ptr[m][i];
1314 #ifdef IFPACK_FLOPCOUNTERS 1315 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
bool StorageOptimized() const
const Epetra_BlockMap & Map() const
const Epetra_BlockMap & RowMap() const
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
int PutScalar(double ScalarConstant)
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
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual int MaxNumEntries() const=0
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int Compute()
Computes the preconditioners.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
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
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
double ** Pointers() const