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