Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Taucs.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_Taucs.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Import.h"
00032 #include "Epetra_CrsMatrix.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_Util.h"
00035 extern "C" {
00036 #include "taucs.h"
00037 }
00038 
00039 using namespace Teuchos;
00040 
00041 
00042 //
00043 //  deaalocFunctorHandleDelete requires a fucntion that takes a ** argument, these
00044 //  functions meet that requirement
00045 // 
00046 void taucs_dccs_free_ptr( taucs_ccs_matrix** taucs_ccs_matrix_in ) { 
00047   taucs_dccs_free( *taucs_ccs_matrix_in ); 
00048 }
00049 
00050 void taucs_supernodal_factor_free_ptr( taucs_ccs_matrix** taucs_ccs_matrix_in ) { 
00051   taucs_supernodal_factor_free( *taucs_ccs_matrix_in ); 
00052 }
00053 
00054 class Amesos_Taucs_Pimpl {
00055 public:
00056 
00057 #define USE_REF_COUNT_PTR_FOR_A_AND_L
00058 #ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
00059    Teuchos::RCP<taucs_ccs_matrix> A_ ;
00060    Teuchos::RCP<taucs_ccs_matrix> L_ ;
00061 #else
00062    taucs_ccs_matrix* A_ ;
00063    taucs_ccs_matrix* L_ ;
00064 
00065   Amesos_Taucs_Pimpl():
00066     A_(0),
00067     L_(0)
00068   {}
00069 
00070   ~Amesos_Taucs_Pimpl(void){
00071   if (A_ != 0)
00072   {
00073     taucs_dccs_free(A_);
00074     A_ = 0;
00075   }
00076 
00077   if (L_ != 0)
00078   {
00079     taucs_supernodal_factor_free(L_);
00080     L_ = 0;
00081   }
00082   }
00083 
00084 #endif
00085 } ;
00086 
00087 //=============================================================================
00088 Amesos_Taucs::Amesos_Taucs(const Epetra_LinearProblem &prob) :
00089   PrivateTaucsData_( rcp( new Amesos_Taucs_Pimpl() ) ),
00090   Matrix_(0),
00091   Problem_(&prob),
00092   MtxConvTime_(-1),
00093   MtxRedistTime_(-1),
00094   VecRedistTime_(-1),
00095   SymFactTime_(-1),
00096   NumFactTime_(-1),
00097   SolveTime_(-1)
00098 { }
00099 
00100 //=============================================================================
00101 Amesos_Taucs::~Amesos_Taucs(void) 
00102 {
00103 
00104   // print out some information if required by the user
00105   if ((verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
00106   if ((verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
00107 
00108 #if 0
00109   if (A_ != 0)
00110   {
00111     taucs_dccs_free(A_);
00112     A_ = 0;
00113   }
00114 
00115   if (L_ != 0)
00116   {
00117     taucs_supernodal_factor_free(L_);
00118     L_ = 0;
00119   }
00120 #endif
00121 }
00122 
00123 //=============================================================================
00124 int Amesos_Taucs::ConvertToSerial() 
00125 {
00126   ResetTimer();
00127 
00128   int NumGlobalRows = Matrix_->NumGlobalRows();
00129 
00130   // create a serial map
00131   int NumMyRows = 0;
00132   if (Comm().MyPID() == 0) 
00133     NumMyRows = NumGlobalRows;
00134 
00135   SerialMap_ = rcp(new Epetra_Map(-1, NumMyRows, 0, Comm()));
00136   if (SerialMap_.get() == 0)
00137     AMESOS_CHK_ERR(-1);
00138 
00139   Importer_ = rcp(new Epetra_Import(SerialMap(),Map()));
00140   if (Importer_.get() == 0)
00141     AMESOS_CHK_ERR(-1);
00142 
00143   SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy, SerialMap(), 0));
00144   if (SerialCrsMatrix_.get() == 0)
00145     AMESOS_CHK_ERR(-1);
00146 
00147   AMESOS_CHK_ERR(SerialCrsMatrix().Import(Matrix(), Importer(), Add));
00148 
00149   AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
00150 
00151   SerialMatrix_ = rcp(SerialCrsMatrix_.get(), false);
00152 
00153   MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
00154 
00155   return 0;
00156 }
00157 
00158 //=============================================================================
00159 int Amesos_Taucs::ConvertToTaucs()
00160 {
00161   ResetTimer();
00162 
00163   if (Comm().MyPID() == 0) 
00164   {
00165     int n = SerialMatrix().NumMyRows();
00166     int nnz = SerialMatrix().NumMyNonzeros();
00167     int nnz_sym = (nnz) / 2 + n;
00168 
00169 #ifndef USE_REF_COUNT_PTR_FOR_A_AND_L
00170     if (PrivateTaucsData_->A_ != 0)
00171     {
00172       taucs_dccs_free(PrivateTaucsData_->A_);
00173       PrivateTaucsData_->A_ = 0;
00174     }
00175 #endif
00176 
00177 #ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
00178     PrivateTaucsData_->A_ = 
00179       rcp( taucs_ccs_create(n, n, nnz_sym, TAUCS_DOUBLE),
00180      deallocFunctorHandleDelete<taucs_ccs_matrix>(taucs_dccs_free_ptr), true );    
00181 #else
00182     PrivateTaucsData_->A_ = taucs_ccs_create(n, n, nnz_sym, TAUCS_DOUBLE) ;
00183 #endif
00184     PrivateTaucsData_->A_->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
00185 
00186     int count = 0;
00187     int MaxNumEntries = SerialMatrix().MaxNumEntries();
00188     std::vector<int>    Indices(MaxNumEntries);
00189     std::vector<double> Values(MaxNumEntries);
00190 
00191     PrivateTaucsData_->A_->colptr[0] = 0;
00192 
00193     for (int i = 0 ; i < n ; ++i)
00194     {
00195       int count2 = 0;
00196       int ierr, NumEntriesThisRow;
00197       ierr = SerialMatrix().ExtractMyRowCopy(i, MaxNumEntries,
00198                                              NumEntriesThisRow,
00199                                              &Values[0], &Indices[0]);
00200       if (ierr < 0)
00201         AMESOS_CHK_ERR(ierr);
00202 
00203       for (int j = 0 ; j < NumEntriesThisRow ; ++j)
00204       {
00205         if (Indices[j] < i)
00206           continue;
00207 
00208         if (Indices[j] == i)
00209           Values[j] += AddToDiag_;
00210 
00211         PrivateTaucsData_->A_->rowind[count] = Indices[j];
00212         PrivateTaucsData_->A_->values.d[count] = Values[j];
00213         ++count;
00214         ++count2;
00215       }
00216 
00217       PrivateTaucsData_->A_->colptr[i + 1] = PrivateTaucsData_->A_->colptr[i] + count2;
00218     }
00219 
00220     if (count > SerialMatrix().NumMyNonzeros())
00221       AMESOS_CHK_ERR(-1); // something wrong here
00222   }
00223 
00224   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
00225 
00226   return 0;
00227 }
00228 
00229 //=============================================================================
00230 int Amesos_Taucs::SetParameters( Teuchos::ParameterList &ParameterList) 
00231 {
00232   // retrive general parameters
00233 
00234   SetStatusParameters( ParameterList );
00235 
00236   SetControlParameters( ParameterList );
00237 
00238   return 0;
00239 }
00240 
00241 //=============================================================================
00242 int Amesos_Taucs::PerformSymbolicFactorization() 
00243 {
00244   ResetTimer();
00245 
00246   if (Comm().MyPID() == 0)
00247   {
00248 
00249 #ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
00250     PrivateTaucsData_->L_ = 
00251       rcp( ((taucs_ccs_matrix*)taucs_ccs_factor_llt_symbolic(&*PrivateTaucsData_->A_)),
00252      deallocFunctorHandleDelete<taucs_ccs_matrix>(taucs_supernodal_factor_free_ptr), true );    
00253 #else
00254     if (PrivateTaucsData_->L_ != 0)
00255     {
00256       taucs_supernodal_factor_free(PrivateTaucsData_->L_);
00257       PrivateTaucsData_->L_ = 0;
00258     }
00259     PrivateTaucsData_->L_ = 
00260       (taucs_ccs_matrix*)taucs_ccs_factor_llt_symbolic(&*PrivateTaucsData_->A_) ;
00261 #endif 
00262   }
00263 
00264   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00265 
00266   return 0;
00267 }
00268 
00269 //=============================================================================
00270 int Amesos_Taucs::PerformNumericFactorization( ) 
00271 {
00272   ResetTimer();
00273 
00274   if (Comm().MyPID() == 0) 
00275   {
00276     taucs_supernodal_factor_free_numeric( &*PrivateTaucsData_->L_);
00277 
00278     int ierr = taucs_ccs_factor_llt_numeric(&*PrivateTaucsData_->A_, &*PrivateTaucsData_->L_);
00279 
00280     if (ierr != 0) 
00281     {
00282       std::cerr << "Amesos_Taucs: error during numeric factorization ("
00283         << ierr << ")" << std::endl;
00284       AMESOS_CHK_ERR(-1);
00285     }
00286   }
00287 
00288   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00289 
00290   return 0;
00291 }
00292 
00293 //=============================================================================
00294 bool Amesos_Taucs::MatrixShapeOK() const 
00295 {
00296   bool OK = true;
00297 
00298   if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
00299       GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) 
00300   {
00301     OK = false;
00302   }
00303   return OK;
00304 }
00305 
00306 //=============================================================================
00307 int Amesos_Taucs::SymbolicFactorization() 
00308 {
00309   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__ << " Entering SymbolicFactorization()" << std::endl ; 
00310   IsSymbolicFactorizationOK_ = false;
00311   IsNumericFactorizationOK_ = false;
00312 
00313   CreateTimer(Comm());
00314 
00315   ++NumSymbolicFact_;
00316 
00317   Matrix_ = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00318   Map_ = &(Matrix_->RowMatrixRowMap());
00319 
00320   // =========================================================== //
00321   // redistribute and create all import/export objects only      //
00322   // if more than one processor is used. Otherwise simply set    //
00323   // dummy pointers to Matrix() and Map(), without giving the    //
00324   // ownership to the smart pointer.                             //
00325   // =========================================================== //
00326 
00327   if (Comm().NumProc() != 1) 
00328     ConvertToSerial();
00329   else
00330   {
00331     SerialMap_ = rcp(const_cast<Epetra_Map*>(&Map()), false);
00332     SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(&Matrix()), false);
00333   }
00334 
00335   // =========================================================== //
00336   // Only on processor zero, convert the matrix into CSR format, //
00337   // as required by TAUCS.                                       //
00338   // =========================================================== //
00339 
00340   ConvertToTaucs();
00341 
00342   PerformSymbolicFactorization();
00343 
00344   IsSymbolicFactorizationOK_ = true;
00345 
00346   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__  << " Leaving SymbolicFactorization()" << std::endl ; 
00347   return(0);
00348 }
00349 
00350 //=============================================================================
00351 int Amesos_Taucs::NumericFactorization() 
00352 {
00353   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__ << " Entering NumericFactorization()" << std::endl ; 
00354   IsNumericFactorizationOK_ = false;
00355 
00356   if (IsSymbolicFactorizationOK_ == false)
00357     AMESOS_CHK_ERR(SymbolicFactorization());
00358 
00359   ++NumNumericFact_;
00360 
00361   // FIXME: this must be checked, now all the matrix is shipped twice here
00362   ConvertToSerial();
00363   ConvertToTaucs();
00364 
00365   PerformNumericFactorization();
00366 
00367   IsNumericFactorizationOK_ = true;
00368 
00369   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
00370          << " Leaving NumericFactorization()" << std::endl ; 
00371   return(0);
00372 }
00373 
00374 //=============================================================================
00375 int Amesos_Taucs::Solve() 
00376 {
00377   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
00378          << " Entering Solve()" << std::endl ; 
00379   if (IsNumericFactorizationOK_ == false)
00380     AMESOS_CHK_ERR(NumericFactorization());
00381 
00382   ++NumSolve_;
00383 
00384   Epetra_MultiVector* X = Problem_->GetLHS();
00385   Epetra_MultiVector* B = Problem_->GetRHS();
00386 
00387   if ((X == 0) || (B == 0))
00388     AMESOS_CHK_ERR(-1); 
00389 
00390   int NumVectors = X->NumVectors();
00391   if (NumVectors != B->NumVectors())
00392     AMESOS_CHK_ERR(-1); 
00393 
00394   // vectors with SerialMap_
00395   RCP<Epetra_MultiVector> SerialB;
00396   RCP<Epetra_MultiVector> SerialX;
00397 
00398   ResetTimer();
00399 
00400   if (Comm().NumProc() == 1) 
00401   {
00402     SerialB = rcp(B,false);
00403     SerialX = rcp(X,false);
00404   } 
00405   else 
00406   {
00407     SerialX = rcp(new Epetra_MultiVector(SerialMap(),NumVectors));
00408     SerialB = rcp(new Epetra_MultiVector(SerialMap(),NumVectors));
00409 
00410     SerialB->Import(*B,Importer(),Insert);
00411   }
00412 
00413   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00414 
00415   ResetTimer();
00416 
00417   if (Comm().MyPID() == 0) 
00418   {
00419     double* SerialXValues;
00420     double* SerialBValues;
00421     int LDA;
00422 
00423     AMESOS_CHK_ERR(SerialX->ExtractView(&SerialXValues,&LDA));
00424 
00425     // FIXME: check LDA
00426     AMESOS_CHK_ERR(SerialB->ExtractView(&SerialBValues,&LDA));
00427 
00428     for (int i = 0 ; i < NumVectors ; ++i)
00429     {
00430       int ierr = taucs_supernodal_solve_llt(&*PrivateTaucsData_->L_, 
00431                                             SerialXValues + i * LDA,
00432                                             SerialBValues + i * LDA);
00433 
00434       if (ierr != TAUCS_SUCCESS)
00435       {
00436         std::cerr << "Error occurred in taucs_ccs_solve()" << std::endl;
00437         AMESOS_CHK_ERR(-1);
00438       }
00439     }
00440   }
00441 
00442   SolveTime_ = AddTime("Total solve time", SolveTime_);
00443 
00444   ResetTimer();
00445 
00446   if (Comm().NumProc() != 1) 
00447     X->Export(*SerialX, Importer(), Insert);
00448 
00449   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00450 
00451   if (ComputeTrueResidual_)
00452     ComputeTrueResidual(Matrix(), *X, *B, false, "Amesos_Taucs");
00453 
00454   if (ComputeVectorNorms_)
00455     ComputeVectorNorms(*X, *B, "Amesos_Taucs");
00456 
00457   if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
00458          << " Leaving Solve()" << std::endl ; 
00459   return(0) ;
00460 }
00461 
00462 // ====================================================================== 
00463 void Amesos_Taucs::PrintStatus() const
00464 {
00465   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00466     return;
00467 
00468   std::string p = "Amesos_Taucs : ";
00469 
00470   if (Matrix_ == 0) return;
00471 
00472   
00473   PrintLine();
00474   int n = Matrix().NumGlobalRows();
00475   int nnz = Matrix().NumGlobalNonzeros();
00476     
00477   std::cout << p << "Matrix has " << n << " rows"
00478        << " and " << nnz << " nonzeros" << std::endl;
00479   if (n > 0) 
00480   { 
00481     std::cout << p << "Nonzero elements per row = "
00482          << 1.0 *  nnz / n << std::endl;
00483     std::cout << p << "Percentage of nonzero elements = "
00484          << 100.0 * nnz /(pow(n,2.0)) << std::endl;
00485   }
00486   PrintLine();
00487 
00488   return;
00489 }
00490 
00491 // ====================================================================== 
00492 void Amesos_Taucs::PrintTiming() const
00493 {
00494   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00495     return;
00496 
00497   double ConTime = GetTime(MtxConvTime_);
00498   double MatTime = GetTime(MtxRedistTime_);
00499   double VecTime = GetTime(VecRedistTime_);
00500   double SymTime = GetTime(SymFactTime_);
00501   double NumTime = GetTime(NumFactTime_);
00502   double SolTime = GetTime(SolveTime_);
00503 
00504   if (NumSymbolicFact_)
00505     SymTime /= NumSymbolicFact_;
00506 
00507   if (NumNumericFact_)
00508     NumTime /= NumNumericFact_;
00509 
00510   if (NumSolve_)
00511     SolTime /= NumSolve_;
00512 
00513   std::string p = "Amesos_Taucs : ";
00514   PrintLine();
00515 
00516   std::cout << p << "Time to convert matrix to Taucs format = "
00517        << ConTime << " (s)" << std::endl;
00518   std::cout << p << "Time to redistribute matrix = "
00519        << MatTime << " (s)" << std::endl;
00520   std::cout << p << "Time to redistribute vectors = "
00521        << VecTime << " (s)" << std::endl;
00522   std::cout << p << "Number of symbolic factorizations = "
00523        << NumSymbolicFact_ << std::endl;
00524   std::cout << p << "Time for sym fact = "
00525        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00526   std::cout << p << "Number of numeric factorizations = "
00527        << NumNumericFact_ << std::endl;
00528   std::cout << p << "Time for num fact = "
00529        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00530   std::cout << p << "Number of solve phases = "
00531        << NumSolve_ << std::endl;
00532   std::cout << p << "Time for solve = "
00533        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00534 
00535   PrintLine();
00536 
00537   return;
00538 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines