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