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) {
588 dtemp += Values[k] * y20_ptr[i];
591 y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
597 for (
int i = 0 ; i < NumMyRows_ ; ++i)
598 y0_ptr[i] = y20_ptr[i];
604 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
608 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
609 &Values[0], &Indices[0]));
611 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
626 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
629 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
630 &Values[0], &Indices[0]));
632 for (
int m = 0 ; m < NumVectors ; ++m) {
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);
649 for (
int m = 0 ; m < NumVectors ; ++m)
650 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
651 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
652 y_ptr[m][i] = y2_ptr[m][i];
657 #ifdef IFPACK_FLOPCOUNTERS 658 ApplyInverseFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
665 int Ifpack_PointRelaxation::
673 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
678 Y2 = Teuchos::rcp( &Y,
false );
680 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
683 Y2->ExtractView(&y2_ptr);
684 Diagonal_->ExtractView(&d_ptr);
686 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
690 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
694 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
695 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
699 double diag = d_ptr[i];
703 for (
int m = 0 ; m < NumVectors ; ++m) {
707 for (
int k = 0; k < NumEntries; ++k) {
710 dtemp += Values[k] * y2_ptr[m][col];
713 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
719 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
720 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
724 double diag = d_ptr[i];
728 for (
int m = 0 ; m < NumVectors ; ++m) {
731 for (
int k = 0; k < NumEntries; ++k) {
734 dtemp += Values[k] * y2_ptr[m][col];
737 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
744 for (
int m = 0 ; m < NumVectors ; ++m)
745 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
746 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
747 y_ptr[m][i] = y2_ptr[m][i];
751 #ifdef IFPACK_FLOPCOUNTERS 752 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
761 int Ifpack_PointRelaxation::
771 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
776 Y2 = Teuchos::rcp( &Y,
false );
778 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
781 Y2->ExtractView(&y2_ptr);
782 Diagonal_->ExtractView(&d_ptr);
784 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
788 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
792 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
795 double diag = d_ptr[i];
797 for (
int m = 0 ; m < NumVectors ; ++m) {
801 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
804 dtemp += Values[k] * y2_ptr[m][col];
807 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
813 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
816 double diag = d_ptr[i];
818 for (
int m = 0 ; m < NumVectors ; ++m) {
821 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
824 dtemp += Values[k] * y2_ptr[m][col];
827 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
835 for (
int m = 0 ; m < NumVectors ; ++m)
836 for (
int i = 0 ; i < NumMyRows_ ; ++i)
837 y_ptr[m][i] = y2_ptr[m][i];
840 #ifdef IFPACK_FLOPCOUNTERS 841 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
851 int Ifpack_PointRelaxation::
861 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
866 Y2 = Teuchos::rcp( &Y,
false );
868 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
871 Y2->ExtractView(&y2_ptr);
872 Diagonal_->ExtractView(&d_ptr);
874 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
878 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
882 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
883 int i=LocalSmoothingIndices_[ii];
886 double diag = d_ptr[i];
888 for (
int m = 0 ; m < NumVectors ; ++m) {
892 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
895 dtemp += Values[k] * y2_ptr[m][col];
898 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
904 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
905 int i=LocalSmoothingIndices_[ii];
908 double diag = d_ptr[i];
910 for (
int m = 0 ; m < NumVectors ; ++m) {
913 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
916 dtemp += Values[k] * y2_ptr[m][col];
919 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
927 for (
int m = 0 ; m < NumVectors ; ++m)
928 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
929 int i=LocalSmoothingIndices_[ii];
930 y_ptr[m][i] = y2_ptr[m][i];
934 #ifdef IFPACK_FLOPCOUNTERS 935 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
942 int Ifpack_PointRelaxation::
945 if (ZeroStartingSolution_)
954 return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
956 return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
958 return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
961 return(ApplyInverseSGS_RowMatrix(X, Y));
965 int Ifpack_PointRelaxation::
970 std::vector<int> Indices(Length);
971 std::vector<double> Values(Length);
973 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
978 Y2 = Teuchos::rcp( &Y,
false );
980 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
983 Y2->ExtractView(&y2_ptr);
984 Diagonal_->ExtractView(&d_ptr);
986 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
990 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
992 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
993 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
997 double diag = d_ptr[i];
999 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1000 &Values[0], &Indices[0]));
1002 for (
int m = 0 ; m < NumVectors ; ++m) {
1006 for (
int k = 0 ; k < NumEntries ; ++k) {
1009 dtemp += Values[k] * y2_ptr[m][col];
1012 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1015 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1016 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1020 double diag = d_ptr[i];
1022 IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1023 &Values[0], &Indices[0]));
1025 for (
int m = 0 ; m < NumVectors ; ++m) {
1028 for (
int k = 0 ; k < NumEntries ; ++k) {
1031 dtemp += Values[k] * y2_ptr[m][col];
1034 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1040 for (
int m = 0 ; m < NumVectors ; ++m)
1041 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1042 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1043 y_ptr[m][i] = y2_ptr[m][i];
1047 #ifdef IFPACK_FLOPCOUNTERS 1048 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1054 int Ifpack_PointRelaxation::
1062 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1067 Y2 = Teuchos::rcp( &Y,
false );
1069 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1072 Y2->ExtractView(&y2_ptr);
1073 Diagonal_->ExtractView(&d_ptr);
1075 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1079 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1081 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1082 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1086 double diag = d_ptr[i];
1090 for (
int m = 0 ; m < NumVectors ; ++m) {
1094 for (
int k = 0; k < NumEntries; ++k) {
1097 dtemp += Values[k] * y2_ptr[m][col];
1100 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1103 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1104 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1108 double diag = d_ptr[i];
1112 for (
int m = 0 ; m < NumVectors ; ++m) {
1115 for (
int k = 0; k < NumEntries; ++k) {
1118 dtemp += Values[k] * y2_ptr[m][col];
1121 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1127 for (
int m = 0 ; m < NumVectors ; ++m)
1128 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1129 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1130 y_ptr[m][i] = y2_ptr[m][i];
1134 #ifdef IFPACK_FLOPCOUNTERS 1135 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1142 int Ifpack_PointRelaxation::
1152 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1157 Y2 = Teuchos::rcp( &Y,
false );
1159 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1162 Y2->ExtractView(&y2_ptr);
1163 Diagonal_->ExtractView(&d_ptr);
1165 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1169 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1171 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
1174 double diag = d_ptr[i];
1179 for (
int m = 0 ; m < NumVectors ; ++m) {
1183 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1186 dtemp += Values[k] * y2_ptr[m][col];
1189 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1193 for (
int i = NumMyRows_ - 1 ; i > -1 ; --i) {
1196 double diag = d_ptr[i];
1201 for (
int m = 0 ; m < NumVectors ; ++m) {
1204 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1207 dtemp += Values[k] * y2_ptr[m][col];
1210 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1216 for (
int m = 0 ; m < NumVectors ; ++m)
1217 for (
int i = 0 ; i < NumMyRows_ ; ++i)
1218 y_ptr[m][i] = y2_ptr[m][i];
1221 #ifdef IFPACK_FLOPCOUNTERS 1222 ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1230 int Ifpack_PointRelaxation::
1240 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1245 Y2 = Teuchos::rcp( &Y,
false );
1247 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1250 Y2->ExtractView(&y2_ptr);
1251 Diagonal_->ExtractView(&d_ptr);
1253 for (
int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1257 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1259 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1260 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1263 double diag = d_ptr[i];
1268 for (
int m = 0 ; m < NumVectors ; ++m) {
1272 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1275 dtemp += Values[k] * y2_ptr[m][col];
1278 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1282 for (
int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1283 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1286 double diag = d_ptr[i];
1291 for (
int m = 0 ; m < NumVectors ; ++m) {
1294 for (
int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1297 dtemp += Values[k] * y2_ptr[m][col];
1300 y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1306 for (
int m = 0 ; m < NumVectors ; ++m)
1307 for (
int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1308 int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1309 y_ptr[m][i] = y2_ptr[m][i];
1313 #ifdef IFPACK_FLOPCOUNTERS 1314 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