Ifpack_Amesos.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #if defined(HAVE_IFPACK_AMESOS) && defined(HAVE_IFPACK_TEUCHOS)
00003 #include "Ifpack_Preconditioner.h"
00004 #include "Ifpack_Amesos.h"
00005 #include "Ifpack_Condest.h"
00006 #include "Epetra_MultiVector.h"
00007 #include "Epetra_Map.h"
00008 #include "Epetra_Comm.h"
00009 #include "Amesos.h"
00010 #include "Epetra_LinearProblem.h"
00011 #include "Epetra_RowMatrix.h"
00012 #include "Epetra_Time.h"
00013 #include "Teuchos_ParameterList.hpp"
00014 
00015 static bool FirstTime = true;
00016 
00017 //==============================================================================
00018 Ifpack_Amesos::Ifpack_Amesos(Epetra_RowMatrix* Matrix) :
00019   Matrix_(Matrix),
00020   Solver_(0),
00021   Problem_(0),
00022   Label_("Amesos_Klu"),
00023   IsInitialized_(false),
00024   IsComputed_(false),
00025   UseTranspose_(false),
00026   NumInitialize_(0),
00027   NumCompute_(0),
00028   NumApplyInverse_(0),
00029   InitializeTime_(0.0),
00030   ComputeTime_(0.0),
00031   ApplyInverseTime_(0.0),
00032   Time_(0),
00033   ComputeFlops_(0),
00034   ApplyInverseFlops_(0),
00035   Condest_(-1.0)
00036 {
00037   Problem_ = new Epetra_LinearProblem;
00038 }
00039 
00040 //==============================================================================
00041 Ifpack_Amesos::Ifpack_Amesos(const Ifpack_Amesos& rhs) :
00042   Matrix_(&rhs.Matrix()),
00043   Solver_(0),
00044   Problem_(0),
00045   Label_(rhs.Label()),
00046   IsInitialized_(false),
00047   IsComputed_(false),
00048   NumInitialize_(rhs.NumInitialize()),
00049   NumCompute_(rhs.NumCompute()),
00050   NumApplyInverse_(rhs.NumApplyInverse()),
00051   InitializeTime_(rhs.InitializeTime()),
00052   ComputeTime_(rhs.ComputeTime()),
00053   ApplyInverseTime_(rhs.ApplyInverseTime()),
00054   Time_(0),
00055   ComputeFlops_(rhs.ComputeFlops()),
00056   ApplyInverseFlops_(rhs.ApplyInverseFlops()),
00057   Condest_(rhs.Condest())
00058 {
00059 
00060   Problem_ = new Epetra_LinearProblem;
00061 
00062   // copy the RHS list in *this.List
00063   Teuchos::ParameterList RHSList(rhs.List());
00064   List_ = RHSList;
00065 
00066   // I do not have a copy constructor for Amesos,
00067   // so Initialize() and Compute() of this object 
00068   // are called if the rhs did so
00069   if (rhs.IsInitialized()) {
00070     IsInitialized_ = true;
00071     Initialize();
00072   }
00073   if (rhs.IsComputed()) {
00074     IsComputed_ = true;
00075     Compute();
00076   }
00077 
00078 }
00079 //==============================================================================
00080 Ifpack_Amesos::~Ifpack_Amesos()
00081 {
00082   if (Problem_)
00083     delete Problem_;
00084 
00085   if (Solver_)
00086     delete Solver_;
00087 
00088   if (Time_)
00089     delete Time_;
00090 }
00091 
00092 //==============================================================================
00093 int Ifpack_Amesos::SetParameters(Teuchos::ParameterList& List)
00094 {
00095 
00096   List_ = List;
00097   Label_ = List.get("amesos: solver type", Label_);
00098   return(0);
00099 }
00100 
00101 //==============================================================================
00102 int Ifpack_Amesos::Initialize()
00103 {
00104 
00105   IsInitialized_ = false;
00106   IsComputed_ = false;
00107 
00108   if (Matrix_ == 0)
00109     IFPACK_CHK_ERR(-1);
00110 
00111 #if 0
00112   // better to avoid strange games with maps, this class should be
00113   // used for Ifpack_LocalFilter'd matrices only
00114   if (Comm().NumProc() != 1) {
00115     cout << "Class Ifpack_Amesos must be used for serial runs;" << endl;
00116     cout << "for parallel runs you should declare objects as:" << endl; 
00117     cout << "Ifpack_AdditiveSchwarz<Ifpack_Amesos> APrec(Matrix)" << endl;
00118     exit(EXIT_FAILURE);
00119   }
00120 #endif
00121 
00122   // only square matrices
00123   if (Matrix_->NumGlobalRows() != Matrix_->NumGlobalCols())
00124     IFPACK_CHK_ERR(-1);
00125 
00126   // at least one nonzero
00127   if (Matrix_->NumMyNonzeros() == 0) 
00128     IFPACK_CHK_ERR(-1);
00129 
00130   Problem_->SetOperator(const_cast<Epetra_RowMatrix*>(Matrix_));
00131 
00132   if (Time_ == 0)
00133     Time_ = new Epetra_Time(Comm());
00134 
00135   // reallocate the solver. 
00136   if (Solver_)
00137     delete Solver_;
00138 
00139   Amesos Factory;
00140   Solver_ = Factory.Create((char*)Label_.c_str(),*Problem_);
00141   
00142   if (Solver_ == 0) 
00143   {
00144     // try to create KLU, it is generally enabled
00145     Solver_ = Factory.Create("Amesos_Klu",*Problem_);
00146   }
00147   if (Solver_ == 0)
00148   {
00149     // finally try to create LAPACK, it is generally enabled
00150     // NOTE: the matrix is dense, so this should only be for
00151     // small problems...
00152     if (FirstTime)
00153     {
00154       cerr << "IFPACK WARNING: In class Ifpack_Amesos:" << endl;
00155       cerr << "IFPACK WARNING: Using LAPACK because other Amesos" << endl;
00156       cerr << "IFPACK WARNING: solvers are not available. LAPACK" << endl;
00157       cerr << "IFPACK WARNING: allocates memory to store the matrix as" << endl;
00158       cerr << "IFPACK WARNING: dense, I hope you have enough memory..." << endl;
00159       cerr << "IFPACK WARNING: (file " << __FILE__ << ", line " << __LINE__ 
00160            << ")" << endl;
00161       FirstTime = false;
00162     }
00163     Solver_ = Factory.Create("Amesos_Lapack",*Problem_);
00164   }
00165   // if empty, give up.
00166   if (Solver_ == 0)
00167     IFPACK_CHK_ERR(-1);
00168 
00169   IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_));
00170   Solver_->SetParameters(List_);
00171   IFPACK_CHK_ERR(Solver_->SymbolicFactorization());
00172 
00173   IsInitialized_ = true;
00174   ++NumInitialize_;
00175   InitializeTime_ += Time_->ElapsedTime();
00176   return(0);
00177 }
00178 
00179 //==============================================================================
00180 int Ifpack_Amesos::Compute()
00181 {
00182 
00183   if (!IsInitialized())
00184     IFPACK_CHK_ERR(Initialize());
00185 
00186   IsComputed_ = false;
00187   Time_->ResetStartTime();
00188 
00189   if (Matrix_ == 0)
00190     IFPACK_CHK_ERR(-1);
00191 
00192   IFPACK_CHK_ERR(Solver_->NumericFactorization());
00193 
00194   IsComputed_ = true;
00195   ++NumCompute_;
00196   ComputeTime_ += Time_->ElapsedTime();
00197   return(0);
00198 }
00199 
00200 //==============================================================================
00201 int Ifpack_Amesos::SetUseTranspose(bool UseTranspose)
00202 {
00203   // store the value in UseTranspose_. If we have the solver, we pass to it
00204   // right away, otherwise we wait till when it is created.
00205   UseTranspose_ = UseTranspose;
00206   if (Solver_ != 0)
00207     IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose));
00208 
00209   return(0);
00210 }
00211 
00212 //==============================================================================
00213 int Ifpack_Amesos::
00214 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00215 {
00216   // check for maps ???
00217   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00218   return(0);
00219 }
00220 
00221 //==============================================================================
00222 int Ifpack_Amesos::
00223 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00224 {
00225 
00226   if (IsComputed() == false)
00227     IFPACK_CHK_ERR(-1);
00228 
00229   if (X.NumVectors() != Y.NumVectors())
00230     IFPACK_CHK_ERR(-1); // wrong input
00231   
00232   Time_->ResetStartTime();
00233 
00234   // AztecOO gives X and Y pointing to the same memory location,
00235   // need to create an auxiliary vector, Xcopy
00236   const Epetra_MultiVector* Xcopy;
00237   if (X.Pointers()[0] == Y.Pointers()[0])
00238     Xcopy = new Epetra_MultiVector(X);
00239   else
00240     Xcopy = &X;
00241     
00242   Problem_->SetLHS(&Y);
00243   Problem_->SetRHS((Epetra_MultiVector*)Xcopy);
00244   IFPACK_CHK_ERR(Solver_->Solve());
00245 
00246   if (Xcopy != &X)
00247     delete Xcopy;
00248 
00249   ++NumApplyInverse_;
00250   ApplyInverseTime_ += Time_->ElapsedTime();
00251 
00252   return(0);
00253 }
00254 
00255 //==============================================================================
00256 double Ifpack_Amesos::NormInf() const
00257 {
00258   return(-1.0);
00259 }
00260 
00261 //==============================================================================
00262 const char* Ifpack_Amesos::Label() const
00263 {
00264   return((char*)Label_.c_str());
00265 }
00266 
00267 //==============================================================================
00268 bool Ifpack_Amesos::UseTranspose() const
00269 {
00270   return(UseTranspose_);
00271 }
00272 
00273 //==============================================================================
00274 bool Ifpack_Amesos::HasNormInf() const
00275 {
00276   return(false);
00277 }
00278 
00279 //==============================================================================
00280 const Epetra_Comm & Ifpack_Amesos::Comm() const
00281 {
00282   return(Matrix_->Comm());
00283 }
00284 
00285 //==============================================================================
00286 const Epetra_Map & Ifpack_Amesos::OperatorDomainMap() const
00287 {
00288   return(Matrix_->OperatorDomainMap());
00289 }
00290 
00291 //==============================================================================
00292 const Epetra_Map & Ifpack_Amesos::OperatorRangeMap() const
00293 {
00294   return(Matrix_->OperatorRangeMap());
00295 }
00296 
00297 //==============================================================================
00298 double Ifpack_Amesos::Condest(const Ifpack_CondestType CT,
00299                               const int MaxIters, const double Tol,
00300                   Epetra_RowMatrix* Matrix)
00301 {
00302 
00303   if (!IsComputed()) // cannot compute right now
00304     return(-1.0);
00305 
00306   if (Condest_ == -1.0)
00307     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00308 
00309   return(Condest_);
00310 }
00311 
00312 //==============================================================================
00313 std::ostream& Ifpack_Amesos::Print(std::ostream& os) const
00314 {
00315   if (!Comm().MyPID()) {
00316     os << endl;
00317     os << "================================================================================" << endl;
00318     os << "Ifpack_Amesos: " << Label () << endl << endl;
00319     os << "Condition number estimate = " << Condest() << endl;
00320     os << "Global number of rows            = " << Matrix_->NumGlobalRows() << endl;
00321     os << endl;
00322     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00323     os << "-----           -------   --------------       ------------     --------" << endl;
00324     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00325        << "  " << std::setw(15) << InitializeTime_ 
00326        << "              0.0              0.0" << endl;
00327     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00328        << "  " << std::setw(15) << ComputeTime_
00329        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00330     if (ComputeTime_ != 0.0) 
00331       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00332     else
00333       os << "  " << std::setw(15) << 0.0 << endl;
00334     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00335        << "  " << std::setw(15) << ApplyInverseTime_
00336        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00337     if (ApplyInverseTime_ != 0.0) 
00338       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00339     else
00340       os << "  " << std::setw(15) << 0.0 << endl;
00341     os << "================================================================================" << endl;
00342     os << endl;
00343   }
00344 
00345   return(os);
00346 }
00347 #endif // HAVE_IFPACK_AMESOS && HAVE_IFPACK_TEUCHOS

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1