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