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     Importer_ = rcp(new Epetra_Import(SerialMap(), Matrix()->RowMatrixRowMap()));
00158     if (Importer_.get() == 0)
00159       AMESOS_CHK_ERR(-1);
00160   }
00161 
00162   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00163   IsSymbolicFactorizationOK_ = true;
00164   ++NumSymbolicFact_;
00165 
00166   return(0);
00167 }
00168 
00169 //=============================================================================
00170 int Amesos_Lapack::NumericFactorization() 
00171 {
00172   IsNumericFactorizationOK_ = false;
00173 
00174   // perform the symbolic (that is, build map and importer) if not done yet.
00175   if (IsSymbolicFactorizationOK_ == false)
00176     AMESOS_CHK_ERR(SymbolicFactorization());
00177 
00178   // Only on processor 0 define the dense matrix.
00179   if (MyPID_ == 0)
00180     AMESOS_CHK_ERR(DenseMatrix_.Shape(NumGlobalRows(),NumGlobalRows()));
00181 
00182   AMESOS_CHK_ERR(DistributedToSerial());
00183   AMESOS_CHK_ERR(SerialToDense());
00184   AMESOS_CHK_ERR(DenseToFactored());
00185 
00186   IsNumericFactorizationOK_ = true;
00187   ++NumNumericFact_;
00188   
00189   return(0);
00190 }
00191 
00192 //=============================================================================
00193 int Amesos_Lapack::Solve() 
00194 {
00195   if (IsNumericFactorizationOK_ == false)
00196     AMESOS_CHK_ERR(NumericFactorization());
00197 
00198   Epetra_MultiVector*       X = Problem_->GetLHS();
00199   const Epetra_MultiVector* B = Problem_->GetRHS();
00200 
00201   if (X == 0) AMESOS_CHK_ERR(-1); // something wrong with user's input
00202   if (B == 0) AMESOS_CHK_ERR(-1); // something wrong with user's input
00203   if (X->NumVectors() != B->NumVectors())
00204     AMESOS_CHK_ERR(-1); // something wrong with user's input
00205 
00206   // Timing are set inside each Solve().
00207   int ierr;
00208   if (NumProcs_ == 1)
00209     ierr = SolveSerial(*X,*B);
00210   else
00211     ierr = SolveDistributed(*X,*B);
00212 
00213   AMESOS_RETURN(ierr);
00214 }
00215 
00216 //=============================================================================
00217 int Amesos_Lapack::SolveSerial(Epetra_MultiVector& X,
00218              const Epetra_MultiVector& B) 
00219 {
00220   ResetTimer();
00221   
00222   int NumVectors = X.NumVectors();
00223 
00224   Epetra_SerialDenseMatrix DenseX(NumGlobalRows(),NumVectors);
00225   Epetra_SerialDenseMatrix DenseB(NumGlobalRows(),NumVectors);
00226 
00227   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00228     for (int j = 0 ; j < NumVectors ; ++j)
00229       DenseB(i,j) = B[j][i];
00230 
00231   DenseSolver_.SetVectors(DenseX,DenseB);
00232   DenseSolver_.SolveWithTranspose(UseTranspose());
00233   AMESOS_CHK_ERR(DenseSolver_.Solve());
00234 
00235   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00236     for (int j = 0 ; j < NumVectors ; ++j)
00237        X[j][i] = DenseX(i,j);
00238 
00239   SolveTime_ = AddTime("Total solve time", SolveTime_);
00240   ++NumSolve_;
00241 
00242   return(0) ;
00243 }
00244 
00245 //=============================================================================
00246 int Amesos_Lapack::SolveDistributed(Epetra_MultiVector& X,
00247             const Epetra_MultiVector& B) 
00248 {
00249   ResetTimer();
00250   
00251   int NumVectors = X.NumVectors();
00252 
00253   // vector that contains RHS, overwritten by solution,
00254   // with all elements on on process 0.
00255   Epetra_MultiVector SerialVector(SerialMap(),NumVectors);
00256   // import off-process data
00257   AMESOS_CHK_ERR(SerialVector.Import(B,Importer(),Insert));
00258 
00259   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00260   ResetTimer();
00261 
00262   if (MyPID_ == 0) {
00263     Epetra_SerialDenseMatrix DenseX(NumGlobalRows(),NumVectors);
00264     Epetra_SerialDenseMatrix DenseB(NumGlobalRows(),NumVectors);
00265 
00266     for (int i = 0 ; i < NumGlobalRows() ; ++i)
00267       for (int j = 0 ; j < NumVectors ; ++j)
00268   DenseB(i,j) = SerialVector[j][i];
00269 
00270     DenseSolver_.SetVectors(DenseX,DenseB);
00271     DenseSolver_.SolveWithTranspose(UseTranspose());
00272     AMESOS_CHK_ERR(DenseSolver_.Solve());
00273 
00274     for (int i = 0 ; i < NumGlobalRows() ; ++i)
00275       for (int j = 0 ; j < NumVectors ; ++j)
00276   SerialVector[j][i] = DenseX(i,j);
00277   }
00278 
00279   SolveTime_ = AddTime("Total solve time", SolveTime_);
00280   ResetTimer();
00281 
00282   AMESOS_CHK_ERR(X.Export(SerialVector,Importer(),Insert));
00283 
00284   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00285   ++NumSolve_;
00286 
00287   return(0) ;
00288 }
00289 
00290 //=============================================================================
00291 int Amesos_Lapack::SerialToDense()
00292 {
00293   if (MyPID_)
00294     return(0);
00295 
00296   ResetTimer();
00297 
00298   for (int i = 0 ; i < NumGlobalRows() ; ++i)
00299     for (int j = 0 ; j < NumGlobalRows() ; ++j)
00300       DenseMatrix_(i,j) = 0.0;
00301 
00302   // allocate storage to extract matrix rows.
00303   int Length = SerialMatrix().MaxNumEntries();
00304   std::vector<double> Values(Length);
00305   std::vector<int>    Indices(Length);
00306 
00307   for (int j = 0 ; j < SerialMatrix().NumMyRows() ; ++j) 
00308   {
00309     int NumEntries;
00310     int ierr = SerialMatrix().ExtractMyRowCopy(j, Length, NumEntries,
00311                                                &Values[0], &Indices[0]);
00312     AMESOS_CHK_ERR(ierr);
00313 
00314     if ( AddZeroToDiag_ ) { 
00315       for (int k = 0 ; k < NumEntries ; ++k) {
00316   DenseMatrix_(j,Indices[k]) = Values[k];
00317       }
00318       DenseMatrix_(j,j) += AddToDiag_;
00319     } else { 
00320       for (int k = 0 ; k < NumEntries ; ++k) {
00321   //      if (fabs(Values[k]) >= Threshold_)       Threshold not used yet - no consistent definition 
00322   //      Lapack would not be the first routine to use a threshold, as it confers no performance
00323   //      advantage
00324   DenseMatrix_(j,Indices[k]) = Values[k];
00325 
00326   //  std::cout << __FILE__ << "::"<<__LINE__ << " j = " << j << " k = " << k << " Values[k = " << Values[k] << " AddToDiag_ = " << AddToDiag_  << std::endl ; 
00327   if (Indices[k] == j) {
00328     DenseMatrix_(j,j) = Values[k] + AddToDiag_;
00329   }
00330       }
00331     }
00332   }
00333 
00334   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
00335   return(0);
00336 }
00337 
00338 //=============================================================================
00339 int Amesos_Lapack::DistributedToSerial()
00340 {
00341   ResetTimer();
00342 
00343   if (NumProcs_ == 1)
00344     SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(Matrix()), false);
00345   else
00346   {
00347     SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy,SerialMap(), 0));
00348     SerialMatrix_    = rcp(&SerialCrsMatrix(), false);
00349     AMESOS_CHK_ERR(SerialCrsMatrix().Import(*Matrix(),Importer(),Insert));
00350     AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
00351   }
00352 
00353   MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
00354 
00355   return(0);
00356 } 
00357 
00358 // ====================================================================== 
00359 int Amesos_Lapack::GEEV(Epetra_Vector& Er, Epetra_Vector& Ei)
00360 {
00361   if (IsSymbolicFactorizationOK_ == false)
00362     AMESOS_CHK_ERR(SymbolicFactorization());
00363 
00364   if (MyPID_ == 0)
00365     AMESOS_CHK_ERR(DenseMatrix_.Shape(NumGlobalRows(),NumGlobalRows()));
00366 
00367   AMESOS_CHK_ERR(DistributedToSerial());
00368   AMESOS_CHK_ERR(SerialToDense());
00369 
00370   Teuchos::RCP<Epetra_Vector> LocalEr;
00371   Teuchos::RCP<Epetra_Vector> LocalEi;
00372 
00373   if (NumProcs_ == 1)
00374   {
00375     LocalEr = Teuchos::rcp(&Er, false);
00376     LocalEi = Teuchos::rcp(&Ei, false);
00377   }
00378   else
00379   {
00380     LocalEr = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
00381     LocalEi = Teuchos::rcp(new Epetra_Vector(*SerialMap_));
00382   }
00383 
00384   if (MyPID_ == 0) 
00385   {
00386     int n = NumGlobalRows();
00387     char jobvl = 'N'; /* V/N to calculate/not calculate left eigenvectors
00388                          of matrix H.*/
00389     char jobvr = 'N'; /* As above, but for right eigenvectors. */
00390     int info = 0;
00391     int ldvl = n;
00392     int ldvr = n;
00393 
00394     Er.PutScalar(0.0);
00395     Ei.PutScalar(0.0);
00396 
00397     Epetra_LAPACK LAPACK;
00398 
00399     std::vector<double> work(1);
00400     int lwork = -1;
00401 
00402     LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
00403                 LocalEr->Values(), LocalEi->Values(), NULL,
00404                 ldvl, NULL, 
00405                 ldvr, &work[0], 
00406                 lwork, &info);
00407 
00408     lwork = (int)work[0];
00409     work.resize(lwork);
00410     LAPACK.GEEV(jobvl, jobvr, n, DenseMatrix_.A(), n, 
00411                 LocalEr->Values(), LocalEi->Values(), NULL,
00412                 ldvl, NULL, 
00413                 ldvr, &work[0], 
00414                 lwork, &info);
00415 
00416     if (info)
00417       AMESOS_CHK_ERR(info);
00418   }
00419 
00420   if (NumProcs_ != 1)
00421   {
00422     // I am not really sure that exporting the results make sense... 
00423     // It is just to be coherent with the other parts of the code.
00424     Er.Export(*LocalEr, Importer(), Insert);
00425     Ei.Export(*LocalEi, Importer(), Insert);
00426   }
00427 
00428   return(0);
00429 }
00430 
00431 //=============================================================================
00432 int Amesos_Lapack::DenseToFactored() 
00433 {
00434   ResetTimer();
00435 
00436   if (MyPID_ == 0) {
00437     AMESOS_CHK_ERR(DenseSolver_.SetMatrix(DenseMatrix_));
00438     int DenseSolverReturn = DenseSolver_.Factor();
00439     if (DenseSolverReturn == 2 ) return NumericallySingularMatrixError ;
00440   }
00441 
00442   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00443   return(0);
00444 }
00445 
00446 // ================================================ ====== ==== ==== == =
00447 void Amesos_Lapack::PrintStatus() const
00448 {
00449   if (MyPID_) return;
00450 
00451   PrintLine();
00452   std::string p = "Amesos_Lapack : ";
00453  
00454   int percentage = 0;
00455   if (NumGlobalRows_ != 0)
00456     percentage = NumGlobalNonzeros_ / (NumGlobalRows_ * NumGlobalRows_);
00457 
00458   std::cout << p << "Matrix has " << NumGlobalRows_ << " rows"
00459        << " and " << NumGlobalNonzeros_ << " nonzeros" << std::endl;
00460   if (NumGlobalRows_ != 0)
00461   {
00462     std::cout << p << "Nonzero elements per row = "
00463          << 1.0 * NumGlobalNonzeros_ / NumGlobalRows_ << std::endl;
00464     std::cout << p << "Percentage of nonzero elements = "
00465          << 100.0 * percentage  << std::endl;
00466   }
00467   std::cout << p << "Use transpose = " << UseTranspose_ << std::endl;
00468 
00469   PrintLine();
00470 
00471   return;
00472 }
00473 
00474 // ================================================ ====== ==== ==== == =
00475 void Amesos_Lapack::PrintTiming() const
00476 {
00477   if (MyPID_ != 0) return;
00478 
00479   double ConTime = GetTime(MtxConvTime_);
00480   double MatTime = GetTime(MtxRedistTime_);
00481   double VecTime = GetTime(VecRedistTime_);
00482   double SymTime = GetTime(SymFactTime_);
00483   double NumTime = GetTime(NumFactTime_);
00484   double SolTime = GetTime(SolveTime_);
00485 
00486   if (NumSymbolicFact_)
00487     SymTime /= NumSymbolicFact_;
00488 
00489   if (NumNumericFact_)
00490     NumTime /= NumNumericFact_;
00491 
00492   if (NumSolve_)
00493     SolTime /= NumSolve_;
00494 
00495   std::string p = "Amesos_Lapack : ";
00496   PrintLine();
00497 
00498   std::cout << p << "Time to convert matrix to Klu format = "
00499        << ConTime << " (s)" << std::endl;
00500   std::cout << p << "Time to redistribute matrix = "
00501        << MatTime << " (s)" << std::endl;
00502   std::cout << p << "Time to redistribute vectors = "
00503        << VecTime << " (s)" << std::endl;
00504   std::cout << p << "Number of symbolic factorizations = "
00505        << NumSymbolicFact_ << std::endl;
00506   std::cout << p << "Time for sym fact = "
00507        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00508   std::cout << p << "Number of numeric factorizations = "
00509        << NumNumericFact_ << std::endl;
00510   std::cout << p << "Time for num fact = "
00511        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00512   std::cout << p << "Number of solve phases = "
00513        << NumSolve_ << std::endl;
00514   std::cout << p << "Time for solve = "
00515        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00516 
00517   PrintLine();
00518 
00519   return;
00520 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines