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_Polynomial.h"
54#include "Ifpack_Condest.h"
55#include "Teuchos_LAPACK.hpp"
56#include "Teuchos_SerialDenseMatrix.hpp"
58#ifdef HAVE_IFPACK_AZTECOO
59#include "Ifpack_DiagPreconditioner.h"
63#ifdef HAVE_IFPACK_EPETRAEXT
64#include "Epetra_CrsMatrix.h"
65#include "EpetraExt_PointToBlockDiagPermute.h"
69#define ABS(x) ((x)>0?(x):-(x))
87 IsInitialized_(false),
96 ApplyInverseTime_(0.0),
98 ApplyInverseFlops_(0.0),
102 UseTranspose_(false),
109 LambdaRealMax_(-1.0),
112 MinDiagonalValue_(0.0),
116 NumGlobalNonzeros_(0),
117 Operator_(Teuchos::rcp(Operator,false)),
118 UseBlockMode_(false),
119 SolveNormalEquations_(false),
121 ZeroStartingSolution_(true)
136 IsInitialized_(false),
138 IsIndefinite_(false),
143 InitializeTime_(0.0),
145 ApplyInverseTime_(0.0),
147 ApplyInverseFlops_(0.0),
151 UseTranspose_(false),
159 LambdaRealMax_(-1.0),
162 MinDiagonalValue_(0.0),
166 NumGlobalNonzeros_(0),
167 Operator_(Teuchos::rcp(Operator,false)),
168 Matrix_(Teuchos::rcp(Operator,false)),
169 UseBlockMode_(false),
170 SolveNormalEquations_(false),
172 ZeroStartingSolution_(true)
180 RealEigRatio_ = List.get(
"polynomial: real eigenvalue ratio", RealEigRatio_);
181 ImagEigRatio_ = List.get(
"polynomial: imag eigenvalue ratio", ImagEigRatio_);
182 LambdaRealMin_ = List.get(
"polynomial: min real part", LambdaRealMin_);
183 LambdaRealMax_ = List.get(
"polynomial: max real part", LambdaRealMax_);
184 LambdaImagMin_ = List.get(
"polynomial: min imag part", LambdaImagMin_);
185 LambdaImagMax_ = List.get(
"polynomial: max imag part", LambdaImagMax_);
186 PolyDegree_ = List.get(
"polynomial: degree",PolyDegree_);
187 LSPointsReal_ = List.get(
"polynomial: real interp points",LSPointsReal_);
188 LSPointsImag_ = List.get(
"polynomial: imag interp points",LSPointsImag_);
189 IsIndefinite_ = List.get(
"polynomial: indefinite",IsIndefinite_);
190 IsComplex_ = List.get(
"polynomial: complex",IsComplex_);
191 MinDiagonalValue_ = List.get(
"polynomial: min diagonal value",
193 ZeroStartingSolution_ = List.get(
"polynomial: zero starting solution",
194 ZeroStartingSolution_);
196 Epetra_Vector* ID = List.get(
"polynomial: operator inv diagonal",
198 EigMaxIters_ = List.get(
"polynomial: eigenvalue max iterations",EigMaxIters_);
200#ifdef HAVE_IFPACK_EPETRAEXT
202 UseBlockMode_ = List.get(
"polynomial: use block mode",UseBlockMode_);
203 if(!List.isParameter(
"polynomial: block list")) UseBlockMode_=
false;
205 BlockList_ = List.get(
"polynomial: block list",BlockList_);
209 Teuchos::ParameterList Blist;
210 Blist=BlockList_.get(
"blockdiagmatrix: list",Blist);
211 std::string dummy(
"invert");
212 std::string ApplyMode=Blist.get(
"apply mode",dummy);
213 if(ApplyMode==std::string(
"multiply")){
214 Blist.set(
"apply mode",
"invert");
215 BlockList_.set(
"blockdiagmatrix: list",Blist);
220 SolveNormalEquations_ = List.get(
"polynomial: solve normal equations",SolveNormalEquations_);
235 return(Operator_->Comm());
241 return(Operator_->OperatorDomainMap());
247 return(Operator_->OperatorRangeMap());
266 IFPACK_CHK_ERR(Operator_->Apply(X,Y));
275 IsInitialized_ =
false;
277 if (Operator_ == Teuchos::null)
280 if (Time_ == Teuchos::null)
285 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
288 NumMyRows_ = Matrix_->NumMyRows();
289 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
290 NumGlobalRows_ = Matrix_->NumGlobalRows64();
291 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
295 if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
296 Operator_->OperatorRangeMap().NumGlobalElements64())
301 InitializeTime_ += Time_->ElapsedTime();
302 IsInitialized_ =
true;
312 Time_->ResetStartTime();
318 if (PolyDegree_ <= 0)
321#ifdef HAVE_IFPACK_EPETRAEXT
323 if(IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
327 if(!CrsMatrix) UseBlockMode_=
false;
330 InvBlockDiagonal_=Teuchos::rcp(
new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
331 if(InvBlockDiagonal_==Teuchos::null) IFPACK_CHK_ERR(-6);
333 ierr=InvBlockDiagonal_->SetParameters(BlockList_);
334 if(ierr) IFPACK_CHK_ERR(-7);
336 ierr=InvBlockDiagonal_->Compute();
337 if(ierr) IFPACK_CHK_ERR(-8);
341 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_)
345 if (InvDiagonal_ == Teuchos::null)
348 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*InvDiagonal_));
352 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
353 double diag = (*InvDiagonal_)[i];
354 if (IFPACK_ABS(diag) < MinDiagonalValue_)
355 (*InvDiagonal_)[i] = MinDiagonalValue_;
357 (*InvDiagonal_)[i] = 1.0 / diag;
362 double lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max;
363 if (LambdaRealMax_ == -1) {
365 GMRES(
Matrix(), *InvDiagonal_, EigMaxIters_, lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max);
366 LambdaRealMin_=lambda_real_min; LambdaImagMin_=lambda_imag_min;
367 LambdaRealMax_=lambda_real_max; LambdaImagMax_=lambda_imag_max;
378 const std::complex<double> zero(0.0,0.0);
381 double lenx = LambdaRealMax_-LambdaRealMin_;
382 int nx = ceil(lenx*((
double) LSPointsReal_));
383 if (nx<2) { nx = 2; }
384 double hx = lenx/((double) nx);
385 std::vector<double> xs;
386 if(abs(lenx)>1.0e-8) {
387 for(
int pt=0; pt<=nx; pt++ ) {
388 xs.push_back(hx*pt+LambdaRealMin_);
392 xs.push_back(LambdaRealMax_);
395 double leny = LambdaImagMax_-LambdaImagMin_;
396 int ny = ceil(leny*((
double) LSPointsImag_));
397 if (ny<2) { ny = 2; }
398 double hy = leny/((double) ny);
399 std::vector<double> ys;
400 if(abs(leny)>1.0e-8) {
401 for(
int pt=0; pt<=ny; pt++ ) {
402 ys.push_back(hy*pt+LambdaImagMin_);
406 ys.push_back(LambdaImagMax_);
409 std::vector< std::complex<double> > cpts;
410 for(
int jj=0; jj<ny; jj++ ) {
411 for(
int ii=0; ii<nx; ii++ ) {
412 std::complex<double> cpt(xs[ii],ys[jj]);
416 cpts.push_back(zero);
418#ifdef HAVE_TEUCHOS_COMPLEX
419 const std::complex<double> one(1.0,0.0);
422 Teuchos::SerialDenseMatrix<int, std::complex<double> > Vmatrix(cpts.size(),PolyDegree_+1);
423 Vmatrix.putScalar(zero);
424 for (
int jj = 0; jj <= PolyDegree_; ++jj) {
425 for (
int ii = 0; ii < static_cast<int> (cpts.size ()) - 1; ++ii) {
427 Vmatrix(ii,jj) = pow(cpts[ii],jj);
430 Vmatrix(ii,jj) = one;
434 Vmatrix(cpts.size()-1,0)=one;
437 Teuchos::SerialDenseMatrix< int,std::complex<double> > RHS(cpts.size(),1);
439 RHS(cpts.size()-1,0)=one;
442 Teuchos::LAPACK< int, std::complex<double> > lapack;
443 const int N = Vmatrix.numCols();
444 Teuchos::Array<double> singularValues(N);
445 Teuchos::Array<double> rwork(1);
446 rwork.resize (std::max (1, 5 * N));
447 std::complex<double> lworkScalar(1.0,0.0);
449 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
450 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
451 &lworkScalar, -1, &info);
452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
453 "_GELSS workspace query returned INFO = "
454 << info <<
" != 0.");
455 const int lwork =
static_cast<int> (real(lworkScalar));
456 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
457 "_GELSS workspace query returned LWORK = "
458 << lwork <<
" < 0.");
460 Teuchos::Array< std::complex<double> > work (std::max (1, lwork));
462 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
463 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
464 &work[0], lwork, &info);
466 coeff_.resize(PolyDegree_+1);
467 std::complex<double> c0=RHS(0,0);
468 for(
int ii=0; ii<=PolyDegree_; ii++) {
473 coeff_[ii]=real(RHS(ii,0)/c0);
480 Teuchos::SerialDenseMatrix< int, double > Vmatrix(xs.size()+1,PolyDegree_+1);
481 Vmatrix.putScalar(0.0);
482 for(
int jj=0; jj<=PolyDegree_; jj++) {
483 for( std::vector<double>::size_type ii=0; ii<xs.size(); ii++) {
485 Vmatrix(ii,jj)=pow(xs[ii],jj);
492 Vmatrix(xs.size(),0)=1.0;
495 Teuchos::SerialDenseMatrix< int, double > RHS(xs.size()+1,1);
497 RHS(xs.size(),0)=1.0;
500 Teuchos::LAPACK< int, double > lapack;
501 const int N = Vmatrix.numCols();
502 Teuchos::Array<double> singularValues(N);
503 Teuchos::Array<double> rwork(1);
504 rwork.resize (std::max (1, 5 * N));
505 double lworkScalar(1.0);
507 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
508 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
509 &lworkScalar, -1, &info);
510 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
511 "_GELSS workspace query returned INFO = "
512 << info <<
" != 0.");
513 const int lwork =
static_cast<int> (lworkScalar);
514 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
515 "_GELSS workspace query returned LWORK = "
516 << lwork <<
" < 0.");
518 Teuchos::Array< double > work (std::max (1, lwork));
520 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
521 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
522 &work[0], lwork, &info);
524 coeff_.resize(PolyDegree_+1);
526 for(
int ii=0; ii<=PolyDegree_; ii++) {
531 coeff_[ii]=RHS(ii,0)/c0;
536#ifdef IFPACK_FLOPCOUNTERS
537 ComputeFlops_ += NumMyRows_;
541 ComputeTime_ += Time_->ElapsedTime();
552 double MyMinVal, MyMaxVal;
553 double MinVal, MaxVal;
556 InvDiagonal_->MinValue(&MyMinVal);
557 InvDiagonal_->MaxValue(&MyMaxVal);
562 if (!
Comm().MyPID()) {
564 os <<
"================================================================================" << endl;
565 os <<
"Ifpack_Polynomial" << endl;
566 os <<
"Degree of polynomial = " << PolyDegree_ << endl;
567 os <<
"Condition number estimate = " <<
Condest() << endl;
568 os <<
"Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
570 os <<
"Minimum value on stored inverse diagonal = " << MinVal << endl;
571 os <<
"Maximum value on stored inverse diagonal = " << MaxVal << endl;
573 if (ZeroStartingSolution_)
574 os <<
"Using zero starting solution" << endl;
576 os <<
"Using input starting solution" << endl;
578 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
579 os <<
"----- ------- -------------- ------------ --------" << endl;
580 os <<
"Initialize() " << std::setw(5) << NumInitialize_
581 <<
" " << std::setw(15) << InitializeTime_
582 <<
" 0.0 0.0" << endl;
583 os <<
"Compute() " << std::setw(5) << NumCompute_
584 <<
" " << std::setw(15) << ComputeTime_
585 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
586 if (ComputeTime_ != 0.0)
587 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
589 os <<
" " << std::setw(15) << 0.0 << endl;
590 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
591 <<
" " << std::setw(15) << ApplyInverseTime_
592 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
593 if (ApplyInverseTime_ != 0.0)
594 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
596 os <<
" " << std::setw(15) << 0.0 << endl;
597 os <<
"================================================================================" << endl;
606Condest(
const Ifpack_CondestType CT,
607 const int MaxIters,
const double Tol,
615 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
621void Ifpack_Polynomial::SetLabel()
623 Label_ =
"IFPACK (Least squares polynomial), degree=" + Ifpack_toString(PolyDegree_);
634 if (PolyDegree_ == 0)
641 Time_->ResetStartTime();
644 if(ZeroStartingSolution_==
true) {
655 Y.
Update(-coeff_[1], Xcopy, 1.0);
656 for (
int ii = 2; ii < static_cast<int> (coeff_.size ()); ++ii) {
658 Operator_->Apply(V,Xcopy);
659 Xcopy.
Multiply(1.0, *InvDiagonal_, Xcopy, 0.0);
661 Y.
Update(-coeff_[ii], Xcopy, 1.0);
666 ApplyInverseTime_ += Time_->ElapsedTime();
674 const int MaximumIterations,
679 double RQ_top, RQ_bottom, norm;
684 if (norm == 0.0) IFPACK_CHK_ERR(-1);
688 for (
int iter = 0; iter < MaximumIterations; ++iter)
690 Operator.
Apply(x, y);
691 IFPACK_CHK_ERR(y.
Multiply(1.0, InvPointDiagonal, y, 0.0));
692 IFPACK_CHK_ERR(y.
Dot(x, &RQ_top));
693 IFPACK_CHK_ERR(x.
Dot(x, &RQ_bottom));
694 lambda_max = RQ_top / RQ_bottom;
695 IFPACK_CHK_ERR(y.
Norm2(&norm));
696 if (norm == 0.0) IFPACK_CHK_ERR(-1);
697 IFPACK_CHK_ERR(x.
Update(1.0 / norm, y, 0.0));
707 const int MaximumIterations,
708 double& lambda_min,
double& lambda_max)
710#ifdef HAVE_IFPACK_AZTECOO
718 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
719 solver.SetAztecOption(AZ_output, AZ_none);
724 solver.SetPrecOperator(&diag);
725 solver.Iterate(MaximumIterations, 1e-10);
727 const double* status = solver.GetAztecStatus();
729 lambda_min = status[AZ_lambda_min];
730 lambda_max = status[AZ_lambda_max];
737 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
738 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
739 cout <<
"in your configure script." << endl;
745#ifdef HAVE_IFPACK_EPETRAEXT
747PowerMethod(
const int MaximumIterations,
double& lambda_max)
750 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
753 double RQ_top, RQ_bottom, norm;
759 if (norm == 0.0) IFPACK_CHK_ERR(-1);
763 for (
int iter = 0; iter < MaximumIterations; ++iter)
765 Operator_->Apply(x, z);
766 InvBlockDiagonal_->ApplyInverse(z,y);
767 if(SolveNormalEquations_){
768 InvBlockDiagonal_->ApplyInverse(y,z);
769 Apply_Transpose(Operator_,z, y);
772 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
773 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
774 lambda_max = RQ_top / RQ_bottom;
775 IFPACK_CHK_ERR(y.Norm2(&norm));
776 if (norm == 0.0) IFPACK_CHK_ERR(-1);
777 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
785#ifdef HAVE_IFPACK_EPETRAEXT
798 if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
800#ifdef HAVE_IFPACK_AZTECOO
808 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
809 solver.SetAztecOption(AZ_output, AZ_none);
811 solver.SetPrecOperator(&*InvBlockDiagonal_);
812 solver.Iterate(MaximumIterations, 1e-10);
814 const double* status = solver.GetAztecStatus();
816 lambda_min = status[AZ_lambda_min];
817 lambda_max = status[AZ_lambda_max];
824 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
825 cout <<
"to use the CG estimator. This may require --enable-aztecoo" << endl;
826 cout <<
"in your configure script." << endl;
838 const int MaximumIterations,
839 double& lambda_real_min,
double& lambda_real_max,
840 double& lambda_imag_min,
double& lambda_imag_max)
842#ifdef HAVE_IFPACK_AZTECOO
849 solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
850 solver.SetAztecOption(AZ_output, AZ_none);
854 solver.SetPrecOperator(&diag);
855 solver.Iterate(MaximumIterations, 1e-10);
856 const double* status = solver.GetAztecStatus();
857 lambda_real_min = status[AZ_lambda_real_min];
858 lambda_real_max = status[AZ_lambda_real_max];
859 lambda_imag_min = status[AZ_lambda_imag_min];
860 lambda_imag_max = status[AZ_lambda_imag_max];
866 cout <<
"You need to configure IFPACK with support for AztecOO" << endl;
867 cout <<
"to use the GMRES estimator. This may require --enable-aztecoo" << endl;
868 cout <<
"in your configure script." << endl;
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
int Scale(double ScalarValue)
int Dot(const Epetra_MultiVector &A, double *Result) const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Norm2(double *Result) const
int PutScalar(double ScalarConstant)
virtual int SetUseTranspose(bool UseTranspose)=0
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Ifpack_Polynomial(const Epetra_Operator *Matrix)
Ifpack_Polynomial constructor with given Epetra_Operator/Epetra_RowMatrix.
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 int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Compute()
Computes the preconditioners.
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 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.
int GMRES(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_real_min, double &lambda_real_max, double &lambda_imag_min, double &lambda_imag_max)
Uses AztecOO's GMRES to estimate the height and width of the spectrum.
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 bool UseTranspose() const
Returns the current UseTranspose setting.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.