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)
179 Time_ = Teuchos::rcp(
new Epetra_Time(
Comm()) );
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 const Epetra_VbrMatrix * VbrMat =
dynamic_cast<const Epetra_VbrMatrix*
>(&
Matrix());
219 if(VbrMat) Diagonal_ = Teuchos::rcp(
new Epetra_Vector(VbrMat->RowMap()) );
220 else Diagonal_ = Teuchos::rcp(
new Epetra_Vector(
Matrix().RowMatrixRowMap()) );
222 if (Diagonal_ == Teuchos::null)
225 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*Diagonal_));
236 if(DoL1Method_ && IsParallel_) {
237 int maxLength =
Matrix().MaxNumEntries();
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))) {
285 Importer_ = Teuchos::rcp(
new Epetra_Import(
Matrix().RowMatrixColMap(),
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);
309 Comm().MinAll(&MyMinVal,&MinVal,1);
310 Comm().MinAll(&MyMaxVal,&MaxVal,1);
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,
368 Epetra_RowMatrix* Matrix_in)
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;
416 if (X.NumVectors() != Y.NumVectors())
419 Time_->ResetStartTime();
423 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
424 if (X.Pointers()[0] == Y.Pointers()[0])
425 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
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::
454 ApplyInverseJacobi(
const Epetra_MultiVector& RHS, Epetra_MultiVector& LHS)
const 456 int NumVectors = LHS.NumVectors();
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;
469 Epetra_MultiVector A_times_LHS(LHS.Map(), NumVectors, zeroOut);
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::
496 ApplyInverseGS(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 498 if (ZeroStartingSolution_)
501 const Epetra_CrsMatrix* CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*Matrix_);
506 if (CrsMatrix->StorageOptimized() && LocalSmoothingIndices_)
507 return(ApplyInverseGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
508 else if (CrsMatrix->StorageOptimized())
509 return(ApplyInverseGS_FastCrsMatrix(CrsMatrix, X, Y));
511 return(ApplyInverseGS_CrsMatrix(CrsMatrix, X, Y));
514 return(ApplyInverseGS_RowMatrix(X, Y));
520 int Ifpack_PointRelaxation::
521 ApplyInverseGS_RowMatrix(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 523 int NumVectors = X.NumVectors();
525 int Length =
Matrix().MaxNumEntries();
526 std::vector<int> Indices(Length);
527 std::vector<double> Values(Length);
529 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
531 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
533 Y2 = Teuchos::rcp( &Y,
false );
536 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
537 X.ExtractView(&x_ptr);
538 Y.ExtractView(&y_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::
667 ApplyInverseGS_CrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 669 int NumVectors = X.NumVectors();
674 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
676 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
679 Y2 = Teuchos::rcp( &Y,
false );
681 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
682 X.ExtractView(&x_ptr);
683 Y.ExtractView(&y_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];
702 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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];
727 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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::
763 ApplyInverseGS_FastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 768 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
770 int NumVectors = X.NumVectors();
772 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
774 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
777 Y2 = Teuchos::rcp( &Y,
false );
779 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
780 X.ExtractView(&x_ptr);
781 Y.ExtractView(&y_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::
853 ApplyInverseGS_LocalFastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 858 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
860 int NumVectors = X.NumVectors();
862 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
864 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
867 Y2 = Teuchos::rcp( &Y,
false );
869 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
870 X.ExtractView(&x_ptr);
871 Y.ExtractView(&y_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::
944 ApplyInverseSGS(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 946 if (ZeroStartingSolution_)
949 const Epetra_CrsMatrix* CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*Matrix_);
954 if (CrsMatrix->StorageOptimized() && LocalSmoothingIndices_)
955 return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
956 else if (CrsMatrix->StorageOptimized())
957 return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
959 return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
962 return(ApplyInverseSGS_RowMatrix(X, Y));
966 int Ifpack_PointRelaxation::
967 ApplyInverseSGS_RowMatrix(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 969 int NumVectors = X.NumVectors();
970 int Length =
Matrix().MaxNumEntries();
971 std::vector<int> Indices(Length);
972 std::vector<double> Values(Length);
974 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
976 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
979 Y2 = Teuchos::rcp( &Y,
false );
981 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
982 X.ExtractView(&x_ptr);
983 Y.ExtractView(&y_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::
1056 ApplyInverseSGS_CrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1058 int NumVectors = X.NumVectors();
1063 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1065 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1068 Y2 = Teuchos::rcp( &Y,
false );
1070 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1071 X.ExtractView(&x_ptr);
1072 Y.ExtractView(&y_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];
1089 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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];
1111 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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::
1144 ApplyInverseSGS_FastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1149 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1151 int NumVectors = X.NumVectors();
1153 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1155 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1158 Y2 = Teuchos::rcp( &Y,
false );
1160 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1161 X.ExtractView(&x_ptr);
1162 Y.ExtractView(&y_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::
1232 ApplyInverseSGS_LocalFastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1237 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1239 int NumVectors = X.NumVectors();
1241 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1243 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1246 Y2 = Teuchos::rcp( &Y,
false );
1248 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1249 X.ExtractView(&x_ptr);
1250 Y.ExtractView(&y_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_);
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
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 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 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.
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.
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.