Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Lapack.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) 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 #include "Amesos_Lapack.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Import.h"
00032 #include "Epetra_RowMatrix.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_Util.h"
00035 #include "Epetra_LAPACK.h"
00036 using namespace Teuchos;
00037 
00038 //=============================================================================
00039 Amesos_Lapack::Amesos_Lapack(const Epetra_LinearProblem &Problem) :
00040   UseTranspose_(false),
00041   Problem_(&Problem),
00042   MtxRedistTime_(-1), 
00043   MtxConvTime_(-1), 
00044   VecRedistTime_(-1),
00045   SymFactTime_(-1), 
00046   NumFactTime_(-1), 
00047   SolveTime_(-1)
00048 
00049 {
00050   ParameterList_ = unsetParameterList();
00051 }
00052 
00053 Teuchos::RCP<Teuchos::ParameterList> Amesos_Lapack::unsetParameterList() {
00054   Teuchos::RCP<Teuchos::ParameterList> PL = rcp( new Teuchos::ParameterList);
00055   setParameterList(PL);
00056   return PL ; 
00057 }
00058 
00059 
00060 void Amesos_Lapack::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& pl) {
00061 
00062   pl_ = pl; 
00063   Teuchos::ParameterList& LAPACKParams = pl_->sublist("Lapack") ;
00064   AddZeroToDiag_ = LAPACKParams.get( "AddZeroToDiag", false );
00065   AddToDiag_ = LAPACKParams.get( "AddToDiag", 0.0 );
00066 
00067   bool Equilibrate = LAPACKParams.get("Equilibrate",true);
00068   DenseSolver_.FactorWithEquilibration(Equilibrate);
00069 }
00070 
00071 
00072 
00073 //=============================================================================
00074 Amesos_Lapack::~Amesos_Lapack(void) 
00075 {
00076   // print out some information if required by the user
00077   if ((verbose_ && PrintTiming_) || (verbose_ == 2)) 
00078     PrintTiming();
00079   if ((verbose_ && PrintStatus_) || (verbose_ == 2)) 
00080     PrintStatus();
00081 }
00082 
00083 //=============================================================================
00084 // Default values are defined in the constructor.
00085 int Amesos_Lapack::SetParameters( Teuchos::ParameterList &ParameterList ) 
00086 {
00087   // retrive general parameters
00088   SetStatusParameters( ParameterList );
00089 
00090   SetControlParameters( ParameterList );
00091 
00092   bool Equilibrate = true;
00093   if (ParameterList.isSublist("Lapack") ) {
00094     const Teuchos::ParameterList& LAPACKParams = ParameterList.sublist("Lapack") ;
00095     if ( LAPACKParams.isParameter("Equilibrate") )
00096       Equilibrate = LAPACKParams.get<bool>("Equilibrate");
00097   }
00098   DenseSolver_.FactorWithEquilibration(Equilibrate);
00099 
00100   return(0);
00101 }
00102 
00103 //=============================================================================
00104 bool Amesos_Lapack::MatrixShapeOK() const 
00105 {
00106   if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
00107       GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints()) {
00108     return(false);
00109   }
00110   else
00111     return(true);
00112 }
00113 
00114 //=============================================================================
00115 // Build SerialMap_ and Importer_. SerialMap_ is equal to RowMatrixRowMap()
00116 // for the serial case, otherwise we actually need to build a new map.
00117 // Importer_ is constructed only for the parallel case. All these objects are
00118 // destroyed and rebuilt every time SymbolicFactorization() is called, to
00119 // allow the possibility of completely different matrices given in input.
00120 //
00121 // NOTE: A possible bug is nasted here. If the user destroys the
00122 // RowMatrixRowMap(), Importer() may not work properly any more. Anyway,
00123 // this bug is due to misuse of the class, since it is supposed that the
00124 // structure of A (and therefore all its components) remain untouched after a
00125 // call to SymbolicFactorization().
00126 // ============================================================================
00127 int Amesos_Lapack::SymbolicFactorization() 
00128 {
00129   IsSymbolicFactorizationOK_ = false;
00130   IsNumericFactorizationOK_ = false;
00131 
00132   if (GetProblem()->GetMatrix() == 0)
00133     AMESOS_CHK_ERR(-1); // problem still not properly set
00134 
00135   CreateTimer(Comm()); // Create timer
00136   ResetTimer();
00137 
00138   MyPID_             = Comm().MyPID();
00139   NumProcs_          = Comm().NumProc();
00140   NumGlobalRows_     = Matrix()->NumGlobalRows();
00141   NumGlobalNonzeros_ = Matrix()->NumGlobalNonzeros();
00142 
00143 
00144   if (NumProcs_ == 1)
00145     SerialMap_ = rcp(const_cast<Epetra_Map*>(&(Matrix()->RowMatrixRowMap())), 
00146                      false);
00147   else
00148   {
00149     int NumElements = 0;
00150     if (MyPID_ == 0)
00151       NumElements = NumGlobalRows();
00152 
00153     SerialMap_ = rcp(new Epetra_Map(-1, NumElements, 0, Comm()));
00154     if (SerialMap_.get() == 0)
00155       AMESOS_CHK_ERR(-1);
00156 
00157     MatrixImporter_ = rcp(new Epetra_Import(SerialMap(), Matrix()->RowMatrixRowMap()));
00158     if (MatrixImporter_.get() == 0)
00159       AMESOS_CHK_ERR(-1);
00160 
00161     const bool switchDomainRangeMaps = (UseTranspose_ != Matrix()->UseTranspose());
00162 
00163     const Epetra_Map &rhsMap = switchDomainRangeMaps ? Matrix()->OperatorDomainMap() : Matrix()->OperatorRangeMap();
00164     RhsExporter_ = rcp(new Epetra_Export(rhsMap, SerialMap()));
00165     if (RhsExporter_.get() == 0)
00166       AMESOS_CHK_ERR(-1);
00167 
00168     const Epetra_Map &solutionMap = switchDomainRangeMaps ? Matrix()->OperatorRangeMap() : Matrix()->OperatorDomainMap();
00169     SolutionImporter_ = rcp(new Epetra_Import(solutionMap, SerialMap()));
00170     if (SolutionImporter_.get() == 0)
00171       AMESOS_CHK_ERR(-1);
00172   }
00173 
00174   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00175   IsSymbolicFactorizationOK_ = true;
00176   ++NumSymbolicFact_;
00177 
00178   return(0);
00179 }
00180 
00181 //=============================================================================
00182 int Amesos_Lapack::NumericFactorization() 
00183 {
00184   IsNumericFactorizationOK_ = false;
00185 
00186   // perform the symbolic (that is, build map and importer) if not done yet.
00187   if (IsSymbolicFactorizationOK_ == false)
00188     AMESOS_CHK_ERR(SymbolicFactorization());
00189 
00190   // Only on processor 0 define the dense matrix.
00191   if (MyPID_ == 0)
00192     AMESOS_CHK_ERR(DenseMatrix_.Shape(NumGlobalRows(),NumGlobalRows()));
00193 
00194   AMESOS_CHK_ERR(DistributedToSerial());
00195   AMESOS_CHK_ERR(SerialToDense());
00196   AMESOS_CHK_ERR(DenseToFactored());
00197 
00198   IsNumericFactorizationOK_ = true;
00199   ++NumNumericFact_;
00200   
00201   return(0);
00202 }
00203 
00204 //=============================================================================
00205 int Amesos_Lapack::Solve() 
00206 {
00207   if (IsNumericFactorizationOK_ == false)
00208     AMESOS_CHK_ERR(NumericFactorization());
00209 
00210   Epetra_MultiVector*       X = Problem_->GetLHS();
00211   const Epetra_MultiVector* B = Problem_->GetRHS();
00212 
00213   if (X == 0) AMESOS_CHK_ERR(-1); // something wrong with user's input
00214   if (B == 0) AMESOS_CHK_ERR(-1); // something wrong with user's input
00215   if (X->NumVectors() != B->NumVectors())
00216     AMESOS_CHK_ERR(-1); // something wrong with user's input
00217 
00218   // Timing are set inside each Solve().
00219   int ierr;
00220   if (NumProcs_ == 1)
00221     ierr = SolveSerial(*X,*B);
00222   else
00223     ierr = SolveDistributed(*X,*B);
00224 
00225   AMESOS_RETURN(ierr);
00226 }
00227 
00228 //=============================================================================
00229 int Amesos_Lapack::SolveSerial(Epetra_MultiVector& X,
00230              const Epetra_MultiVector& B) 
00231 {
00232   ResetTimer();
00233   
00234   int NumVectors = X.NumVectors();
00235 
00236   Epetra_SerialDenseMatrix DenseX(NumGlobalRows(),NumVectors);
00237   Epetra_SerialDenseMatrix DenseB(NumGlobalRows(),NumVectors);
00238 
00239   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00240     for (int j = 0 ; j < NumVectors ; ++j)
00241       DenseB(i,j) = B[j][i];
00242 
00243   DenseSolver_.SetVectors(DenseX,DenseB);
00244   DenseSolver_.SolveWithTranspose(UseTranspose());
00245   AMESOS_CHK_ERR(DenseSolver_.Solve());
00246 
00247   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00248     for (int j = 0 ; j < NumVectors ; ++j)
00249        X[j][i] = DenseX(i,j);
00250 
00251   SolveTime_ = AddTime("Total solve time", SolveTime_);
00252   ++NumSolve_;
00253 
00254   return(0) ;
00255 }
00256 
00257 //=============================================================================
00258 int Amesos_Lapack::SolveDistributed(Epetra_MultiVector& X,
00259             const Epetra_MultiVector& B) 
00260 {
00261   ResetTimer();
00262   
00263   int NumVectors = X.NumVectors();
00264 
00265   // vector that contains RHS, overwritten by solution,
00266   // with all elements on on process 0.
00267   Epetra_MultiVector SerialVector(SerialMap(),NumVectors);
00268   // import off-process data
00269   AMESOS_CHK_ERR(SerialVector.Export(B,RhsExporter(),Insert));
00270 
00271   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00272   ResetTimer();
00273 
00274   if (MyPID_ == 0) {
00275     Epetra_SerialDenseMatrix DenseX(NumGlobalRows(),NumVectors);
00276     Epetra_SerialDenseMatrix DenseB(NumGlobalRows(),NumVectors);
00277 
00278     for (int i = 0 ; i < NumGlobalRows() ; ++i)
00279       for (int j = 0 ; j < NumVectors ; ++j)
00280   DenseB(i,j) = SerialVector[j][i];
00281 
00282     DenseSolver_.SetVectors(DenseX,DenseB);
00283     DenseSolver_.SolveWithTranspose(UseTranspose());
00284     AMESOS_CHK_ERR(DenseSolver_.Solve());
00285 
00286     for (int i = 0 ; i < NumGlobalRows() ; ++i)
00287       for (int j = 0 ; j < NumVectors ; ++j)
00288   SerialVector[j][i] = DenseX(i,j);
00289   }
00290 
00291   SolveTime_ = AddTime("Total solve time", SolveTime_);
00292   ResetTimer();
00293 
00294   AMESOS_CHK_ERR(X.Import(SerialVector,SolutionImporter(),Insert));
00295 
00296   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00297   ++NumSolve_;
00298 
00299   return(0) ;
00300 }
00301 
00302 //=============================================================================
00303 int Amesos_Lapack::SerialToDense()
00304 {
00305   if (MyPID_)
00306     return(0);
00307 
00308   ResetTimer();
00309 
00310   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00311     for (int j = 0 ; j < NumGlobalRows() ; ++j)
00312       DenseMatrix_(i,j) = 0.0;
00313 
00314   // allocate storage to extract matrix rows.
00315   int Length = SerialMatrix().MaxNumEntries();
00316   std::vector<double> Values(Length);
00317   std::vector<int>    Indices(Length);
00318 
00319   for (int j = 0 ; j < SerialMatrix().NumMyRows() ; ++j) 
00320   {
00321     int NumEntries;
00322     int ierr = SerialMatrix().ExtractMyRowCopy(j, Length, NumEntries,
00323                                                &Values[0], &Indices[0]);
00324     AMESOS_CHK_ERR(ierr);
00325 
00326     if ( AddZeroToDiag_ ) { 
00327       for (int k = 0 ; k < NumEntries ; ++k) {
00328   DenseMatrix_(j,Indices[k]) = Values[k];
00329       }
00330       DenseMatrix_(j,j) += AddToDiag_;
00331     } else { 
00332       for (int k = 0 ; k < NumEntries ; ++k) {
00333   //      if (fabs(Values[k]) >= Threshold_)       Threshold not used yet - no consistent definition 
00334   //      Lapack would not be the first routine to use a threshold, as it confers no performance
00335   //      advantage
00336   DenseMatrix_(j,Indices[k]) = Values[k];
00337 
00338   //  std::cout << __FILE__ << "::"<<__LINE__ << " j = " << j << " k = " << k << " Values[k = " << Values[k] << " AddToDiag_ = " << AddToDiag_  << std::endl ; 
00339   if (Indices[k] == j) {
00340     DenseMatrix_(j,j) = Values[k] + AddToDiag_;
00341   }
00342       }
00343     }
00344   }
00345 
00346   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
00347   return(0);
00348 }
00349 
00350 //=============================================================================
00351 int Amesos_Lapack::DistributedToSerial()
00352 {
00353   ResetTimer();
00354 
00355   if (NumProcs_ == 1)
00356     SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(Matrix()), false);
00357   else
00358   {
00359     SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy,SerialMap(), 0));
00360     SerialMatrix_    = rcp(&SerialCrsMatrix(), false);
00361     AMESOS_CHK_ERR(SerialCrsMatrix().Import(*Matrix(),MatrixImporter(),Insert));
00362     AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
00363   }
00364 
00365   MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
00366 
00367   return(0);
00368 } 
00369 
00370 // ====================================================================== 
00371 int Amesos_Lapack::GEEV(Epetra_Vector& Er, Epetra_Vector& Ei)
00372 {
00373   if (IsSymbolicFactorizationOK_ == false)
00374     AMESOS_CHK_ERR(SymbolicFactorization());
00375 
00376   if (MyPID_ == 0)
00377     AMESOS_CHK_ERR(DenseMatrix_.Shape(NumGlobalRows(),NumGlobalRows()));
00378 
00379   AMESOS_CHK_ERR(DistributedToSerial());
00380   AMESOS_CHK_ERR(SerialToDense());
00381 
00382   Teuchos::RCP<Epetra_Vector> LocalEr;
00383   Teuchos::RCP<Epetra_Vector> LocalEi;
00384 
00385   if (NumProcs_ == 1)
00386   {
00387     LocalEr = Teuchos::rcp(&Er, false);
00388     LocalEi = Teuchos::rcp(&Ei, false);
00389   }
00390   else
00391   {
00392     LocalEr = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
00393     LocalEi = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
00394   }
00395 
00396   if (MyPID_ == 0) 
00397   {
00398     int n = NumGlobalRows();
00399     char jobvl = 'N'; /* V/N to calculate/not calculate left eigenvectors
00400                          of matrix H.*/
00401     char jobvr = 'N'; /* As above, but for right eigenvectors. */
00402     int info = 0;
00403     int ldvl = n;
00404     int ldvr = n;
00405 
00406     Er.PutScalar(0.0);
00407     Ei.PutScalar(0.0);
00408 
00409     Epetra_LAPACK LAPACK;
00410 
00411     std::vector<double> work(1);
00412     int lwork = -1;
00413 
00414     LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
00415                 LocalEr->Values(), LocalEi->Values(), NULL,
00416                 ldvl, NULL, 
00417                 ldvr, &work[0], 
00418                 lwork, &info);
00419 
00420     lwork = (int)work[0];
00421     work.resize(lwork);
00422     LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
00423                 LocalEr->Values(), LocalEi->Values(), NULL,
00424                 ldvl, NULL, 
00425                 ldvr, &work[0], 
00426                 lwork, &info);
00427 
00428     if (info)
00429       AMESOS_CHK_ERR(info);
00430   }
00431 
00432   if (NumProcs_ != 1)
00433   {
00434     // I am not really sure that exporting the results make sense... 
00435     // It is just to be coherent with the other parts of the code.
00436     Er.Import(*LocalEr, Epetra_Import(Er.Map(), SerialMap()), Insert);
00437     Ei.Import(*LocalEi, Epetra_Import(Ei.Map(), SerialMap()), Insert);
00438   }
00439 
00440   return(0);
00441 }
00442 
00443 //=============================================================================
00444 int Amesos_Lapack::DenseToFactored() 
00445 {
00446   ResetTimer();
00447 
00448   if (MyPID_ == 0) {
00449     AMESOS_CHK_ERR(DenseSolver_.SetMatrix(DenseMatrix_));
00450     int DenseSolverReturn = DenseSolver_.Factor();
00451     if (DenseSolverReturn == 2 ) return NumericallySingularMatrixError ;
00452   }
00453 
00454   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00455   return(0);
00456 }
00457 
00458 // ================================================ ====== ==== ==== == =
00459 void Amesos_Lapack::PrintStatus() const
00460 {
00461   if (MyPID_) return;
00462 
00463   PrintLine();
00464   std::string p = "Amesos_Lapack : ";
00465  
00466   int percentage = 0;
00467   if (NumGlobalRows_ != 0)
00468     percentage = NumGlobalNonzeros_ / (NumGlobalRows_ * NumGlobalRows_);
00469 
00470   std::cout << p << "Matrix has " << NumGlobalRows_ << " rows"
00471        << " and " << NumGlobalNonzeros_ << " nonzeros" << std::endl;
00472   if (NumGlobalRows_ != 0)
00473   {
00474     std::cout << p << "Nonzero elements per row = "
00475          << 1.0 * NumGlobalNonzeros_ / NumGlobalRows_ << std::endl;
00476     std::cout << p << "Percentage of nonzero elements = "
00477          << 100.0 * percentage  << std::endl;
00478   }
00479   std::cout << p << "Use transpose = " << UseTranspose_ << std::endl;
00480 
00481   PrintLine();
00482 
00483   return;
00484 }
00485 
00486 // ================================================ ====== ==== ==== == =
00487 void Amesos_Lapack::PrintTiming() const
00488 {
00489   if (MyPID_ != 0) return;
00490 
00491   double ConTime = GetTime(MtxConvTime_);
00492   double MatTime = GetTime(MtxRedistTime_);
00493   double VecTime = GetTime(VecRedistTime_);
00494   double SymTime = GetTime(SymFactTime_);
00495   double NumTime = GetTime(NumFactTime_);
00496   double SolTime = GetTime(SolveTime_);
00497 
00498   if (NumSymbolicFact_)
00499     SymTime /= NumSymbolicFact_;
00500 
00501   if (NumNumericFact_)
00502     NumTime /= NumNumericFact_;
00503 
00504   if (NumSolve_)
00505     SolTime /= NumSolve_;
00506 
00507   std::string p = "Amesos_Lapack : ";
00508   PrintLine();
00509 
00510   std::cout << p << "Time to convert matrix to Klu format = "
00511        << ConTime << " (s)" << std::endl;
00512   std::cout << p << "Time to redistribute matrix = "
00513        << MatTime << " (s)" << std::endl;
00514   std::cout << p << "Time to redistribute vectors = "
00515        << VecTime << " (s)" << std::endl;
00516   std::cout << p << "Number of symbolic factorizations = "
00517        << NumSymbolicFact_ << std::endl;
00518   std::cout << p << "Time for sym fact = "
00519        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00520   std::cout << p << "Number of numeric factorizations = "
00521        << NumNumericFact_ << std::endl;
00522   std::cout << p << "Time for num fact = "
00523        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00524   std::cout << p << "Number of solve phases = "
00525        << NumSolve_ << std::endl;
00526   std::cout << p << "Time for solve = "
00527        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00528 
00529   PrintLine();
00530 
00531   return;
00532 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines