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" 55 #ifdef HAVE_IFPACK_AZTECOO 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),
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),
152 UseBlockMode_(
false),
153 SolveNormalEquations_(
false),
155 ZeroStartingSolution_(
true)
172 Epetra_Vector* ID = List.get(
"chebyshev: operator inv diagonal",
176 #ifdef HAVE_IFPACK_EPETRAEXT 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);
228 Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 233 if (X.NumVectors() != Y.NumVectors())
256 if (
Time_ == Teuchos::null)
257 Time_ = Teuchos::rcp(
new Epetra_Time(
Comm()) );
261 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
271 if (
Operator_->OperatorDomainMap().NumGlobalElements64() !=
272 Operator_->OperatorRangeMap().NumGlobalElements64())
288 Time_->ResetStartTime();
297 #ifdef HAVE_IFPACK_EPETRAEXT 300 const Epetra_CrsMatrix *CrsMatrix =
dynamic_cast<const Epetra_CrsMatrix*
>(&*
Matrix_);
309 if (InvBlockDiagonal_ == Teuchos::null)
312 ierr = InvBlockDiagonal_->SetParameters(BlockList_);
316 ierr = InvBlockDiagonal_->Compute();
322 double lambda_max = 0;
347 double diag = (*InvDiagonal_)[i];
364 #ifdef IFPACK_FLOPCOUNTERS 382 double MyMinVal, MyMaxVal;
383 double MinVal, MaxVal;
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;
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;
412 <<
" 0.0 0.0" << endl;
419 os <<
" " << std::setw(15) << 0.0 << endl;
426 os <<
" " << std::setw(15) << 0.0 << endl;
427 os <<
"================================================================================" << endl;
437 const int MaxIters,
const double Tol,
438 Epetra_RowMatrix* Matrix_in)
453 std::ostringstream oss;
454 oss <<
"\"Ifpack Chebyshev polynomial\": {" 456 <<
", Computed: " << (
IsComputed() ?
"true" :
"false")
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 497 IBD = &*InvBlockDiagonal_;
506 #ifdef HAVE_IFPACK_EPETRAEXT 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;
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;
563 #ifdef HAVE_IFPACK_EPETRAEXT 565 Temp.Update(oneOverTheta, X, -oneOverTheta, V, 0.0);
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 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;
640 double *xPointer = xPtr[0];
641 for (
int k = 0; k < degreeMinusOne; ++k) {
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);
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) {
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);
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);
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);
746 for (
int iter = 0; iter < MaximumIterations; ++iter)
748 Operator.Apply(x, y);
752 lambda_max = RQ_top / RQ_bottom;
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)
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);
823 for (
int iter = 0; iter < MaximumIterations; ++iter)
826 InvBlockDiagonal_->ApplyInverse(z,y);
828 InvBlockDiagonal_->ApplyInverse(y,z);
834 lambda_max = RQ_top / RQ_bottom;
845 #ifdef HAVE_IFPACK_EPETRAEXT 847 CG(
const int MaximumIterations,
848 double& lambda_min,
double& lambda_max,
const unsigned int * RngSeed)
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.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool UseBlockMode_
Use Block Preconditioning.
virtual int Compute()
Computes the preconditioners.
int NumMyNonzeros_
Number of local nonzeros.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int NumCompute_
Contains the number of successful call to Compute().
double EigRatio_
Contains the ratio such that [LambdaMax_ / EigRatio_, LambdaMax_] is the interval of interest for the...
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
std::string Label_
Contains the label of this object.
void Apply_Transpose(Teuchos::RCP< const Epetra_Operator > Operator_, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
int EigMaxIters_
Max number of iterations to use in eigenvalue estimation (if automatic).
long long NumGlobalNonzeros_
Number of global nonzeros.
long long NumGlobalRows_
Number of global rows.
double LambdaMin_
Contains an approximation to the smallest eigenvalue.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
double ComputeFlops_
Contains the number of flops for Compute().
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
int PolyDegree_
Contains the degree of Chebyshev polynomial.
bool IsRowMatrix_
If true, the Operator_ is an Epetra_RowMatrix.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
Teuchos::RefCountPtr< const Epetra_Operator > Operator_
Pointers to the matrix to be preconditioned as an Epetra_Operator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
double Condest_
Contains the estimated condition number.
bool SolveNormalEquations_
Run on the normal equations.
int NumMyRows_
Number of local rows.
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.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual void SetLabel()
Sets the label.
bool IsInitialized_
If true, the preconditioner has been computed successfully.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
double LambdaMax_
Contains an approximation to the largest eigenvalue.
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.
Teuchos::RefCountPtr< Epetra_Vector > InvDiagonal_
Contains the inverse of diagonal elements of Matrix.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
double ComputeTime_
Contains the time for all successful calls to Compute().
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.
#define IFPACK_CHK_ERR(ifpack_err)
bool IsComputed_
If true, the preconditioner has been computed successfully.
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.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned as an Epetra_RowMatrix.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double MinDiagonalValue_
Contains the minimum value on the diagonal.