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