43 #include "Ifpack_ConfigDefs.h" 45 #include "Epetra_Operator.h" 46 #include "Epetra_RowMatrix.h" 47 #include "Epetra_Comm.h" 48 #include "Epetra_Map.h" 49 #include "Epetra_MultiVector.h" 50 #include "Epetra_Vector.h" 51 #include "Epetra_Time.h" 52 #include "Ifpack_Chebyshev.h" 54 #include "Ifpack_Condest.h" 55 #ifdef HAVE_IFPACK_AZTECOO 56 #include "Ifpack_DiagPreconditioner.h" 60 #ifdef HAVE_IFPACK_EPETRAEXT 61 #include "Epetra_CrsMatrix.h" 62 #include "EpetraExt_PointToBlockDiagPermute.h" 66 #define ABS(x) ((x)>0?(x):-(x)) 69 inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,
const Epetra_MultiVector &X,Epetra_MultiVector &Y){
70 Epetra_Operator * Operator=
const_cast<Epetra_Operator*
>(&*Operator_);
71 Operator->SetUseTranspose(
true);
73 Operator->SetUseTranspose(
false);
84 IsInitialized_(false),
91 ApplyInverseTime_(0.0),
93 ApplyInverseFlops_(0.0),
102 MinDiagonalValue_(0.0),
106 NumGlobalNonzeros_(0),
107 Operator_(
Teuchos::rcp(Operator,false)),
108 UseBlockMode_(false),
109 SolveNormalEquations_(false),
111 ZeroStartingSolution_(true)
126 IsInitialized_(false),
131 InitializeTime_(0.0),
133 ApplyInverseTime_(0.0),
135 ApplyInverseFlops_(0.0),
137 UseTranspose_(false),
145 MinDiagonalValue_(0.0),
149 NumGlobalNonzeros_(0),
150 Operator_(
Teuchos::rcp(Operator,false)),
151 Matrix_(
Teuchos::rcp(Operator,false)),
152 UseBlockMode_(false),
153 SolveNormalEquations_(false),
155 ZeroStartingSolution_(true)
163 EigRatio_ = List.get(
"chebyshev: ratio eigenvalue", EigRatio_);
164 LambdaMin_ = List.get(
"chebyshev: min eigenvalue", LambdaMin_);
165 LambdaMax_ = List.get(
"chebyshev: max eigenvalue", LambdaMax_);
166 PolyDegree_ = List.get(
"chebyshev: degree",PolyDegree_);
167 MinDiagonalValue_ = List.get(
"chebyshev: min diagonal value",
169 ZeroStartingSolution_ = List.get(
"chebyshev: zero starting solution",
170 ZeroStartingSolution_);
172 Epetra_Vector* ID = List.get(
"chebyshev: operator inv diagonal",
174 EigMaxIters_ = List.get(
"chebyshev: eigenvalue max iterations",EigMaxIters_);
176 #ifdef HAVE_IFPACK_EPETRAEXT 178 UseBlockMode_ = List.get(
"chebyshev: use block mode",UseBlockMode_);
179 if(!List.isParameter(
"chebyshev: block list")) UseBlockMode_=
false;
181 BlockList_ = List.get(
"chebyshev: block list",BlockList_);
185 Teuchos::ParameterList Blist;
186 Blist=BlockList_.get(
"blockdiagmatrix: list",Blist);
187 std::string dummy(
"invert");
188 std::string ApplyMode=Blist.get(
"apply mode",dummy);
189 if(ApplyMode==std::string(
"multiply")){
190 Blist.set(
"apply mode",
"invert");
191 BlockList_.set(
"blockdiagmatrix: list",Blist);
196 SolveNormalEquations_ = List.get(
"chebyshev: solve normal equations",SolveNormalEquations_);
200 InvDiagonal_ = Teuchos::rcp(
new Epetra_Vector(*ID) );
211 return(Operator_->Comm());
217 return(Operator_->OperatorDomainMap());
223 return(Operator_->OperatorRangeMap());
228 Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 233 if (X.NumVectors() != Y.NumVectors())
242 IFPACK_CHK_ERR(Operator_->Apply(X,Y));
251 IsInitialized_ =
false;
253 if (Operator_ == Teuchos::null)
256 if (Time_ == Teuchos::null)
257 Time_ = Teuchos::rcp(
new Epetra_Time(
Comm()) );
261 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
264 NumMyRows_ = Matrix_->NumMyRows();
265 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
266 NumGlobalRows_ = Matrix_->NumGlobalRows64();
267 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
271 if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
272 Operator_->OperatorRangeMap().NumGlobalElements64())
277 InitializeTime_ += Time_->ElapsedTime();
278 IsInitialized_ =
true;
288 Time_->ResetStartTime();
294 if (PolyDegree_ <= 0)
297 #ifdef HAVE_IFPACK_EPETRAEXT 299 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
300 const Epetra_CrsMatrix *CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*Matrix_);
304 UseBlockMode_ =
false;
308 InvBlockDiagonal_ = Teuchos::rcp(
new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
309 if (InvBlockDiagonal_ == Teuchos::null)
312 ierr = InvBlockDiagonal_->SetParameters(BlockList_);
316 ierr = InvBlockDiagonal_->Compute();
322 double lambda_max = 0;
323 if (LambdaMax_ == -1) {
325 LambdaMax_ = lambda_max;
329 if (ABS(LambdaMax_-1) < 1e-6)
330 LambdaMax_ = LambdaMin_ = 1.0;
332 LambdaMin_ = LambdaMax_/EigRatio_;
336 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_) {
337 InvDiagonal_ = Teuchos::rcp(
new Epetra_Vector(
Matrix().Map()));
339 if (InvDiagonal_ == Teuchos::null)
342 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*InvDiagonal_));
346 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
347 double diag = (*InvDiagonal_)[i];
348 if (IFPACK_ABS(diag) < MinDiagonalValue_)
349 (*InvDiagonal_)[i] = MinDiagonalValue_;
351 (*InvDiagonal_)[i] = 1.0 / diag;
355 if (LambdaMax_ == -1) {
357 LambdaMax_=lambda_max;
359 if (ABS(LambdaMax_-1) < 1e-6) LambdaMax_=LambdaMin_=1.0;
360 else LambdaMin_=LambdaMax_/EigRatio_;
364 #ifdef IFPACK_FLOPCOUNTERS 365 ComputeFlops_ += NumMyRows_;
369 ComputeTime_ += Time_->ElapsedTime();
382 double MyMinVal, MyMaxVal;
383 double MinVal, MaxVal;
386 InvDiagonal_->MinValue(&MyMinVal);
387 InvDiagonal_->MaxValue(&MyMaxVal);
388 Comm().MinAll(&MyMinVal,&MinVal,1);
389 Comm().MinAll(&MyMaxVal,&MaxVal,1);
392 if (!
Comm().MyPID()) {
394 os <<
"================================================================================" << endl;
395 os <<
"Ifpack_Chebyshev" << endl;
396 os <<
"Degree of polynomial = " << PolyDegree_ << endl;
397 os <<
"Condition number estimate = " <<
Condest() << endl;
398 os <<
"Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
400 os <<
"Minimum value on stored inverse diagonal = " << MinVal << endl;
401 os <<
"Maximum value on stored inverse diagonal = " << MaxVal << endl;
403 if (ZeroStartingSolution_)
404 os <<
"Using zero starting solution" << endl;
406 os <<
"Using input starting solution" << endl;
408 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
409 os <<
"----- ------- -------------- ------------ --------" << endl;
410 os <<
"Initialize() " << std::setw(5) << NumInitialize_
411 <<
" " << std::setw(15) << InitializeTime_
412 <<
" 0.0 0.0" << endl;
413 os <<
"Compute() " << std::setw(5) << NumCompute_
414 <<
" " << std::setw(15) << ComputeTime_
415 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
416 if (ComputeTime_ != 0.0)
417 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
419 os <<
" " << std::setw(15) << 0.0 << endl;
420 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
421 <<
" " << std::setw(15) << ApplyInverseTime_
422 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
423 if (ApplyInverseTime_ != 0.0)
424 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
426 os <<
" " << std::setw(15) << 0.0 << endl;
427 os <<
"================================================================================" << endl;
437 const int MaxIters,
const double Tol,
438 Epetra_RowMatrix* Matrix_in)
445 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
451 void Ifpack_Chebyshev::SetLabel()
453 std::ostringstream oss;
454 oss <<
"\"Ifpack Chebyshev polynomial\": {" 456 <<
", Computed: " << (
IsComputed() ?
"true" :
"false")
457 <<
", degree: " << PolyDegree_
458 <<
", lambdaMax: " << LambdaMax_
459 <<
", alpha: " << EigRatio_
460 <<
", lambdaMin: " << LambdaMin_
472 if (PolyDegree_ == 0)
475 int nVec = X.NumVectors();
476 int len = X.MyLength();
477 if (nVec != Y.NumVectors())
480 Time_->ResetStartTime();
484 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
485 if (X.Pointers()[0] == Y.Pointers()[0])
486 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
488 Xcopy = Teuchos::rcp( &X,
false );
490 double **xPtr = 0, **yPtr = 0;
491 Xcopy->ExtractView(&xPtr);
492 Y.ExtractView(&yPtr);
494 #ifdef HAVE_IFPACK_EPETRAEXT 495 EpetraExt_PointToBlockDiagPermute* IBD=0;
497 IBD = &*InvBlockDiagonal_;
503 invDiag = InvDiagonal_->Values();
505 if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
506 #ifdef HAVE_IFPACK_EPETRAEXT 508 IBD->ApplyInverse(*Xcopy, Y);
512 double *yPointer = yPtr[0], *xPointer = xPtr[0];
513 for (
int i = 0; i < len; ++i)
514 yPointer[i] = xPointer[i] * invDiag[i];
517 for (
int i = 0; i < len; ++i) {
518 double coeff = invDiag[i];
519 for (
int k = 0; k < nVec; ++k)
520 yPtr[k][i] = xPtr[k][i] * coeff;
529 double alpha = LambdaMax_ / EigRatio_;
530 double beta = 1.1 * LambdaMax_;
531 double delta = 2.0 / (beta - alpha);
532 double theta = 0.5 * (beta + alpha);
533 double s1 = theta * delta;
540 bool zeroOut =
false;
541 Epetra_MultiVector V(X.Map(), X.NumVectors(), zeroOut);
542 Epetra_MultiVector W(X.Map(), X.NumVectors(), zeroOut);
543 #ifdef HAVE_IFPACK_EPETRAEXT 544 Epetra_MultiVector Temp(X.Map(), X.NumVectors(), zeroOut);
547 double *vPointer = V.Values(), *wPointer = W.Values();
549 double oneOverTheta = 1.0/theta;
552 if (SolveNormalEquations_) {
553 Apply_Transpose(Operator_, Y, V);
559 if (ZeroStartingSolution_ ==
false) {
560 Operator_->Apply(Y, V);
563 #ifdef HAVE_IFPACK_EPETRAEXT 565 Temp.Update(oneOverTheta, X, -oneOverTheta, V, 0.0);
566 IBD->ApplyInverse(Temp, W);
570 if (SolveNormalEquations_){
571 IBD->ApplyInverse(W, Temp);
572 Apply_Transpose(Operator_, Temp, W);
578 double *xPointer = xPtr[0];
579 for (
int i = 0; i < len; ++i)
580 wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
583 for (
int i = 0; i < len; ++i) {
584 double coeff = invDiag[i]*oneOverTheta;
585 double *wi = wPointer + i, *vi = vPointer + i;
586 for (
int k = 0; k < nVec; ++k) {
587 *wi = (xPtr[k][i] - (*vi)) * coeff;
588 wi = wi + len; vi = vi + len;
593 Y.Update(1.0, W, 1.0);
597 #ifdef HAVE_IFPACK_EPETRAEXT 599 IBD->ApplyInverse(X, W);
603 if (SolveNormalEquations_) {
604 IBD->ApplyInverse(W, Temp);
605 Apply_Transpose(Operator_, Temp, W);
608 W.Scale(oneOverTheta);
609 Y.Update(1.0, W, 0.0);
614 double *xPointer = xPtr[0];
615 for (
int i = 0; i < len; ++i)
616 wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
618 memcpy(yPtr[0], wPointer, len*
sizeof(
double));
621 for (
int i = 0; i < len; ++i) {
622 double coeff = invDiag[i]*oneOverTheta;
623 double *wi = wPointer + i;
624 for (
int k = 0; k < nVec; ++k) {
625 *wi = xPtr[k][i] * coeff;
630 for (
int k = 0; k < nVec; ++k)
631 memcpy(yPtr[k], wPointer + k*len, len*
sizeof(
double));
636 double rhok = 1.0/s1, rhokp1;
637 double dtemp1, dtemp2;
638 int degreeMinusOne = PolyDegree_ - 1;
640 double *xPointer = xPtr[0];
641 for (
int k = 0; k < degreeMinusOne; ++k) {
642 Operator_->Apply(Y, V);
643 rhokp1 = 1.0 / (2.0*s1 - rhok);
644 dtemp1 = rhokp1 * rhok;
645 dtemp2 = 2.0 * rhokp1 * delta;
652 #ifdef HAVE_IFPACK_EPETRAEXT 655 V.Update(dtemp2, X, -dtemp2);
656 IBD->ApplyInverse(V, Temp);
660 if (SolveNormalEquations_) {
661 IBD->ApplyInverse(V, Temp);
662 Apply_Transpose(Operator_, Temp, V);
665 W.Update(1.0, Temp, 1.0);
669 for (
int i = 0; i < len; ++i)
670 wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
673 Y.Update(1.0, W, 1.0);
677 for (
int k = 0; k < degreeMinusOne; ++k) {
678 Operator_->Apply(Y, V);
679 rhokp1 = 1.0 / (2.0*s1 - rhok);
680 dtemp1 = rhokp1 * rhok;
681 dtemp2 = 2.0 * rhokp1 * delta;
688 #ifdef HAVE_IFPACK_EPETRAEXT 691 V.Update(dtemp2, X, -dtemp2);
692 IBD->ApplyInverse(V, Temp);
696 if (SolveNormalEquations_) {
697 IBD->ApplyInverse(V,Temp);
698 Apply_Transpose(Operator_,Temp,V);
701 W.Update(1.0, Temp, 1.0);
705 for (
int i = 0; i < len; ++i) {
706 double coeff = invDiag[i]*dtemp2;
707 double *wi = wPointer + i, *vi = vPointer + i;
708 for (
int j = 0; j < nVec; ++j) {
709 *wi += (xPtr[j][i] - (*vi)) * coeff;
710 wi = wi + len; vi = vi + len;
715 Y.Update(1.0, W, 1.0);
722 ApplyInverseTime_ += Time_->ElapsedTime();
730 const Epetra_Vector& InvPointDiagonal,
731 const int MaximumIterations,
732 double& lambda_max,
const unsigned int * RngSeed)
736 double RQ_top, RQ_bottom, norm;
737 Epetra_Vector x(Operator.OperatorDomainMap());
738 Epetra_Vector y(Operator.OperatorRangeMap());
739 if(RngSeed) x.SetSeed(*RngSeed);
742 if (norm == 0.0) IFPACK_CHK_ERR(-1);
746 for (
int iter = 0; iter < MaximumIterations; ++iter)
748 Operator.Apply(x, y);
749 IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
750 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
751 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
752 lambda_max = RQ_top / RQ_bottom;
753 IFPACK_CHK_ERR(y.Norm2(&norm));
754 if (norm == 0.0) IFPACK_CHK_ERR(-1);
755 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
763 CG(
const Epetra_Operator& Operator,
764 const Epetra_Vector& InvPointDiagonal,
765 const int MaximumIterations,
766 double& lambda_min,
double& lambda_max,
const unsigned int * RngSeed)
768 #ifdef HAVE_IFPACK_AZTECOO 769 Epetra_Vector x(Operator.OperatorDomainMap());
770 Epetra_Vector y(Operator.OperatorRangeMap());
771 if(RngSeed) x.SetSeed(*RngSeed);
775 Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
777 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
778 solver.SetAztecOption(AZ_output, AZ_none);
781 Operator.OperatorRangeMap(),
783 solver.SetPrecOperator(&diag);
784 solver.Iterate(MaximumIterations, 1e-10);
786 const double* status = solver.GetAztecStatus();
788 lambda_min = status[AZ_lambda_min];
789 lambda_max = status[AZ_lambda_max];
796 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
797 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
798 cout <<
"in your configure script." << endl;
804 #ifdef HAVE_IFPACK_EPETRAEXT 806 PowerMethod(
const int MaximumIterations,
double& lambda_max,
const unsigned int * RngSeed)
809 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
812 double RQ_top, RQ_bottom, norm;
813 Epetra_Vector x(Operator_->OperatorDomainMap());
814 Epetra_Vector y(Operator_->OperatorRangeMap());
815 Epetra_Vector z(Operator_->OperatorRangeMap());
816 if(RngSeed) x.SetSeed(*RngSeed);
819 if (norm == 0.0) IFPACK_CHK_ERR(-1);
823 for (
int iter = 0; iter < MaximumIterations; ++iter)
825 Operator_->Apply(x, z);
826 InvBlockDiagonal_->ApplyInverse(z,y);
827 if(SolveNormalEquations_){
828 InvBlockDiagonal_->ApplyInverse(y,z);
829 Apply_Transpose(Operator_,z, y);
832 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
833 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
834 lambda_max = RQ_top / RQ_bottom;
835 IFPACK_CHK_ERR(y.Norm2(&norm));
836 if (norm == 0.0) IFPACK_CHK_ERR(-1);
837 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
845 #ifdef HAVE_IFPACK_EPETRAEXT 847 CG(
const int MaximumIterations,
848 double& lambda_min,
double& lambda_max,
const unsigned int * RngSeed)
856 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
858 #ifdef HAVE_IFPACK_AZTECOO 859 Epetra_Vector x(Operator_->OperatorDomainMap());
860 Epetra_Vector y(Operator_->OperatorRangeMap());
861 if(RngSeed) x.SetSeed(*RngSeed);
864 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
867 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
868 solver.SetAztecOption(AZ_output, AZ_none);
870 solver.SetPrecOperator(&*InvBlockDiagonal_);
871 solver.Iterate(MaximumIterations, 1e-10);
873 const double* status = solver.GetAztecStatus();
875 lambda_min = status[AZ_lambda_min];
876 lambda_max = status[AZ_lambda_max];
883 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
884 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
885 cout <<
"in your configure script." << endl;
891 #endif // HAVE_IFPACK_EPETRAEXT virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual int Compute()
Computes the preconditioners.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
static int CG(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_min, double &lambda_max, const unsigned int *RngSeed=0)
Uses AztecOO's CG to estimate lambda_min and lambda_max.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
static int PowerMethod(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &LambdaMax, const unsigned int *RngSeed=0)
Simple power method to compute lambda_max.