Amesos Package Browser (Single Doxygen Collection) Development
Amesos_CSparse.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_CSparse.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 
00036 
00037 using namespace Teuchos;
00038 
00039 //=============================================================================
00040 Amesos_CSparse::Amesos_CSparse(const Epetra_LinearProblem &prob) :
00041   UseTranspose_(false),
00042   Problem_(&prob),
00043   MtxConvTime_(-1),
00044   MtxRedistTime_(-1),
00045   VecRedistTime_(-1),
00046   SymFactTime_(-1),
00047   NumFactTime_(-1),
00048   SolveTime_(-1)
00049 {
00050     // Initialize data structures for CSparse
00051 }
00052 
00053 //=============================================================================
00054 Amesos_CSparse::~Amesos_CSparse() 
00055 {
00056     // Clean up
00057 
00058     //AMESOS_CHK_ERRV(CheckError(error));
00059     if ((verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
00060     if ((verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
00061 }
00062 
00063 //=============================================================================
00064 int Amesos_CSparse::ConvertToSerial() 
00065 {
00066     ResetTimer();
00067   const Epetra_RowMatrix* StdIndexMatrix_ ; 
00068   Epetra_CrsMatrix* CrsMatrixA_;
00069 
00070     if  ( Reindex_ ) {
00071         if ( Matrix_ == 0 ) AMESOS_CHK_ERR(-7);
00072 #ifdef HAVE_AMESOS_EPETRAEXT
00073         const Epetra_Map& OriginalMap = Matrix_->RowMatrixRowMap();
00074 
00075         const Epetra_Map &OriginalDomainMap = 
00076         UseTranspose()?GetProblem()->GetOperator()->OperatorRangeMap():
00077         GetProblem()->GetOperator()->OperatorDomainMap();
00078         const Epetra_Map &OriginalRangeMap = 
00079         UseTranspose()?GetProblem()->GetOperator()->OperatorDomainMap():
00080         GetProblem()->GetOperator()->OperatorRangeMap();
00081 
00082         StdIndex_ = rcp( new Amesos_StandardIndex( OriginalMap  ) );
00083         StdIndexDomain_ = rcp( new Amesos_StandardIndex( OriginalDomainMap  ) );
00084         StdIndexRange_ = rcp( new Amesos_StandardIndex( OriginalRangeMap  ) );
00085 
00086         CrsMatrixA_ = dynamic_cast<Epetra_CrsMatrix *>(Problem_->GetOperator());
00087         StdIndexMatrix_ = StdIndex_->StandardizeIndex( CrsMatrixA_ );
00088 #else
00089         std::cerr << "Amesos_CSparse requires EpetraExt to reindex matrices." << std::endl 
00090         AMESOS_CHK_ERR(-8);
00091 #endif
00092     } else { 
00093         StdIndexMatrix_ = Matrix_ ;
00094     }
00095 
00096     int NumGlobalRows = Matrix_->NumGlobalRows();
00097 
00098     // create a serial map
00099     int NumMyRows = 0;
00100     if (Comm().MyPID() == 0) 
00101         NumMyRows = NumGlobalRows;
00102 
00103     SerialMap_ = rcp(new Epetra_Map(-1, NumMyRows, 0, Comm()));
00104     if (SerialMap_.get() == 0)
00105     AMESOS_CHK_ERR(-1);
00106 
00107     Importer_ = rcp(new Epetra_Import(SerialMap(),StdIndexMatrix_->RowMatrixRowMap()));
00108     if (Importer_.get() == 0)
00109     AMESOS_CHK_ERR(-1);
00110 
00111     SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy, SerialMap(), 0));
00112     if (SerialCrsMatrix_.get() == 0)
00113     AMESOS_CHK_ERR(-1);
00114 
00115     AMESOS_CHK_ERR(SerialCrsMatrix().Import(*StdIndexMatrix_, Importer(), Add));
00116 
00117     AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
00118 
00119     SerialMatrix_ = rcp(SerialCrsMatrix_.get(), false);
00120 
00121     MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
00122 
00123     return 0;
00124 }
00125 
00126 //=============================================================================
00127 int Amesos_CSparse::ConvertToCSparse()
00128 {
00129   ResetTimer();
00130 
00131   if (Comm().MyPID() == 0) 
00132   {
00133     csMatrix.p = (ptrdiff_t *) malloc((SerialMatrix().NumMyRows()+1)*sizeof(ptrdiff_t));
00134     csMatrix.i = (ptrdiff_t *) malloc(SerialMatrix().NumMyNonzeros()*sizeof(ptrdiff_t));
00135     csMatrix.x = (double *) malloc(SerialMatrix().NumMyNonzeros()*sizeof(double));
00136     csMatrix.nzmax = SerialMatrix().NumMyNonzeros();
00137     csMatrix.m = SerialMatrix().NumMyRows();
00138     csMatrix.n = SerialMatrix().NumMyRows();
00139     csMatrix.nz = -1;
00140 
00141     int MaxNumEntries = SerialMatrix().MaxNumEntries();
00142     std::vector<int>    Indices(MaxNumEntries);
00143     std::vector<double> Values(MaxNumEntries);
00144 
00145     csMatrix.p[0] = 0;
00146     int count = 0;
00147 
00148     for (int i = 0 ; i < SerialMatrix().NumMyRows() ; ++i)
00149     {
00150       int ierr, NumEntriesThisRow;
00151       ierr = SerialMatrix().ExtractMyRowCopy(i, MaxNumEntries, 
00152                                              NumEntriesThisRow, 
00153                                              &Values[0], &Indices[0]);
00154       if (ierr < 0)
00155         AMESOS_CHK_ERR(ierr);
00156 
00157       csMatrix.p[i + 1] = csMatrix.p[i] + NumEntriesThisRow;
00158 
00159       for (int j = 0 ; j < NumEntriesThisRow ; ++j)
00160       {
00161         if (Indices[j] == i) 
00162           Values[j] += AddToDiag_;
00163 
00164         csMatrix.i[count] = Indices[j];
00165         csMatrix.x[count] = Values[j];
00166         ++count;
00167       }
00168     }
00169     
00170     if (count != SerialMatrix().NumMyNonzeros())
00171       AMESOS_CHK_ERR(-1); // something wrong here*/
00172 
00173     // TODO: Avoid this transpose once we have cs_sptsolve()
00174     csTranMatrix = cs_transpose(&csMatrix, 1);
00175 
00176     free(csMatrix.p);
00177     free(csMatrix.i);
00178     free(csMatrix.x);
00179   }
00180 
00181   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
00182 
00183   return 0;
00184 }
00185 
00186 //=============================================================================
00187 int Amesos_CSparse::SetParameters( Teuchos::ParameterList &ParameterList) 
00188 {
00189   // retrive general parameters
00190 
00191   SetStatusParameters( ParameterList );
00192 
00193   SetControlParameters( ParameterList );
00194 
00195   // retrive CSparse's specific parameters
00196 
00197   if (ParameterList.isSublist("CSparse")) 
00198   {
00199       // set solver specific parameters here.
00200 
00201   }
00202   
00203   return 0;
00204 }
00205 
00206 //=============================================================================
00207 int Amesos_CSparse::PerformSymbolicFactorization() 
00208 {
00209   ResetTimer();
00210 
00211   if (Comm().MyPID() == 0) 
00212   {
00213     // ============================================================== //
00214     // Setup CSparse parameteres and call symbolic factorization.
00215     // ============================================================== //
00216     int error = 0;
00217 
00218     // Call Csparse here.
00219     csSymbolic = cs_sqr(2, csTranMatrix, 0);
00220 
00221     AMESOS_CHK_ERR(CheckError(error));
00222   }
00223 
00224   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00225 
00226   return 0;
00227 }
00228 
00229 //=============================================================================
00230 int Amesos_CSparse::PerformNumericFactorization( ) 
00231 {
00232   ResetTimer();
00233 
00234   if (Comm().MyPID() == 0) 
00235   {
00236     int error = 0;
00237     // Call Csparse here.
00238     csNumeric = cs_lu(csTranMatrix, csSymbolic, 1e-15);
00239 
00240     AMESOS_CHK_ERR(CheckError(error));
00241   }
00242 
00243   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00244 
00245   return 0;
00246 }
00247 
00248 //=============================================================================
00249 bool Amesos_CSparse::MatrixShapeOK() const 
00250 {
00251   bool OK = true;
00252 
00253   if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
00254       GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) 
00255   {
00256     OK = false;
00257   }
00258   return OK;
00259 }
00260 
00261 //=============================================================================
00262 int Amesos_CSparse::SymbolicFactorization() 
00263 {
00264   IsSymbolicFactorizationOK_ = false;
00265   IsNumericFactorizationOK_ = false;
00266 
00267   CreateTimer(Comm());
00268 
00269   ++NumSymbolicFact_;
00270 
00271   Matrix_ = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00272   Map_ = &(Matrix_->RowMatrixRowMap());
00273 
00274   // =========================================================== //
00275   // redistribute and create all import/export objects only      //
00276   // if more than one processor is used. Otherwise simply set    //
00277   // dummy pointers to Matrix() and Map(), without giving the    //
00278   // ownership to the smart pointer.                             //
00279   // =========================================================== //
00280 
00281   if (Comm().NumProc() != 1) 
00282     ConvertToSerial();
00283   else
00284   {
00285     SerialMap_ = rcp(const_cast<Epetra_Map*>(&Map()), false);
00286     SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(&Matrix()), false);
00287   }
00288 
00289   // =========================================================== //
00290   // Only on processor zero, convert the matrix into CSR format, //
00291   // as required by CSparse.                                     //
00292   // =========================================================== //
00293 
00294   ConvertToCSparse();
00295 
00296   PerformSymbolicFactorization();
00297 
00298   IsSymbolicFactorizationOK_ = true;
00299 
00300   return(0);
00301 }
00302 
00303 //=============================================================================
00304 int Amesos_CSparse::NumericFactorization() 
00305 {
00306   IsNumericFactorizationOK_ = false;
00307 
00308   if (IsSymbolicFactorizationOK_ == false)
00309     AMESOS_CHK_ERR(SymbolicFactorization());
00310 
00311   ++NumNumericFact_;
00312 
00313   // FIXME: this must be checked, now all the matrix is shipped twice here
00314   ConvertToSerial();
00315   ConvertToCSparse();
00316 
00317   PerformNumericFactorization();
00318 
00319   IsNumericFactorizationOK_ = true;
00320 
00321   return(0);
00322 }
00323 
00324 //=============================================================================
00325 int Amesos_CSparse::Solve() 
00326 {
00327   Epetra_MultiVector* vecX = 0 ;
00328   Epetra_MultiVector* vecB = 0 ;
00329 
00330 #ifdef HAVE_AMESOS_EPETRAEXT
00331   Teuchos::RCP<Epetra_MultiVector> vecX_rcp;
00332   Teuchos::RCP<Epetra_MultiVector> vecB_rcp;
00333 #endif
00334 
00335   if (IsNumericFactorizationOK_ == false)
00336     AMESOS_CHK_ERR(NumericFactorization());
00337 
00338   Epetra_MultiVector* X = Problem_->GetLHS();
00339   Epetra_MultiVector* B = Problem_->GetRHS();
00340 
00341   if ((X == 0) || (B == 0))
00342     AMESOS_CHK_ERR(-1); 
00343 
00344   int NumVectors = X->NumVectors();
00345   if (NumVectors != B->NumVectors())
00346     AMESOS_CHK_ERR(-1); 
00347 
00348     if ( Reindex_ ) { 
00349 #ifdef HAVE_AMESOS_EPETRAEXT
00350       vecX_rcp = StdIndexDomain_->StandardizeIndex( *X ) ;
00351       vecB_rcp = StdIndexRange_->StandardizeIndex( *B ) ;
00352 
00353       vecX = &*vecX_rcp;
00354       vecB = &*vecB_rcp;
00355 #else
00356       AMESOS_CHK_ERR( -13 ) ; // Amesos_CSparse can't handle non-standard indexing without EpetraExt 
00357 #endif
00358     } else {
00359       vecX = X ;
00360       vecB = B ;
00361     } 
00362     
00363   // vectors with SerialMap_
00364   Epetra_MultiVector* SerialB;
00365   Epetra_MultiVector* SerialX;
00366 
00367   ResetTimer();
00368 
00369   SerialX = new Epetra_MultiVector(SerialMap(),NumVectors);
00370   SerialB = new Epetra_MultiVector(SerialMap(),NumVectors);
00371 
00372   SerialB->Import(*vecB,Importer(),Insert);
00373 
00374   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00375 
00376   ResetTimer();
00377 
00378   if (Comm().MyPID() == 0) 
00379   {
00380     double* SerialXValues;
00381     double* SerialBValues;
00382     int LDA;
00383 
00384     AMESOS_CHK_ERR(SerialX->ExtractView(&SerialXValues,&LDA));
00385 
00386     // FIXME: check LDA
00387     AMESOS_CHK_ERR(SerialB->ExtractView(&SerialBValues,&LDA));
00388 
00389     int error = 0;
00390     int n = SerialMatrix().NumMyRows();
00391 
00392     double *x = (double *) malloc(n * sizeof(double));
00393     // TODO: OMP here ??
00394     for (int i = 0 ; i < NumVectors ; ++i)
00395     {
00396         // Call Csparse here
00397         cs_ipvec(csNumeric->pinv, SerialBValues+i*n, x, n);
00398         cs_lsolve(csNumeric->L, x);
00399         cs_usolve(csNumeric->U, x);
00400         cs_ipvec(csSymbolic->q, x, SerialXValues+i*n, n);
00401     }
00402     free(x);
00403 
00404     AMESOS_CHK_ERR(CheckError(error));
00405   }
00406 
00407   SolveTime_ = AddTime("Total solve time", SolveTime_);
00408 
00409   //  Copy X back to the original vector
00410 
00411   ResetTimer();
00412 
00413   vecX->Export(*SerialX, Importer(), Insert);
00414   delete SerialB;
00415   delete SerialX;
00416 
00417   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00418 
00419   if (ComputeTrueResidual_)
00420     ComputeTrueResidual(Matrix(), *X, *B, UseTranspose(), "Amesos_CSparse");
00421 
00422   if (ComputeVectorNorms_)
00423     ComputeVectorNorms(*X, *B, "Amesos_CSparse");
00424 
00425   ++NumSolve_;
00426 
00427   return(0) ;
00428 }
00429 
00430 // ====================================================================== 
00431 void Amesos_CSparse::PrintStatus() const
00432 {
00433   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00434     return;
00435 
00436   std::string p = "Amesos_CSparse : ";
00437   PrintLine();
00438 
00439   PrintLine();
00440 
00441   return;
00442 }
00443 
00444 // ====================================================================== 
00445 void Amesos_CSparse::PrintTiming() const
00446 {
00447   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00448     return;
00449 
00450   double ConTime = GetTime(MtxConvTime_);
00451   double MatTime = GetTime(MtxRedistTime_);
00452   double VecTime = GetTime(VecRedistTime_);
00453   double SymTime = GetTime(SymFactTime_);
00454   double NumTime = GetTime(NumFactTime_);
00455   double SolTime = GetTime(SolveTime_);
00456 
00457   if (NumSymbolicFact_)
00458     SymTime /= NumSymbolicFact_;
00459 
00460   if (NumNumericFact_)
00461     NumTime /= NumNumericFact_;
00462 
00463   if (NumSolve_)
00464     SolTime /= NumSolve_;
00465 
00466   std::string p = "Amesos_CSparse : ";
00467   PrintLine();
00468 
00469   std::cout << p << "Time to convert matrix to Csparse format = "
00470        << ConTime << " (s)" << std::endl;
00471   std::cout << p << "Time to redistribute matrix = "
00472        << MatTime << " (s)" << std::endl;
00473   std::cout << p << "Time to redistribute vectors = "
00474        << VecTime << " (s)" << std::endl;
00475   std::cout << p << "Number of symbolic factorizations = "
00476        << NumSymbolicFact_ << std::endl;
00477   std::cout << p << "Time for sym fact = "
00478        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00479   std::cout << p << "Number of numeric factorizations = "
00480        << NumNumericFact_ << std::endl;
00481   std::cout << p << "Time for num fact = "
00482        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00483   std::cout << p << "Number of solve phases = "
00484        << NumSolve_ << std::endl;
00485   std::cout << p << "Time for solve = "
00486        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00487 
00488   PrintLine();
00489 
00490   return;
00491 }
00492 
00493 // ====================================================================== 
00494 int Amesos_CSparse::CheckError(const int error) const
00495 {
00496   if (!error)
00497     return 0;
00498   
00499   std::cerr << "Amesos: CSparse returned error code " << error << std::endl;
00500 
00501   AMESOS_RETURN(error);
00502 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines