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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #include "Ifpack_ConfigDefs.h"
00044 #include "Ifpack_Preconditioner.h"
00045 #include "Ifpack_Amesos.h"
00046 #include "Ifpack_Condest.h"
00047 #include "Epetra_MultiVector.h"
00048 #include "Epetra_Map.h"
00049 #include "Epetra_Comm.h"
00050 #include "Amesos.h"
00051 #include "Epetra_LinearProblem.h"
00052 #include "Epetra_RowMatrix.h"
00053 #include "Epetra_Time.h"
00054 #include "Teuchos_ParameterList.hpp"
00055 
00056 static bool FirstTime = true;
00057 
00058 //==============================================================================
00059 Ifpack_Amesos::Ifpack_Amesos(Epetra_RowMatrix* Matrix_in) :
00060   Matrix_(Teuchos::rcp( Matrix_in, false )),
00061   Label_("Amesos_Klu"),
00062   IsEmpty_(false),
00063   IsInitialized_(false),
00064   IsComputed_(false),
00065   UseTranspose_(false),
00066   NumInitialize_(0),
00067   NumCompute_(0),
00068   NumApplyInverse_(0),
00069   InitializeTime_(0.0),
00070   ComputeTime_(0.0),
00071   ApplyInverseTime_(0.0),
00072   ComputeFlops_(0),
00073   ApplyInverseFlops_(0),
00074   Condest_(-1.0)
00075 {
00076   Problem_ = Teuchos::rcp( new Epetra_LinearProblem );
00077 }
00078 
00079 //==============================================================================
00080 Ifpack_Amesos::Ifpack_Amesos(const Ifpack_Amesos& rhs) :
00081   Matrix_(Teuchos::rcp( &rhs.Matrix(), false )),
00082   Label_(rhs.Label()),
00083   IsEmpty_(false),
00084   IsInitialized_(false),
00085   IsComputed_(false),
00086   NumInitialize_(rhs.NumInitialize()),
00087   NumCompute_(rhs.NumCompute()),
00088   NumApplyInverse_(rhs.NumApplyInverse()),
00089   InitializeTime_(rhs.InitializeTime()),
00090   ComputeTime_(rhs.ComputeTime()),
00091   ApplyInverseTime_(rhs.ApplyInverseTime()),
00092   ComputeFlops_(rhs.ComputeFlops()),
00093   ApplyInverseFlops_(rhs.ApplyInverseFlops()),
00094   Condest_(rhs.Condest())
00095 {
00096 
00097   Problem_ = Teuchos::rcp( new Epetra_LinearProblem );
00098 
00099   // copy the RHS list in *this.List
00100   Teuchos::ParameterList RHSList(rhs.List());
00101   List_ = RHSList;
00102 
00103   // I do not have a copy constructor for Amesos,
00104   // so Initialize() and Compute() of this object 
00105   // are called if the rhs did so
00106   if (rhs.IsInitialized()) {
00107     IsInitialized_ = true;
00108     Initialize();
00109   }
00110   if (rhs.IsComputed()) {
00111     IsComputed_ = true;
00112     Compute();
00113   }
00114 
00115 }
00116 //==============================================================================
00117 int Ifpack_Amesos::SetParameters(Teuchos::ParameterList& List_in)
00118 {
00119 
00120   List_ = List_in;
00121   Label_ = List_in.get("amesos: solver type", Label_);
00122   return(0);
00123 }
00124 
00125 //==============================================================================
00126 int Ifpack_Amesos::Initialize()
00127 {
00128 
00129   IsEmpty_ = false;
00130   IsInitialized_ = false;
00131   IsComputed_ = false;
00132 
00133   if (Matrix_ == Teuchos::null)
00134     IFPACK_CHK_ERR(-1);
00135 
00136 #if 0
00137   // better to avoid strange games with maps, this class should be
00138   // used for Ifpack_LocalFilter'd matrices only
00139   if (Comm().NumProc() != 1) {
00140     cout << "Class Ifpack_Amesos must be used for serial runs;" << endl;
00141     cout << "for parallel runs you should declare objects as:" << endl; 
00142     cout << "Ifpack_AdditiveSchwarz<Ifpack_Amesos> APrec(Matrix)" << endl;
00143     exit(EXIT_FAILURE);
00144   }
00145 #endif
00146 
00147   // only square matrices
00148   if (Matrix_->NumGlobalRows64() != Matrix_->NumGlobalCols64())
00149     IFPACK_CHK_ERR(-1);
00150 
00151   // if the matrix has a dimension of 0, this is an empty preconditioning object.
00152   if (Matrix_->NumGlobalRows64() == 0) {
00153     IsEmpty_ = true;
00154     IsInitialized_ = true;
00155     ++NumInitialize_;
00156     return(0);
00157   }
00158 
00159   Problem_->SetOperator(const_cast<Epetra_RowMatrix*>(Matrix_.get()));
00160 
00161   // create timer, which also starts it.
00162   if (Time_ == Teuchos::null)
00163     Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
00164 
00165   Amesos Factory;
00166   Solver_ = Teuchos::rcp( Factory.Create((char*)Label_.c_str(),*Problem_) );
00167   
00168   if (Solver_ == Teuchos::null) 
00169   {
00170     // try to create KLU, it is generally enabled
00171     Label_ = "Amesos_Klu";
00172     Solver_ = Teuchos::rcp( Factory.Create("Amesos_Klu",*Problem_) );
00173   }
00174   if (Solver_ == Teuchos::null)
00175   {
00176     // finally try to create LAPACK, it is generally enabled
00177     // NOTE: the matrix is dense, so this should only be for
00178     // small problems...
00179     if (FirstTime)
00180     {
00181       cerr << "IFPACK WARNING: In class Ifpack_Amesos:" << endl;
00182       cerr << "IFPACK WARNING: Using LAPACK because other Amesos" << endl;
00183       cerr << "IFPACK WARNING: solvers are not available. LAPACK" << endl;
00184       cerr << "IFPACK WARNING: allocates memory to store the matrix as" << endl;
00185       cerr << "IFPACK WARNING: dense, I hope you have enough memory..." << endl;
00186       cerr << "IFPACK WARNING: (file " << __FILE__ << ", line " << __LINE__ 
00187            << ")" << endl;
00188       FirstTime = false;
00189     }
00190     Label_ = "Amesos_Lapack";
00191     Solver_ = Teuchos::rcp( Factory.Create("Amesos_Lapack",*Problem_) );
00192   }
00193   // if empty, give up.
00194   if (Solver_ == Teuchos::null)
00195     IFPACK_CHK_ERR(-1);
00196 
00197   IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_));
00198   Solver_->SetParameters(List_);
00199   IFPACK_CHK_ERR(Solver_->SymbolicFactorization());
00200 
00201   IsInitialized_ = true;
00202   ++NumInitialize_;
00203   InitializeTime_ += Time_->ElapsedTime();
00204   return(0);
00205 }
00206 
00207 //==============================================================================
00208 int Ifpack_Amesos::Compute()
00209 {
00210 
00211   if (!IsInitialized())
00212     IFPACK_CHK_ERR(Initialize());
00213 
00214   if (IsEmpty_) {
00215     IsComputed_ = true;
00216     ++NumCompute_;
00217     return(0);
00218   }
00219 
00220   IsComputed_ = false;
00221   Time_->ResetStartTime();
00222 
00223   if (Matrix_ == Teuchos::null)
00224     IFPACK_CHK_ERR(-1);
00225 
00226   IFPACK_CHK_ERR(Solver_->NumericFactorization());
00227 
00228   IsComputed_ = true;
00229   ++NumCompute_;
00230   ComputeTime_ += Time_->ElapsedTime();
00231   return(0);
00232 }
00233 
00234 //==============================================================================
00235 int Ifpack_Amesos::SetUseTranspose(bool UseTranspose_in)
00236 {
00237   // store the value in UseTranspose_. If we have the solver, we pass to it
00238   // right away, otherwise we wait till when it is created.
00239   UseTranspose_ = UseTranspose_in;
00240   if (Solver_ != Teuchos::null)
00241     IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_in));
00242 
00243   return(0);
00244 }
00245 
00246 //==============================================================================
00247 int Ifpack_Amesos::
00248 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00249 {
00250   // check for maps ???
00251   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00252   return(0);
00253 }
00254 
00255 //==============================================================================
00256 int Ifpack_Amesos::
00257 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00258 {
00259   if (IsEmpty_) {
00260     ++NumApplyInverse_;
00261     return(0);
00262   }
00263 
00264   if (IsComputed() == false)
00265     IFPACK_CHK_ERR(-1);
00266 
00267   if (X.NumVectors() != Y.NumVectors())
00268     IFPACK_CHK_ERR(-1); // wrong input
00269   
00270   Time_->ResetStartTime();
00271 
00272   // AztecOO gives X and Y pointing to the same memory location,
00273   // need to create an auxiliary vector, Xcopy
00274   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00275   if (X.Pointers()[0] == Y.Pointers()[0])
00276     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00277   else
00278     Xcopy = Teuchos::rcp( &X, false );
00279     
00280   Problem_->SetLHS(&Y);
00281   Problem_->SetRHS((Epetra_MultiVector*)Xcopy.get());
00282   IFPACK_CHK_ERR(Solver_->Solve());
00283 
00284   ++NumApplyInverse_;
00285   ApplyInverseTime_ += Time_->ElapsedTime();
00286 
00287   return(0);
00288 }
00289 
00290 //==============================================================================
00291 double Ifpack_Amesos::NormInf() const
00292 {
00293   return(-1.0);
00294 }
00295 
00296 //==============================================================================
00297 const char* Ifpack_Amesos::Label() const
00298 {
00299   return((char*)Label_.c_str());
00300 }
00301 
00302 //==============================================================================
00303 bool Ifpack_Amesos::UseTranspose() const
00304 {
00305   return(UseTranspose_);
00306 }
00307 
00308 //==============================================================================
00309 bool Ifpack_Amesos::HasNormInf() const
00310 {
00311   return(false);
00312 }
00313 
00314 //==============================================================================
00315 const Epetra_Comm & Ifpack_Amesos::Comm() const
00316 {
00317   return(Matrix_->Comm());
00318 }
00319 
00320 //==============================================================================
00321 const Epetra_Map & Ifpack_Amesos::OperatorDomainMap() const
00322 {
00323   return(Matrix_->OperatorDomainMap());
00324 }
00325 
00326 //==============================================================================
00327 const Epetra_Map & Ifpack_Amesos::OperatorRangeMap() const
00328 {
00329   return(Matrix_->OperatorRangeMap());
00330 }
00331 
00332 //==============================================================================
00333 double Ifpack_Amesos::Condest(const Ifpack_CondestType CT,
00334                               const int MaxIters, const double Tol,
00335                   Epetra_RowMatrix* Matrix_in)
00336 {
00337 
00338   if (!IsComputed()) // cannot compute right now
00339     return(-1.0);
00340 
00341   if (Condest_ == -1.0)
00342     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00343 
00344   return(Condest_);
00345 }
00346 
00347 //==============================================================================
00348 std::ostream& Ifpack_Amesos::Print(std::ostream& os) const
00349 {
00350   if (!Comm().MyPID()) {
00351     os << endl;
00352     os << "================================================================================" << endl;
00353     os << "Ifpack_Amesos: " << Label () << endl << endl;
00354     os << "Condition number estimate = " << Condest() << endl;
00355     os << "Global number of rows            = " << Matrix_->NumGlobalRows64() << endl;
00356     os << endl;
00357     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00358     os << "-----           -------   --------------       ------------     --------" << endl;
00359     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00360        << "  " << std::setw(15) << InitializeTime_ 
00361        << "              0.0              0.0" << endl;
00362     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00363        << "  " << std::setw(15) << ComputeTime_
00364        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00365     if (ComputeTime_ != 0.0) 
00366       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00367     else
00368       os << "  " << std::setw(15) << 0.0 << endl;
00369     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00370        << "  " << std::setw(15) << ApplyInverseTime_
00371        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00372     if (ApplyInverseTime_ != 0.0) 
00373       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00374     else
00375       os << "  " << std::setw(15) << 0.0 << endl;
00376     os << "================================================================================" << endl;
00377     os << endl;
00378   }
00379 
00380   return(os);
00381 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends