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());
257 if (X.NumVectors() != Y.NumVectors())
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;
607 const int MaxIters,
const double Tol,
615 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
621 void Ifpack_Polynomial::SetLabel()
623 Label_ =
"IFPACK (Least squares polynomial), degree=" +
Ifpack_toString(PolyDegree_);
634 if (PolyDegree_ == 0)
637 int nVec = X.NumVectors();
638 if (nVec != Y.NumVectors())
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
747 PowerMethod(
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;
832 #endif // HAVE_IFPACK_EPETRAEXT
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;
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.
virtual int SetUseTranspose(bool UseTranspose)=0
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
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
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 Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual int Compute()
Computes the preconditioners.
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 SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
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_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.