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))
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_);
211 return(Operator_->Comm());
217 return(Operator_->OperatorDomainMap());
223 return(Operator_->OperatorRangeMap());
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)
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_){
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_) {
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);
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,
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])
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;
543 #ifdef HAVE_IFPACK_EPETRAEXT
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();
731 const int MaximumIterations,
732 double& lambda_max,
const unsigned int * RngSeed)
736 double RQ_top, RQ_bottom, norm;
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));
765 const int MaximumIterations,
766 double& lambda_min,
double& lambda_max,
const unsigned int * RngSeed)
768 #ifdef HAVE_IFPACK_AZTECOO
771 if(RngSeed) x.SetSeed(*RngSeed);
777 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
778 solver.SetAztecOption(AZ_output, AZ_none);
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;
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
848 double& ,
double& ,
const unsigned int * )
856 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
858 #ifdef HAVE_IFPACK_AZTECOO
861 if(RngSeed) x.SetSeed(*RngSeed);
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 int Compute()
Computes the preconditioners.
virtual int SetUseTranspose(bool UseTranspose)=0
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
virtual const Epetra_Map & OperatorDomainMap() const =0
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
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 UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
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.