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 #include "Teuchos_LAPACK.hpp" 56 #include "Teuchos_SerialDenseMatrix.hpp" 58 #ifdef HAVE_IFPACK_AZTECOO 63 #ifdef HAVE_IFPACK_EPETRAEXT 64 #include "Epetra_CrsMatrix.h" 65 #include "EpetraExt_PointToBlockDiagPermute.h" 69 #define ABS(x) ((x)>0?(x):-(x)) 72 inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,
const Epetra_MultiVector &X,Epetra_MultiVector &Y){
73 Epetra_Operator * Operator=
const_cast<Epetra_Operator*
>(&*Operator_);
74 Operator->SetUseTranspose(
true);
76 Operator->SetUseTranspose(
false);
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),
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),
169 UseBlockMode_(
false),
170 SolveNormalEquations_(
false),
172 ZeroStartingSolution_(
true)
196 Epetra_Vector* ID = List.get(
"polynomial: operator inv diagonal",
200 #ifdef HAVE_IFPACK_EPETRAEXT 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);
252 Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 257 if (X.NumVectors() != Y.NumVectors())
280 if (
Time_ == Teuchos::null)
281 Time_ = Teuchos::rcp(
new Epetra_Time(
Comm()) );
285 if (
Matrix().NumGlobalRows64() !=
Matrix().NumGlobalCols64())
295 if (
Operator_->OperatorDomainMap().NumGlobalElements64() !=
296 Operator_->OperatorRangeMap().NumGlobalElements64())
312 Time_->ResetStartTime();
321 #ifdef HAVE_IFPACK_EPETRAEXT 324 const Epetra_CrsMatrix *CrsMatrix=
dynamic_cast<const Epetra_CrsMatrix*
>(&*
Matrix_);
330 InvBlockDiagonal_=Teuchos::rcp(
new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
333 ierr=InvBlockDiagonal_->SetParameters(BlockList_);
336 ierr=InvBlockDiagonal_->Compute();
353 double diag = (*InvDiagonal_)[i];
362 double lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max;
378 const std::complex<double> zero(0.0,0.0);
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++ ) {
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++ ) {
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);
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);
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);
467 std::complex<double> c0=
RHS(0,0);
480 Teuchos::SerialDenseMatrix< int, double > Vmatrix(xs.size()+1,
PolyDegree_+1);
481 Vmatrix.putScalar(0.0);
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);
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);
536 #ifdef IFPACK_FLOPCOUNTERS 552 double MyMinVal, MyMaxVal;
553 double MinVal, MaxVal;
558 Comm().MinAll(&MyMinVal,&MinVal,1);
559 Comm().MinAll(&MyMaxVal,&MaxVal,1);
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;
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;
582 <<
" 0.0 0.0" << endl;
589 os <<
" " << std::setw(15) << 0.0 << endl;
596 os <<
" " << std::setw(15) << 0.0 << endl;
597 os <<
"================================================================================" << endl;
607 const int MaxIters,
const double Tol,
608 Epetra_RowMatrix* Matrix_in)
637 int nVec = X.NumVectors();
638 if (nVec != Y.NumVectors())
641 Time_->ResetStartTime();
643 Epetra_MultiVector Xcopy(X);
655 Y.Update(-
coeff_[1], Xcopy, 1.0);
656 for (
int ii = 2; ii < static_cast<int> (
coeff_.size ()); ++ii) {
657 const Epetra_MultiVector V(Xcopy);
661 Y.Update(-
coeff_[ii], Xcopy, 1.0);
673 const Epetra_Vector& InvPointDiagonal,
674 const int MaximumIterations,
679 double RQ_top, RQ_bottom, norm;
680 Epetra_Vector x(Operator.OperatorDomainMap());
681 Epetra_Vector y(Operator.OperatorRangeMap());
688 for (
int iter = 0; iter < MaximumIterations; ++iter)
690 Operator.Apply(x, y);
694 lambda_max = RQ_top / RQ_bottom;
705 CG(
const Epetra_Operator& Operator,
706 const Epetra_Vector& InvPointDiagonal,
707 const int MaximumIterations,
708 double& lambda_min,
double& lambda_max)
710 #ifdef HAVE_IFPACK_AZTECOO 711 Epetra_Vector x(Operator.OperatorDomainMap());
712 Epetra_Vector y(Operator.OperatorRangeMap());
716 Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
718 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
719 solver.SetAztecOption(AZ_output, AZ_none);
722 Operator.OperatorRangeMap(),
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 747 PowerMethod(
const int MaximumIterations,
double& lambda_max)
753 double RQ_top, RQ_bottom, norm;
754 Epetra_Vector x(
Operator_->OperatorDomainMap());
755 Epetra_Vector y(
Operator_->OperatorRangeMap());
756 Epetra_Vector z(
Operator_->OperatorRangeMap());
763 for (
int iter = 0; iter < MaximumIterations; ++iter)
766 InvBlockDiagonal_->ApplyInverse(z,y);
768 InvBlockDiagonal_->ApplyInverse(y,z);
774 lambda_max = RQ_top / RQ_bottom;
785 #ifdef HAVE_IFPACK_EPETRAEXT 787 CG(
const int MaximumIterations,
788 double& lambda_min,
double& lambda_max)
800 #ifdef HAVE_IFPACK_AZTECOO 801 Epetra_Vector x(
Operator_->OperatorDomainMap());
802 Epetra_Vector y(
Operator_->OperatorRangeMap());
805 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*
Matrix_), &x, &y);
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;
832 #endif // HAVE_IFPACK_EPETRAEXT 836 GMRES(
const Epetra_Operator& Operator,
837 const Epetra_Vector& InvPointDiagonal,
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 843 Epetra_Vector x(
Operator_->OperatorDomainMap());
844 Epetra_Vector y(
Operator_->OperatorRangeMap());
847 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*
Matrix_), &x, &y);
849 solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
850 solver.SetAztecOption(AZ_output, AZ_none);
852 Operator.OperatorRangeMap(),
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;
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.
bool IsInitialized_
If true, the preconditioner has been computed successfully.
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
Teuchos::RefCountPtr< Epetra_Vector > InvDiagonal_
Contains the inverse of diagonal elements of Matrix.
int NumCompute_
Contains the number of successful call to Compute().
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
double Condest_
Contains the estimated condition number.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
long long NumGlobalNonzeros_
Number of global nonzeros.
int LSPointsReal_
Contains the number of discretization points of the least squares problem.
double RealEigRatio_
Contains the ratio such that the rectangular domain in the complex plane is [-LambdaRealMax_ / EigRat...
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
long long NumGlobalRows_
Number of global rows.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
double ComputeFlops_
Contains the number of flops for Compute().
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned as an Epetra_RowMatrix.
std::string Label_
Contains the label of this object.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
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.
void Apply_Transpose(Teuchos::RCP< const Epetra_Operator > Operator_, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
std::vector< double > coeff_
coefficients of the polynomial
int PolyDegree_
Contains the degree of the least squares polynomial.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Teuchos::RefCountPtr< const Epetra_Operator > Operator_
Pointers to the matrix to be preconditioned as an Epetra_Operator.
double InitializeTime_
Contains the time for all successful calls to Initialize().
double MinDiagonalValue_
Contains the minimum value on the diagonal.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
Ifpack_Polynomial(const Epetra_Operator *Matrix)
Ifpack_Polynomial 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.
int NumMyNonzeros_
Number of local nonzeros.
double ComputeTime_
Contains the time for all successful calls to Compute().
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
bool IsComputed_
If true, the preconditioner has been computed successfully.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
int NumMyRows_
Number of local rows.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
bool IsIndefinite_
If true, have to compute polynomial for a spectrum with negative eigenvalues.
int EigMaxIters_
Max number of iterations to use in eigenvalue estimation (if automatic).
bool SolveNormalEquations_
Run on the normal equations.
static int PowerMethod(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &LambdaMax)
Simple power method to compute lambda_max.
bool UseBlockMode_
Use Block Preconditioning.
bool IsRowMatrix_
If true, the Operator_ is an Epetra_RowMatrix.
double LambdaRealMin_
Bounds on the spectrum.
int NumInitialize_
Contains the number of successful calls to Initialize().
#define IFPACK_CHK_ERR(ifpack_err)
bool IsComplex_
If true, have to compute polynomial for a spectrum with nonzero imaginary part.
virtual void SetLabel()
Sets the label.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.