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) {
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::
666 ApplyInverseGS_CrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 668 int NumVectors = X.NumVectors();
673 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
675 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
678 Y2 = Teuchos::rcp( &Y,
false );
680 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
681 X.ExtractView(&x_ptr);
682 Y.ExtractView(&y_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];
701 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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];
726 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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::
762 ApplyInverseGS_FastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 767 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
769 int NumVectors = X.NumVectors();
771 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
773 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
776 Y2 = Teuchos::rcp( &Y,
false );
778 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
779 X.ExtractView(&x_ptr);
780 Y.ExtractView(&y_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::
852 ApplyInverseGS_LocalFastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 857 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
859 int NumVectors = X.NumVectors();
861 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
863 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
866 Y2 = Teuchos::rcp( &Y,
false );
868 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
869 X.ExtractView(&x_ptr);
870 Y.ExtractView(&y_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::
943 ApplyInverseSGS(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 945 if (ZeroStartingSolution_)
948 const Epetra_CrsMatrix* CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*Matrix_);
953 if (CrsMatrix->StorageOptimized() && LocalSmoothingIndices_)
954 return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
955 else if (CrsMatrix->StorageOptimized())
956 return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
958 return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
961 return(ApplyInverseSGS_RowMatrix(X, Y));
965 int Ifpack_PointRelaxation::
966 ApplyInverseSGS_RowMatrix(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 968 int NumVectors = X.NumVectors();
969 int Length =
Matrix().MaxNumEntries();
970 std::vector<int> Indices(Length);
971 std::vector<double> Values(Length);
973 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
975 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
978 Y2 = Teuchos::rcp( &Y,
false );
980 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
981 X.ExtractView(&x_ptr);
982 Y.ExtractView(&y_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::
1055 ApplyInverseSGS_CrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1057 int NumVectors = X.NumVectors();
1062 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1064 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1067 Y2 = Teuchos::rcp( &Y,
false );
1069 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1070 X.ExtractView(&x_ptr);
1071 Y.ExtractView(&y_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];
1088 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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];
1110 IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
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::
1143 ApplyInverseSGS_FastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1148 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1150 int NumVectors = X.NumVectors();
1152 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1154 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1157 Y2 = Teuchos::rcp( &Y,
false );
1159 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1160 X.ExtractView(&x_ptr);
1161 Y.ExtractView(&y_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::
1231 ApplyInverseSGS_LocalFastCrsMatrix(
const Epetra_CrsMatrix* A,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1236 IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1238 int NumVectors = X.NumVectors();
1240 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1242 Y2 = Teuchos::rcp(
new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1245 Y2 = Teuchos::rcp( &Y,
false );
1247 double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1248 X.ExtractView(&x_ptr);
1249 Y.ExtractView(&y_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_);
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.