Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Paraklete.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_Paraklete.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 #include "Amesos_Support.h"
00036 extern "C" {
00037 #include "amesos_paraklete_decl.h"
00038 }
00039 
00040 // The following hack allows assert to work even though it has been turned off by paraklete.h bug #1952 
00041 
00042 #undef  _ASSERT_H
00043 #undef  assert
00044 #undef __ASSERT_VOID_CAST
00045 #undef NDEBUG
00046 #include <assert.h>
00047 
00048 namespace {
00049 
00050 using Teuchos::RCP;
00051 
00052 template<class T, class DeleteFunctor>
00053 class DeallocFunctorDeleteWithCommon
00054 {
00055 public:
00056   DeallocFunctorDeleteWithCommon(
00057          const RCP<paraklete_common>  &common
00058          ,DeleteFunctor                        deleteFunctor
00059          )
00060     : common_(common), deleteFunctor_(deleteFunctor)
00061     {}
00062   typedef T ptr_t;
00063   void free( T* ptr ) {
00064     if(ptr) deleteFunctor_(&ptr,&*common_);
00065   }
00066 private:
00067   Teuchos::RCP<paraklete_common> common_;
00068   DeleteFunctor deleteFunctor_;
00069   DeallocFunctorDeleteWithCommon(); // Not defined and not to be called!
00070 };
00071 
00072 template<class T, class DeleteFunctor>
00073 DeallocFunctorDeleteWithCommon<T,DeleteFunctor>
00074 deallocFunctorDeleteWithCommon(
00075              const RCP<paraklete_common>  &common
00076              ,DeleteFunctor                        deleteFunctor
00077              )
00078 {
00079   return DeallocFunctorDeleteWithCommon<T,DeleteFunctor>(common,deleteFunctor);
00080 }
00081 
00082 
00083 } // namespace 
00084 
00085 
00086 
00087 class Amesos_Paraklete_Pimpl {
00088 public:
00089 
00090 
00091   cholmod_sparse pk_A_ ;
00092   Teuchos::RCP<paraklete_symbolic> LUsymbolic_ ;
00093   Teuchos::RCP<paraklete_numeric> LUnumeric_ ;
00094   Teuchos::RCP<paraklete_common> common_;
00095 
00096 
00097 
00098 } ;
00099 
00100 //=============================================================================
00101 Amesos_Paraklete::Amesos_Paraklete(const Epetra_LinearProblem &prob ) :
00102   PrivateParakleteData_( rcp( new Amesos_Paraklete_Pimpl() ) ),
00103   ParakleteComm_(0),
00104   CrsMatrixA_(0),
00105   TrustMe_(false),
00106   UseTranspose_(false),
00107   Problem_(&prob),
00108   MtxConvTime_(-1),
00109   MtxRedistTime_(-1),
00110   VecRedistTime_(-1),
00111   SymFactTime_(-1),
00112   NumFactTime_(-1),
00113   SolveTime_(-1),
00114   OverheadTime_(-1)
00115 {
00116 
00117   // MS // move declaration of Problem_ above because I need it
00118   // MS // set up before calling Comm()
00119   MaxProcesses_ = -3;    // Use all processes unless the user requests otherwise
00120   Teuchos::ParameterList ParamList ;
00121   SetParameters( ParamList ) ;
00122 }
00123 
00124 //=============================================================================
00125 Amesos_Paraklete::~Amesos_Paraklete(void) {
00126 
00127 
00128     if(ParakleteComm_)  {
00129       MPI_Comm_free(&ParakleteComm_);
00130       ParakleteComm_ = 0 ; 
00131     }
00132 
00133   // print out some information if required by the user
00134   if( (verbose_ && PrintTiming_) || verbose_ == 2 ) PrintTiming();
00135   if( (verbose_ && PrintStatus_) || verbose_ == 2 ) PrintStatus();
00136 }
00137 
00138 //=============================================================================
00139 int Amesos_Paraklete::ExportToSerial() 
00140 {
00141   if (  AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros()) { 
00142     std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization().  " ;
00143     AMESOS_CHK_ERR( -2 );
00144   }
00145   if ( UseDataInPlace_ == 0 ) {
00146     assert ( RowMatrixA_ != 0 ) ; 
00147     if (  AddZeroToDiag_==0 && numentries_ != RowMatrixA_->NumGlobalNonzeros()) { 
00148       std::cerr << " The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization().  " ;
00149       AMESOS_CHK_ERR( -2 );
00150     }
00151 
00152 #ifdef HAVE_AMESOS_EPETRAEXT
00153     if ( transposer_.get() != 0  ) {
00154       //      int OriginalTracebackMode = CrsMatrixA_->GetTracebackMode();
00155       //      CrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) );
00156       transposer_->fwd() ;
00157       //      CrsMatrixA_->SetTracebackMode( OriginalTracebackMode );
00158     }
00159 #endif
00160     assert ( ImportToSerial_.get() != 0 ) ; 
00161     AMESOS_CHK_ERR(SerialCrsMatrixA_->Import(*StdIndexMatrix_, 
00162                *ImportToSerial_, Insert ));
00163 
00164     if (AddZeroToDiag_ ) { 
00165       int OriginalTracebackMode = SerialCrsMatrixA_->GetTracebackMode() ; 
00166       SerialCrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) ) ; // ExportToSerial is called both by PerformSymbolicFactorization() and PerformNumericFactorization().  When called by the latter, the call to insertglobalvalues is both unnecessary and illegal.  Fortunately, Epetra allows us to ignore the error message by setting the traceback mode to 0.
00167 
00168       //
00169       //  Add 0.0 to each diagonal entry to avoid empty diagonal entries;
00170       //
00171       double zero = 0.0;
00172       for ( int i = 0 ; i < SerialMap_->NumGlobalElements(); i++ ) 
00173   if ( SerialCrsMatrixA_->LRID(i) >= 0 ) 
00174     SerialCrsMatrixA_->InsertGlobalValues( i, 1, &zero, &i ) ;
00175       SerialCrsMatrixA_->SetTracebackMode( OriginalTracebackMode ) ; 
00176     }
00177     AMESOS_CHK_ERR(SerialCrsMatrixA_->FillComplete());
00178     //AMESOS_CHK_ERR(SerialCrsMatrixA_->OptimizeStorage());
00179 
00180     if( !AddZeroToDiag_ && numentries_ != SerialMatrix_->NumGlobalNonzeros()) {
00181       std::cerr << " Amesos_Paraklete cannot handle this matrix.  " ;
00182       if ( Reindex_ ) {
00183   std::cerr << "Unknown error" << std::endl ; 
00184   AMESOS_CHK_ERR( -5 );
00185       } else {
00186   std::cerr << " Try setting the Reindex parameter to true. " << std::endl ; 
00187 #ifndef HAVE_AMESOS_EPETRAEXT
00188   std::cerr << " You will need to rebuild the Amesos library with the EpetraExt library to use the reindex feature" << std::endl ; 
00189   std::cerr << " To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
00190 #endif
00191   AMESOS_CHK_ERR( -3 );
00192       }
00193     }
00194   }
00195   return 0;
00196 }
00197 //=============================================================================
00198 //
00199 // CreateLocalMatrixAndExporters() is called only by SymbolicFactorization() 
00200 // for CrsMatrix objects.  All objects should be created here.  No assumptions 
00201 // are made about the input operator.  I.e. it can be completely different from 
00202 // the matrix at the time of the previous call to SymbolicFactorization(). 
00203 //
00204 
00205 int Amesos_Paraklete::CreateLocalMatrixAndExporters() 
00206 {
00207   ResetTimer(0);
00208 
00209   RowMatrixA_ = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00210   if (RowMatrixA_ == 0) AMESOS_CHK_ERR(-1);
00211   
00212   const Epetra_Map &OriginalMatrixMap = RowMatrixA_->RowMatrixRowMap() ;
00213   const Epetra_Map &OriginalDomainMap = 
00214     UseTranspose()?GetProblem()->GetOperator()->OperatorRangeMap():
00215     GetProblem()->GetOperator()->OperatorDomainMap();
00216   const Epetra_Map &OriginalRangeMap = 
00217     UseTranspose()?GetProblem()->GetOperator()->OperatorDomainMap():
00218     GetProblem()->GetOperator()->OperatorRangeMap();
00219   
00220   NumGlobalElements_ = RowMatrixA_->NumGlobalRows();
00221   numentries_ = RowMatrixA_->NumGlobalNonzeros();
00222   assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
00223   
00224   //
00225   //  Create a serial matrix
00226   //
00227   assert( NumGlobalElements_ == OriginalMatrixMap.NumGlobalElements() ) ;
00228   int NumMyElements_ = 0 ;
00229   if (MyPID_==0) NumMyElements_ = NumGlobalElements_;
00230   
00231   //
00232   //  UseDataInPlace_ is set to 1 (true) only if everything is perfectly
00233   //  normal.  Anything out of the ordinary reverts to the more expensive
00234   //  path. 
00235   //
00236   UseDataInPlace_ = ( OriginalMatrixMap.NumMyElements() ==
00237           OriginalMatrixMap.NumGlobalElements() )?1:0;
00238   if ( ! OriginalRangeMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ; 
00239   if ( ! OriginalDomainMap.SameAs( OriginalMatrixMap ) ) UseDataInPlace_ = 0 ; 
00240   if ( AddZeroToDiag_ ) UseDataInPlace_ = 0 ; 
00241   Comm().Broadcast( &UseDataInPlace_, 1, 0 ) ;
00242   
00243   UseDataInPlace_ = 0 ; // bug - remove this someday.
00244     
00245   //
00246   //  Reindex matrix if necessary (and possible - i.e. CrsMatrix)
00247   //
00248   //  For now, since I don't know how to determine if we need to reindex the matrix, 
00249   //  we allow the user to control re-indexing through the "Reindex" parameter.
00250   //
00251   CrsMatrixA_ = dynamic_cast<Epetra_CrsMatrix *>(Problem_->GetOperator());
00252   Epetra_CrsMatrix* CcsMatrixA = 0 ; 
00253     
00254   //
00255   //  CcsMatrixA points to a matrix in Compressed Column Format 
00256   //  i.e. the format needed by Paraklete.  If we are solving  
00257   //  A^T x = b, CcsMatrixA = CrsMatrixA_, otherwise we must 
00258   //  transpose the matrix.
00259   //
00260   if (UseTranspose()) {
00261     CcsMatrixA = CrsMatrixA_ ;
00262   } else {
00263     if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR( -7 );   //  Amesos_Paraklete only supports CrsMatrix objects in the non-transpose case
00264 #ifdef HAVE_AMESOS_EPETRAEXT
00265     bool MakeDataContiguous = true;
00266     transposer_ = rcp ( new EpetraExt::RowMatrix_Transpose( MakeDataContiguous ));
00267     
00268     int OriginalTracebackMode = CrsMatrixA_->GetTracebackMode();
00269     CrsMatrixA_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) );
00270     
00271     CcsMatrixA = &(dynamic_cast<Epetra_CrsMatrix&>(((*transposer_)(*CrsMatrixA_))));
00272     CrsMatrixA_->SetTracebackMode( OriginalTracebackMode );
00273 #else
00274     std::cerr << "Amesos_Paraklete requires the EpetraExt library to solve non-transposed problems. " << std::endl ; 
00275     std::cerr << " To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
00276     AMESOS_CHK_ERR( -3 );
00277 #endif
00278   }
00279   
00280   
00281   if  ( Reindex_ ) {
00282     if ( CcsMatrixA == 0 ) AMESOS_CHK_ERR(-4);
00283 #ifdef HAVE_AMESOS_EPETRAEXT
00284     const Epetra_Map& OriginalMap = CcsMatrixA->RowMap();
00285     StdIndex_ = rcp( new Amesos_StandardIndex( OriginalMap  ) );
00286     //const Epetra_Map& OriginalColMap = CcsMatrixA->RowMap();
00287     StdIndexDomain_ = rcp( new Amesos_StandardIndex( OriginalDomainMap  ) );
00288     StdIndexRange_ = rcp( new Amesos_StandardIndex( OriginalRangeMap  ) );
00289     
00290     StdIndexMatrix_ = StdIndex_->StandardizeIndex( CcsMatrixA );
00291 #else
00292     std::cerr << "Amesos_Paraklete requires EpetraExt to reindex matrices." << std::endl 
00293    <<  " Please rebuild with the EpetraExt library by adding --enable-epetraext to your configure invocation" << std::endl ;
00294     AMESOS_CHK_ERR(-4);
00295 #endif
00296   } else { 
00297     if ( UseTranspose() ){ 
00298       StdIndexMatrix_ = RowMatrixA_ ;
00299     } else {
00300       StdIndexMatrix_ = CcsMatrixA ;
00301     } 
00302   }
00303   
00304   //
00305   //  At this point, StdIndexMatrix_ points to a matrix with 
00306   //  standard indexing.  
00307   //
00308 
00309   //
00310   //  Convert Original Matrix to Serial (if it is not already)
00311   //
00312   if (UseDataInPlace_ == 1) {
00313     SerialMatrix_ = StdIndexMatrix_;
00314   } else {
00315     SerialMap_ = rcp(new Epetra_Map(NumGlobalElements_, NumMyElements_, 0, Comm()));
00316     
00317     if ( Problem_->GetRHS() ) 
00318       NumVectors_ = Problem_->GetRHS()->NumVectors() ; 
00319     else
00320       NumVectors_ = 1 ; 
00321     SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
00322     SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
00323 
00324     ImportToSerial_ = rcp(new Epetra_Import ( *SerialMap_, StdIndexMatrix_->RowMatrixRowMap() ) );
00325     if (ImportToSerial_.get() == 0) AMESOS_CHK_ERR(-5);
00326 
00327     SerialCrsMatrixA_ = rcp( new Epetra_CrsMatrix(Copy, *SerialMap_, 0) ) ;
00328     SerialMatrix_ = &*SerialCrsMatrixA_ ;
00329   }
00330   
00331   MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
00332   
00333   return(0);
00334 }
00335 
00336 //=============================================================================
00337 //
00338 //  See also pre and post conditions in Amesos_Paraklete.h
00339 //  Preconditions:
00340 //    firsttime specifies that this is the first time that 
00341 //    ConertToParakleteCrs has been called, i.e. in symbolic factorization.  
00342 //    No data allocation should happen unless firsttime=true.
00343 //    SerialMatrix_ points to the matrix to be factored and solved
00344 //    NumGlobalElements_ has been set to the dimension of the matrix
00345 //    numentries_ has been set to the number of non-zeros in the matrix
00346 //      (i.e. CreateLocalMatrixAndExporters() has been callded)
00347 //
00348 //  Postconditions:
00349 //    pk_A_ contains the matrix as Paraklete needs it
00350 //
00351 //
00352 int Amesos_Paraklete::ConvertToParakleteCRS(bool firsttime)
00353 {
00354   ResetTimer(0);
00355   
00356   //
00357   //  Convert matrix to the form that Klu expects (Ap, Ai, VecAval)
00358   //
00359 
00360   if (MyPID_==0) {
00361     assert( NumGlobalElements_ == SerialMatrix_->NumGlobalRows());
00362     assert( NumGlobalElements_ == SerialMatrix_->NumGlobalCols());
00363 
00364     if ( ! AddZeroToDiag_ ) {
00365       assert( numentries_ == SerialMatrix_->NumGlobalNonzeros()) ;
00366     } else {
00367       numentries_ = SerialMatrix_->NumGlobalNonzeros() ;
00368     }
00369 
00370     Epetra_CrsMatrix *CrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(SerialMatrix_);
00371     bool StorageOptimized = ( CrsMatrix != 0 && CrsMatrix->StorageOptimized() );
00372 
00373     if ( AddToDiag_ != 0.0 ) StorageOptimized = false ;
00374 
00375     if ( firsttime ) { 
00376       Ap.resize( NumGlobalElements_+1 );
00377       Ai.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
00378       if ( ! StorageOptimized ) { 
00379   VecAval.resize( EPETRA_MAX( NumGlobalElements_, numentries_) ) ;
00380   Aval = &VecAval[0];
00381       }
00382     }
00383 
00384     double *RowValues;
00385     int *ColIndices;
00386     int NumEntriesThisRow;
00387 
00388     if( StorageOptimized ) {
00389       if ( firsttime ) {
00390   Ap[0] = 0;
00391   for ( int MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
00392     if( CrsMatrix->
00393         ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
00394         ColIndices ) != 0 ) 
00395       AMESOS_CHK_ERR( -6 );
00396           for ( int j=0; j < NumEntriesThisRow; j++ ) {
00397             Ai[Ap[MyRow]+j] = ColIndices[j];
00398           }
00399     if ( MyRow == 0 ) {
00400       Aval = RowValues ;
00401     }
00402     Ap[MyRow+1] = Ap[MyRow] + NumEntriesThisRow ;
00403   }
00404       }
00405     } else { 
00406       
00407       int Ai_index = 0 ;
00408       int MyRow;
00409       
00410       int MaxNumEntries_ = SerialMatrix_->MaxNumEntries();
00411       if ( firsttime && CrsMatrix == 0 ) {
00412   ColIndicesV_.resize(MaxNumEntries_);
00413   RowValuesV_.resize(MaxNumEntries_);
00414       }
00415       
00416       for ( MyRow = 0; MyRow <NumGlobalElements_; MyRow++ ) {
00417   if ( CrsMatrix != 0 ) {
00418     if( CrsMatrix->
00419         ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
00420         ColIndices ) != 0 ) 
00421       AMESOS_CHK_ERR( -6 );
00422 
00423   } else {
00424     if( SerialMatrix_->
00425         ExtractMyRowCopy( MyRow, MaxNumEntries_,
00426               NumEntriesThisRow, &RowValuesV_[0],
00427               &ColIndicesV_[0] ) != 0 ) 
00428       AMESOS_CHK_ERR( -6 );
00429 
00430     RowValues =  &RowValuesV_[0];
00431     ColIndices = &ColIndicesV_[0];
00432   }
00433   
00434   if ( firsttime ) {
00435     Ap[MyRow] = Ai_index ;
00436     for ( int j = 0; j < NumEntriesThisRow; j++ ) {
00437       Ai[Ai_index] = ColIndices[j] ;
00438       //  VecAval[Ai_index] = RowValues[j] ;      //  We have to do this because of the hacks to get around bug #1502 
00439       if (ColIndices[j] == MyRow) {
00440         VecAval[Ai_index] += AddToDiag_;    
00441       }
00442       Ai_index++;
00443     }
00444   } else { 
00445     for ( int j = 0; j < NumEntriesThisRow; j++ ) {
00446       VecAval[Ai_index] = RowValues[j] ;     
00447       if (ColIndices[j] == MyRow) {
00448         VecAval[Ai_index] += AddToDiag_;   
00449       }
00450       Ai_index++;
00451     }
00452   }
00453       }
00454       Ap[MyRow] = Ai_index ;
00455     }
00456   PrivateParakleteData_->pk_A_.nrow = NumGlobalElements_ ; 
00457   cholmod_sparse& pk_A =  PrivateParakleteData_->pk_A_ ;
00458   pk_A.nrow = NumGlobalElements_ ; 
00459   pk_A.ncol = NumGlobalElements_ ; 
00460   pk_A.nzmax = numentries_ ; 
00461   pk_A.p = &Ap[0] ;
00462   pk_A.i = &Ai[0] ;
00463   pk_A.nz = 0; 
00464   if ( firsttime ) { 
00465     pk_A.x = 0 ; 
00466     pk_A.xtype = CHOLMOD_PATTERN ; 
00467   }
00468   else
00469   {
00470     pk_A.x = Aval ; 
00471     pk_A.xtype = CHOLMOD_REAL ; 
00472   }
00473 
00474   pk_A.z = 0 ; 
00475   pk_A.stype = 0 ; //  symmetric 
00476   pk_A.itype = CHOLMOD_LONG ; 
00477   pk_A.dtype = CHOLMOD_DOUBLE ; 
00478   pk_A.sorted = 0 ; 
00479   pk_A.packed = 1 ; 
00480   }
00481 
00482   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
00483 
00484   return 0;
00485 }
00486 
00487 //=============================================================================
00488 int Amesos_Paraklete::SetParameters( Teuchos::ParameterList &ParameterList ) {
00489   
00490   // ========================================= //
00491   // retrive PARAKLETE's parameters from list.       //
00492   // default values defined in the constructor //
00493   // ========================================= //
00494   
00495   // retrive general parameters
00496   
00497   SetStatusParameters( ParameterList );
00498   SetControlParameters( ParameterList );
00499   
00500   
00501   
00502   if (ParameterList.isParameter("TrustMe") ) 
00503     TrustMe_ = ParameterList.get<bool>( "TrustMe" );
00504   
00505 #if 0
00506   
00507   unused for now
00508     
00509     if (ParameterList.isSublist("Paraklete") ) {
00510       Teuchos::ParameterList ParakleteParams = ParameterList.sublist("Paraklete") ;
00511     }
00512 #endif
00513   
00514   return 0;
00515 }
00516 
00517 
00518 //=============================================================================
00519 int Amesos_Paraklete::PerformSymbolicFactorization() 
00520 {
00521   ResetTimer(0);
00522   
00523   if ( IamInGroup_ ) 
00524     PrivateParakleteData_->LUsymbolic_ =
00525       rcp( amesos_paraklete_analyze ( &PrivateParakleteData_->pk_A_, &*PrivateParakleteData_->common_ )
00526      ,deallocFunctorDeleteWithCommon<paraklete_symbolic>(PrivateParakleteData_->common_,
00527                      amesos_paraklete_free_symbolic)
00528      ,true
00529      );
00530   
00531   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_, 0);
00532   
00533   return 0;
00534 }
00535 
00536 //=============================================================================
00537 int Amesos_Paraklete::PerformNumericFactorization( ) 
00538 {
00539   // Changed this; it was "if (!TrustMe)...
00540   // The behavior is not intuitive. Maybe we should introduce a new
00541   // parameter, FastSolvers or something like that, that does not perform
00542   // any AddTime, ResetTimer, GetTime.
00543   
00544   ResetTimer(0);
00545   
00546   //bool factor_with_pivoting = true ;
00547   
00548   if ( IamInGroup_ ) 
00549     PrivateParakleteData_->LUnumeric_ =
00550       rcp( amesos_paraklete_factorize ( &PrivateParakleteData_->pk_A_,
00551                                  &*PrivateParakleteData_->LUsymbolic_, 
00552                                  &*PrivateParakleteData_->common_ ) 
00553            ,deallocFunctorDeleteWithCommon<paraklete_numeric>(PrivateParakleteData_->common_,amesos_paraklete_free_numeric)
00554      ,true
00555      );
00556 
00557   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
00558   
00559   return 0;
00560 }
00561 
00562 //=============================================================================
00563 bool Amesos_Paraklete::MatrixShapeOK() const {
00564   bool OK = true;
00565   
00566   // Comment by Tim:  The following code seems suspect.  The variable "OK"
00567   // is not set if the condition is true.
00568   // Does the variable "OK" default to true?
00569   if ( GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
00570        GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) {
00571     OK = false;
00572   }
00573   return OK;
00574 }
00575 
00576 
00577 extern "C" {
00578   void my_handler (int status, char *file, int line, char *msg)
00579   {
00580     printf ("Error handler: %s %d %d: %s\n", file, line, status, msg) ;
00581     if (status != CHOLMOD_OK)
00582       {
00583   fprintf (stderr, "\n\n*********************************************\n");
00584   fprintf (stderr, "**** Test failure: %s %d %d %s\n", file, line,
00585      status, msg) ;
00586   fprintf (stderr, "*********************************************\n\n");
00587   fflush (stderr) ;
00588   fflush (stdout) ;
00589       }
00590   }
00591 }  
00592 
00593 //=============================================================================
00594 int Amesos_Paraklete::SymbolicFactorization() 
00595 {
00596   MyPID_    = Comm().MyPID();
00597   NumProcs_ = Comm().NumProc();
00598   
00599   IsSymbolicFactorizationOK_ = false;
00600   IsNumericFactorizationOK_ = false;
00601   
00602 #ifdef HAVE_AMESOS_EPETRAEXT
00603   transposer_ = static_cast<Teuchos::ENull>( 0 ); 
00604 #endif
00605   
00606   CreateTimer(Comm(), 2);
00607   
00608   ResetTimer(1);
00609   
00610   // "overhead" time for the following method is considered here
00611   AMESOS_CHK_ERR( CreateLocalMatrixAndExporters() ) ;
00612   assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
00613   
00614   
00615   SetMaxProcesses(MaxProcesses_, *RowMatrixA_);
00616   
00617   //
00618   //  Perform checks in SymbolicFactorization(), but none in 
00619   //  NumericFactorization() or Solve()
00620   //
00621   assert( ! TrustMe_ ) ;
00622   if ( TrustMe_ ) { 
00623     if ( CrsMatrixA_ == 0 ) AMESOS_CHK_ERR(10 );
00624     if( UseDataInPlace_ != 1 ) AMESOS_CHK_ERR( 10 ) ;
00625     if( Reindex_ )  AMESOS_CHK_ERR( 10 ) ;
00626     if( ! Problem_->GetLHS() )  AMESOS_CHK_ERR( 10 ) ;
00627     if( ! Problem_->GetRHS() )  AMESOS_CHK_ERR( 10 ) ;
00628     if( ! Problem_->GetLHS()->NumVectors() ) AMESOS_CHK_ERR( 10 ) ;
00629     if( ! Problem_->GetRHS()->NumVectors() ) AMESOS_CHK_ERR( 10 ) ; 
00630     SerialB_ = Problem_->GetRHS() ;
00631     SerialX_ = Problem_->GetLHS() ;
00632     NumVectors_ = SerialX_->NumVectors();
00633     if (MyPID_ == 0) {
00634       AMESOS_CHK_ERR(SerialX_->ExtractView(&SerialXBvalues_,&SerialXlda_ ));
00635       AMESOS_CHK_ERR(SerialB_->ExtractView(&SerialBvalues_,&SerialXlda_ ));
00636     }
00637   }
00638   
00639   
00640   PrivateParakleteData_->common_ = rcp(new paraklete_common());
00641   
00642   const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
00643   assert (MpiComm != 0);
00644   
00645   MPI_Comm PK_Comm;
00646   //
00647   //  Create an MPI group with MaxProcesses_ processes
00648   //
00649   if ( MaxProcesses_ != Comm().NumProc()) {
00650     if(ParakleteComm_)  {
00651       MPI_Comm_free(&ParakleteComm_);
00652       ParakleteComm_ = 0 ; 
00653     }
00654     std::vector<int> ProcsInGroup(MaxProcesses_);
00655     IamInGroup_ = false; 
00656     for (int i = 0 ; i < MaxProcesses_ ; ++i) {
00657       ProcsInGroup[i] = i;
00658       if ( Comm().MyPID() == i ) IamInGroup_ = true; 
00659     }
00660     
00661     MPI_Group OrigGroup, ParakleteGroup;
00662     MPI_Comm_group(MpiComm->GetMpiComm(), &OrigGroup);
00663     MPI_Group_incl(OrigGroup, MaxProcesses_, &ProcsInGroup[0], &ParakleteGroup);
00664     MPI_Comm_create(MpiComm->GetMpiComm(), ParakleteGroup, &ParakleteComm_);
00665     PK_Comm = ParakleteComm_ ; 
00666   } else {
00667     IamInGroup_ = true; 
00668     PK_Comm = MpiComm->GetMpiComm() ;
00669   }
00670   
00671   paraklete_common& pk_common =  *PrivateParakleteData_->common_ ;
00672   cholmod_common *cm = &(pk_common.cm) ;
00673   amesos_cholmod_l_start (cm) ;
00674   PK_DEBUG_INIT ("pk", cm) ;
00675   pk_common.nproc = MaxProcesses_ ;
00676   pk_common.myid = Comm().MyPID() ; 
00677   //pk_common.mpicomm = PK_Comm ; 
00678   cm->print = 1 ;
00679   cm->precise = TRUE ;
00680   cm->error_handler = my_handler ;
00681   
00682   pk_common.tol_diag = 0.001 ;
00683   pk_common.tol_offdiag = 0.1 ;
00684   pk_common.growth = 2. ;
00685 
00686   
00687 
00688   // "overhead" time for the following two methods is considered here
00689   AMESOS_CHK_ERR( ExportToSerial() );
00690 
00691   AMESOS_CHK_ERR( ConvertToParakleteCRS(true) );
00692 
00693   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00694 
00695   // All this time if PARAKLETE time
00696   AMESOS_CHK_ERR( PerformSymbolicFactorization() );
00697 
00698   NumSymbolicFact_++;
00699 
00700   IsSymbolicFactorizationOK_ = true;
00701   
00702   return 0;
00703 }
00704 
00705 //=============================================================================
00706 int Amesos_Paraklete::NumericFactorization() 
00707 {
00708   if ( true || !TrustMe_ ) 
00709   { 
00710     IsNumericFactorizationOK_ = false;
00711     if (IsSymbolicFactorizationOK_ == false)
00712       AMESOS_CHK_ERR(SymbolicFactorization());
00713     
00714     ResetTimer(1); // "overhead" time
00715 
00716     Epetra_CrsMatrix *CrsMatrixA = dynamic_cast<Epetra_CrsMatrix *>(RowMatrixA_);
00717     if ( false && CrsMatrixA == 0 )   // hack to get around Bug #1502 
00718       AMESOS_CHK_ERR( CreateLocalMatrixAndExporters() ) ;
00719     assert( NumGlobalElements_ == RowMatrixA_->NumGlobalCols() );
00720     if ( AddZeroToDiag_ == 0 ) 
00721       assert( numentries_ == RowMatrixA_->NumGlobalNonzeros() );
00722 
00723     AMESOS_CHK_ERR( ExportToSerial() );
00724     
00725     if (  false && CrsMatrixA == 0 ) {  // continuation of hack to avoid bug #1502
00726       AMESOS_CHK_ERR( ConvertToParakleteCRS(true) );
00727     }  else {
00728       AMESOS_CHK_ERR( ConvertToParakleteCRS(false) );
00729     }
00730 
00731     OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00732   }
00733 
00734   // this time is all for PARAKLETE
00735   AMESOS_CHK_ERR( PerformNumericFactorization() );
00736 
00737   NumNumericFact_++;
00738 
00739   IsNumericFactorizationOK_ = true;
00740   
00741   return 0;
00742 }
00743 
00744 //=============================================================================
00745 int Amesos_Paraklete::Solve() 
00746 {
00747   Epetra_MultiVector* vecX=0;
00748   Epetra_MultiVector* vecB=0;
00749 
00750   if ( !TrustMe_ ) { 
00751 
00752     SerialB_ = Problem_->GetRHS() ;
00753     SerialX_ = Problem_->GetLHS() ;
00754 
00755     Epetra_MultiVector* OrigVecX ;
00756     Epetra_MultiVector* OrigVecB ;
00757 
00758     if (IsNumericFactorizationOK_ == false)
00759       AMESOS_CHK_ERR(NumericFactorization());
00760     
00761     ResetTimer(1);
00762 
00763     //
00764     //  Reindex the LHS and RHS 
00765     //
00766     OrigVecX = Problem_->GetLHS() ;
00767     OrigVecB = Problem_->GetRHS() ;
00768     
00769     if ( Reindex_ ) { 
00770 #ifdef HAVE_AMESOS_EPETRAEXT
00771       vecX = StdIndexDomain_->StandardizeIndex( OrigVecX ) ;
00772       vecB = StdIndexRange_->StandardizeIndex( OrigVecB ) ;
00773 #else
00774       AMESOS_CHK_ERR( -13 ) ; // Amesos_Paraklete can't handle non-standard indexing without EpetraExt 
00775 #endif
00776     } else {
00777       vecX = OrigVecX ;
00778       vecB = OrigVecB ;
00779     } 
00780     
00781     if ((vecX == 0) || (vecB == 0))
00782       AMESOS_CHK_ERR(-1); // something wrong in input
00783     
00784     //  Extract Serial versions of X and B
00785     
00786     ResetTimer(0);
00787 
00788     //  Copy B to the serial version of B
00789     //
00790     if (false && UseDataInPlace_ == 1) {
00791       SerialB_ = vecB;
00792       SerialX_ = vecX;
00793       NumVectors_ = Problem_->GetRHS()->NumVectors() ; 
00794     } else {
00795       assert (UseDataInPlace_ == 0);
00796       
00797       if( NumVectors_ != Problem_->GetRHS()->NumVectors() ) {
00798   NumVectors_ = Problem_->GetRHS()->NumVectors() ; 
00799   SerialXextract_ = rcp( new Epetra_MultiVector(*SerialMap_,NumVectors_));
00800   SerialBextract_ = rcp (new Epetra_MultiVector(*SerialMap_,NumVectors_));
00801       }
00802       if (NumVectors_ != vecB->NumVectors())
00803   AMESOS_CHK_ERR(-1); // internal error 
00804       
00805       ImportRangeToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecB->Map() ) );
00806       if ( SerialBextract_->Import(*vecB,*ImportRangeToSerial_,Insert) )
00807   AMESOS_CHK_ERR( -1 ) ; // internal error
00808       
00809       SerialB_ = &*SerialBextract_ ;
00810       SerialX_ = &*SerialXextract_ ;
00811     }
00812     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00813 
00814     //  Call PARAKLETE to perform the solve
00815     
00816     ResetTimer(0);
00817     if (MyPID_ == 0) {
00818       AMESOS_CHK_ERR(SerialB_->ExtractView(&SerialBvalues_,&SerialXlda_ ));
00819       AMESOS_CHK_ERR(SerialX_->ExtractView(&SerialXBvalues_,&SerialXlda_ ));
00820       if (SerialXlda_ != NumGlobalElements_)
00821   AMESOS_CHK_ERR(-1);
00822     }
00823 
00824     OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00825   }
00826   if ( MyPID_ == 0) {
00827     if ( NumVectors_ == 1 ) {
00828       for ( int i = 0 ; i < NumGlobalElements_ ; i++ ) 
00829   SerialXBvalues_[i] = SerialBvalues_[i] ;
00830     } else {
00831       SerialX_->Scale(1.0, *SerialB_ ) ;    // X = B (Klu overwrites B with X)
00832     }
00833   }
00834   if ( IamInGroup_ ) 
00835     for (int i = 0; i < NumVectors_ ; i++ ) { 
00836       amesos_paraklete_solve(  &*PrivateParakleteData_->LUnumeric_, &*PrivateParakleteData_->LUsymbolic_,
00837       &SerialXBvalues_[i*SerialXlda_],  &*PrivateParakleteData_->common_ );
00838     }
00839   SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
00840 
00841   //  Copy X back to the original vector
00842 
00843   ResetTimer(0);
00844   ResetTimer(1);
00845 
00846   if (UseDataInPlace_ == 0) {
00847     ImportDomainToSerial_ = rcp(new Epetra_Import ( *SerialMap_, vecX->Map() ) );
00848     vecX->Export( *SerialX_, *ImportDomainToSerial_, Insert ) ;
00849 
00850   } // otherwise we are already in place.
00851 
00852   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00853 
00854 #if 0
00855   //
00856   //  ComputeTrueResidual causes TestOptions to fail on my linux box //  Bug #1417
00857   if (ComputeTrueResidual_)
00858     ComputeTrueResidual(*SerialMatrix_, *vecX, *vecB, UseTranspose(), "Amesos_Paraklete");
00859 #endif
00860 
00861   if (ComputeVectorNorms_)
00862     ComputeVectorNorms(*vecX, *vecB, "Amesos_Paraklete");
00863 
00864   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00865 
00866   ++NumSolve_;
00867 
00868   return(0) ;
00869 }
00870 
00871 // ================================================ ====== ==== ==== == =
00872 
00873 void Amesos_Paraklete::PrintStatus() const
00874 {
00875 
00876   if (MyPID_) return;
00877 
00878   PrintLine();
00879 
00880   std::cout << "Amesos_Paraklete : Matrix has " << NumGlobalElements_ << " rows"
00881        << " and " << numentries_ << " nonzeros" << std::endl;
00882   std::cout << "Amesos_Paraklete : Nonzero elements per row = "
00883        << 1.0*numentries_/NumGlobalElements_ << std::endl;
00884   std::cout << "Amesos_Paraklete : Percentage of nonzero elements = "
00885        << 100.0*numentries_/(pow(double(NumGlobalElements_),double(2.0))) << std::endl;
00886   std::cout << "Amesos_Paraklete : Use transpose = " << UseTranspose_ << std::endl;
00887 
00888   PrintLine();
00889 
00890   return;
00891 
00892 }
00893 
00894 // ================================================ ====== ==== ==== == =
00895 
00896 void Amesos_Paraklete::PrintTiming() const
00897 {
00898   if (MyPID_) return;
00899 
00900   double ConTime = GetTime(MtxConvTime_);
00901   double MatTime = GetTime(MtxRedistTime_);
00902   double VecTime = GetTime(VecRedistTime_);
00903   double SymTime = GetTime(SymFactTime_);
00904   double NumTime = GetTime(NumFactTime_);
00905   double SolTime = GetTime(SolveTime_);
00906   double OveTime = GetTime(OverheadTime_);
00907 
00908   if (NumSymbolicFact_)
00909     SymTime /= NumSymbolicFact_;
00910 
00911   if (NumNumericFact_)
00912     NumTime /= NumNumericFact_;
00913 
00914   if (NumSolve_)
00915     SolTime /= NumSolve_;
00916 
00917   std::string p = "Amesos_Paraklete : ";
00918   PrintLine();
00919 
00920   std::cout << p << "Time to convert matrix to Paraklete format = "
00921        << ConTime << " (s)" << std::endl;
00922   std::cout << p << "Time to redistribute matrix = "
00923        << MatTime << " (s)" << std::endl;
00924   std::cout << p << "Time to redistribute vectors = "
00925        << VecTime << " (s)" << std::endl;
00926   std::cout << p << "Number of symbolic factorizations = "
00927        << NumSymbolicFact_ << std::endl;
00928   std::cout << p << "Time for sym fact = "
00929        << SymTime * NumSymbolicFact_ << " (s), avg = " << SymTime << " (s)" << std::endl;
00930   std::cout << p << "Number of numeric factorizations = "
00931        << NumNumericFact_ << std::endl;
00932   std::cout << p << "Time for num fact = "
00933        << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
00934   std::cout << p << "Number of solve phases = "
00935        << NumSolve_ << std::endl;
00936   std::cout << p << "Time for solve = "
00937        << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
00938 
00939   double tt = SymTime * NumSymbolicFact_ + NumTime * NumNumericFact_ + SolTime * NumSolve_;
00940   if (tt != 0)
00941   {
00942     std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
00943     std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
00944     std::cout << p << "(the above time does not include PARAKLETE time)" << std::endl;
00945     std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
00946   }
00947 
00948   PrintLine();
00949 
00950   return;
00951 }
00952 
00953 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines