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;
324 LambdaMax_ = lambda_max;
327 if (ABS(LambdaMax_-1) < 1e-6)
328 LambdaMax_ = LambdaMin_ = 1.0;
330 LambdaMin_ = LambdaMax_/EigRatio_;
334 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_) {
335 InvDiagonal_ = Teuchos::rcp(
new Epetra_Vector(
Matrix().Map()));
337 if (InvDiagonal_ == Teuchos::null)
340 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*InvDiagonal_));
344 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
345 double diag = (*InvDiagonal_)[i];
346 if (IFPACK_ABS(diag) < MinDiagonalValue_)
347 (*InvDiagonal_)[i] = MinDiagonalValue_;
349 (*InvDiagonal_)[i] = 1.0 / diag;
353 if (LambdaMax_ == -1) {
355 LambdaMax_=lambda_max;
357 if (ABS(LambdaMax_-1) < 1e-6) LambdaMax_=LambdaMin_=1.0;
358 else LambdaMin_=LambdaMax_/EigRatio_;
362 #ifdef IFPACK_FLOPCOUNTERS 363 ComputeFlops_ += NumMyRows_;
367 ComputeTime_ += Time_->ElapsedTime();
380 double MyMinVal, MyMaxVal;
381 double MinVal, MaxVal;
384 InvDiagonal_->MinValue(&MyMinVal);
385 InvDiagonal_->MaxValue(&MyMaxVal);
386 Comm().MinAll(&MyMinVal,&MinVal,1);
387 Comm().MinAll(&MyMaxVal,&MaxVal,1);
390 if (!
Comm().MyPID()) {
392 os <<
"================================================================================" << endl;
393 os <<
"Ifpack_Chebyshev" << endl;
394 os <<
"Degree of polynomial = " << PolyDegree_ << endl;
395 os <<
"Condition number estimate = " <<
Condest() << endl;
396 os <<
"Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
398 os <<
"Minimum value on stored inverse diagonal = " << MinVal << endl;
399 os <<
"Maximum value on stored inverse diagonal = " << MaxVal << endl;
401 if (ZeroStartingSolution_)
402 os <<
"Using zero starting solution" << endl;
404 os <<
"Using input starting solution" << endl;
406 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
407 os <<
"----- ------- -------------- ------------ --------" << endl;
408 os <<
"Initialize() " << std::setw(5) << NumInitialize_
409 <<
" " << std::setw(15) << InitializeTime_
410 <<
" 0.0 0.0" << endl;
411 os <<
"Compute() " << std::setw(5) << NumCompute_
412 <<
" " << std::setw(15) << ComputeTime_
413 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
414 if (ComputeTime_ != 0.0)
415 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
417 os <<
" " << std::setw(15) << 0.0 << endl;
418 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
419 <<
" " << std::setw(15) << ApplyInverseTime_
420 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
421 if (ApplyInverseTime_ != 0.0)
422 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
424 os <<
" " << std::setw(15) << 0.0 << endl;
425 os <<
"================================================================================" << endl;
435 const int MaxIters,
const double Tol,
436 Epetra_RowMatrix* Matrix_in)
443 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
449 void Ifpack_Chebyshev::SetLabel()
451 std::ostringstream oss;
452 oss <<
"\"Ifpack Chebyshev polynomial\": {" 454 <<
", Computed: " << (
IsComputed() ?
"true" :
"false")
455 <<
", degree: " << PolyDegree_
456 <<
", lambdaMax: " << LambdaMax_
457 <<
", alpha: " << EigRatio_
458 <<
", lambdaMin: " << LambdaMin_
470 if (PolyDegree_ == 0)
473 int nVec = X.NumVectors();
474 int len = X.MyLength();
475 if (nVec != Y.NumVectors())
478 Time_->ResetStartTime();
482 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
483 if (X.Pointers()[0] == Y.Pointers()[0])
484 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
486 Xcopy = Teuchos::rcp( &X,
false );
488 double **xPtr = 0, **yPtr = 0;
489 Xcopy->ExtractView(&xPtr);
490 Y.ExtractView(&yPtr);
492 #ifdef HAVE_IFPACK_EPETRAEXT 493 EpetraExt_PointToBlockDiagPermute* IBD=0;
495 IBD = &*InvBlockDiagonal_;
501 invDiag = InvDiagonal_->Values();
503 if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
504 #ifdef HAVE_IFPACK_EPETRAEXT 506 IBD->ApplyInverse(*Xcopy, Y);
510 double *yPointer = yPtr[0], *xPointer = xPtr[0];
511 for (
int i = 0; i < len; ++i)
512 yPointer[i] = xPointer[i] * invDiag[i];
515 for (
int i = 0; i < len; ++i) {
516 double coeff = invDiag[i];
517 for (
int k = 0; k < nVec; ++k)
518 yPtr[k][i] = xPtr[k][i] * coeff;
527 double alpha = LambdaMax_ / EigRatio_;
528 double beta = 1.1 * LambdaMax_;
529 double delta = 2.0 / (beta - alpha);
530 double theta = 0.5 * (beta + alpha);
531 double s1 = theta * delta;
538 bool zeroOut =
false;
539 Epetra_MultiVector V(X.Map(), X.NumVectors(), zeroOut);
540 Epetra_MultiVector W(X.Map(), X.NumVectors(), zeroOut);
541 #ifdef HAVE_IFPACK_EPETRAEXT 542 Epetra_MultiVector Temp(X.Map(), X.NumVectors(), zeroOut);
545 double *vPointer = V.Values(), *wPointer = W.Values();
547 double oneOverTheta = 1.0/theta;
550 if (SolveNormalEquations_) {
551 Apply_Transpose(Operator_, Y, V);
557 if (ZeroStartingSolution_ ==
false) {
558 Operator_->Apply(Y, V);
561 #ifdef HAVE_IFPACK_EPETRAEXT 563 Temp.Update(oneOverTheta, X, -oneOverTheta, V, 0.0);
564 IBD->ApplyInverse(Temp, W);
568 if (SolveNormalEquations_){
569 IBD->ApplyInverse(W, Temp);
570 Apply_Transpose(Operator_, Temp, W);
576 double *xPointer = xPtr[0];
577 for (
int i = 0; i < len; ++i)
578 wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
581 for (
int i = 0; i < len; ++i) {
582 double coeff = invDiag[i]*oneOverTheta;
583 double *wi = wPointer + i, *vi = vPointer + i;
584 for (
int k = 0; k < nVec; ++k) {
585 *wi = (xPtr[k][i] - (*vi)) * coeff;
586 wi = wi + len; vi = vi + len;
591 Y.Update(1.0, W, 1.0);
595 #ifdef HAVE_IFPACK_EPETRAEXT 597 IBD->ApplyInverse(X, W);
601 if (SolveNormalEquations_) {
602 IBD->ApplyInverse(W, Temp);
603 Apply_Transpose(Operator_, Temp, W);
606 W.Scale(oneOverTheta);
607 Y.Update(1.0, W, 0.0);
612 double *xPointer = xPtr[0];
613 for (
int i = 0; i < len; ++i)
614 wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
616 memcpy(yPtr[0], wPointer, len*
sizeof(
double));
619 for (
int i = 0; i < len; ++i) {
620 double coeff = invDiag[i]*oneOverTheta;
621 double *wi = wPointer + i;
622 for (
int k = 0; k < nVec; ++k) {
623 *wi = xPtr[k][i] * coeff;
628 for (
int k = 0; k < nVec; ++k)
629 memcpy(yPtr[k], wPointer + k*len, len*
sizeof(
double));
634 double rhok = 1.0/s1, rhokp1;
635 double dtemp1, dtemp2;
636 int degreeMinusOne = PolyDegree_ - 1;
638 double *xPointer = xPtr[0];
639 for (
int k = 0; k < degreeMinusOne; ++k) {
640 Operator_->Apply(Y, V);
641 rhokp1 = 1.0 / (2.0*s1 - rhok);
642 dtemp1 = rhokp1 * rhok;
643 dtemp2 = 2.0 * rhokp1 * delta;
650 #ifdef HAVE_IFPACK_EPETRAEXT 653 V.Update(dtemp2, X, -dtemp2);
654 IBD->ApplyInverse(V, Temp);
658 if (SolveNormalEquations_) {
659 IBD->ApplyInverse(V, Temp);
660 Apply_Transpose(Operator_, Temp, V);
663 W.Update(1.0, Temp, 1.0);
667 for (
int i = 0; i < len; ++i)
668 wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
671 Y.Update(1.0, W, 1.0);
675 for (
int k = 0; k < degreeMinusOne; ++k) {
676 Operator_->Apply(Y, V);
677 rhokp1 = 1.0 / (2.0*s1 - rhok);
678 dtemp1 = rhokp1 * rhok;
679 dtemp2 = 2.0 * rhokp1 * delta;
686 #ifdef HAVE_IFPACK_EPETRAEXT 689 V.Update(dtemp2, X, -dtemp2);
690 IBD->ApplyInverse(V, Temp);
694 if (SolveNormalEquations_) {
695 IBD->ApplyInverse(V,Temp);
696 Apply_Transpose(Operator_,Temp,V);
699 W.Update(1.0, Temp, 1.0);
703 for (
int i = 0; i < len; ++i) {
704 double coeff = invDiag[i]*dtemp2;
705 double *wi = wPointer + i, *vi = vPointer + i;
706 for (
int j = 0; j < nVec; ++j) {
707 *wi += (xPtr[j][i] - (*vi)) * coeff;
708 wi = wi + len; vi = vi + len;
713 Y.Update(1.0, W, 1.0);
720 ApplyInverseTime_ += Time_->ElapsedTime();
728 const Epetra_Vector& InvPointDiagonal,
729 const int MaximumIterations,
734 double RQ_top, RQ_bottom, norm;
735 Epetra_Vector x(Operator.OperatorDomainMap());
736 Epetra_Vector y(Operator.OperatorRangeMap());
739 if (norm == 0.0) IFPACK_CHK_ERR(-1);
743 for (
int iter = 0; iter < MaximumIterations; ++iter)
745 Operator.Apply(x, y);
746 IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
747 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
748 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
749 lambda_max = RQ_top / RQ_bottom;
750 IFPACK_CHK_ERR(y.Norm2(&norm));
751 if (norm == 0.0) IFPACK_CHK_ERR(-1);
752 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
760 CG(
const Epetra_Operator& Operator,
761 const Epetra_Vector& InvPointDiagonal,
762 const int MaximumIterations,
763 double& lambda_min,
double& lambda_max)
765 #ifdef HAVE_IFPACK_AZTECOO 766 Epetra_Vector x(Operator.OperatorDomainMap());
767 Epetra_Vector y(Operator.OperatorRangeMap());
771 Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
773 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
774 solver.SetAztecOption(AZ_output, AZ_none);
777 Operator.OperatorRangeMap(),
779 solver.SetPrecOperator(&diag);
780 solver.Iterate(MaximumIterations, 1e-10);
782 const double* status = solver.GetAztecStatus();
784 lambda_min = status[AZ_lambda_min];
785 lambda_max = status[AZ_lambda_max];
792 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
793 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
794 cout <<
"in your configure script." << endl;
800 #ifdef HAVE_IFPACK_EPETRAEXT 802 PowerMethod(
const int MaximumIterations,
double& lambda_max)
805 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
808 double RQ_top, RQ_bottom, norm;
809 Epetra_Vector x(Operator_->OperatorDomainMap());
810 Epetra_Vector y(Operator_->OperatorRangeMap());
811 Epetra_Vector z(Operator_->OperatorRangeMap());
814 if (norm == 0.0) IFPACK_CHK_ERR(-1);
818 for (
int iter = 0; iter < MaximumIterations; ++iter)
820 Operator_->Apply(x, z);
821 InvBlockDiagonal_->ApplyInverse(z,y);
822 if(SolveNormalEquations_){
823 InvBlockDiagonal_->ApplyInverse(y,z);
824 Apply_Transpose(Operator_,z, y);
827 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
828 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
829 lambda_max = RQ_top / RQ_bottom;
830 IFPACK_CHK_ERR(y.Norm2(&norm));
831 if (norm == 0.0) IFPACK_CHK_ERR(-1);
832 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
840 #ifdef HAVE_IFPACK_EPETRAEXT 842 CG(
const int MaximumIterations,
843 double& lambda_min,
double& lambda_max)
848 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
850 #ifdef HAVE_IFPACK_AZTECOO 851 Epetra_Vector x(Operator_->OperatorDomainMap());
852 Epetra_Vector y(Operator_->OperatorRangeMap());
855 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
858 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
859 solver.SetAztecOption(AZ_output, AZ_none);
861 solver.SetPrecOperator(&*InvBlockDiagonal_);
862 solver.Iterate(MaximumIterations, 1e-10);
864 const double* status = solver.GetAztecStatus();
866 lambda_min = status[AZ_lambda_min];
867 lambda_max = status[AZ_lambda_max];
874 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
875 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
876 cout <<
"in your configure script." << endl;
static int CG(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_min, double &lambda_max)
Uses AztecOO's CG to estimate lambda_min and lambda_max.
virtual int Compute()
Computes the preconditioners.
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 bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
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.
static int PowerMethod(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &LambdaMax)
Simple power method to compute lambda_max.
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.
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.