30 #include "Epetra_Map.h" 31 #include "Epetra_Import.h" 32 #include "Epetra_CrsMatrix.h" 33 #include "Epetra_Vector.h" 34 #include "Epetra_Util.h" 47 template<
class T,
class DeleteFunctor>
48 class DeallocFunctorDeleteWithCommon
51 DeallocFunctorDeleteWithCommon(
53 ,DeleteFunctor deleteFunctor
55 : common_(common), deleteFunctor_(deleteFunctor)
59 if(ptr) deleteFunctor_(&ptr,&*common_);
63 DeleteFunctor deleteFunctor_;
64 DeallocFunctorDeleteWithCommon();
67 template<
class T,
class DeleteFunctor>
68 DeallocFunctorDeleteWithCommon<T,DeleteFunctor>
69 deallocFunctorDeleteWithCommon(
71 ,DeleteFunctor deleteFunctor
74 return DeallocFunctorDeleteWithCommon<T,DeleteFunctor>(common,deleteFunctor);
123 std::cerr <<
" The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
129 std::cerr <<
" The number of non zero entries in the matrix has changed since the last call to SymbolicFactorization(). " ;
145 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 147 for (
int i = 0 ; i <
SerialMap_->NumGlobalElements(); i++ ) {
154 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 156 for (
long long i = 0 ; i <
SerialMap_->NumGlobalElements64(); i++ ) {
163 throw "Amesos_Klu::ExportToSerial: ERROR, GlobalIndices type unknown.";
171 std::cerr <<
" Amesos_Klu cannot handle this matrix. " ;
173 std::cerr <<
"Unknown error" << std::endl ;
176 std::cerr <<
" Try setting the Reindex parameter to true. " << std::endl ;
177 #ifndef HAVE_AMESOS_EPETRAEXT 178 std::cerr <<
" You will need to rebuild the Amesos library with the EpetraExt library to use the reindex feature" << std::endl ;
179 std::cerr <<
" To rebuild Amesos with EpetraExt, add --enable-epetraext to your configure invocation" << std::endl ;
203 const Epetra_Map &OriginalMatrixMap =
RowMatrixA_->RowMatrixRowMap() ;
204 const Epetra_Map &OriginalDomainMap =
206 GetProblem()->GetOperator()->OperatorDomainMap();
207 const Epetra_Map &OriginalRangeMap =
209 GetProblem()->GetOperator()->OperatorRangeMap();
219 int NumMyElements_ = 0 ;
227 OriginalMatrixMap.NumGlobalElements64() )?1:0;
228 if ( ! OriginalRangeMap.SameAs( OriginalMatrixMap ) )
UseDataInPlace_ = 0 ;
229 if ( ! OriginalDomainMap.SameAs( OriginalMatrixMap ) )
UseDataInPlace_ = 0 ;
243 #ifdef HAVE_AMESOS_EPETRAEXT 244 const Epetra_Map& OriginalMap =
CrsMatrixA_->RowMap();
252 std::cerr <<
"Amesos_Klu requires EpetraExt to reindex matrices." << std::endl
253 <<
" Please rebuild with the EpetraExt library by adding --enable-epetraext to your configure invocation" << std::endl ;
271 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 272 if(OriginalMatrixMap.GlobalIndicesInt()) {
273 int indexBase = OriginalMatrixMap.IndexBase();
278 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 279 if(OriginalMatrixMap.GlobalIndicesLongLong()) {
280 long long indexBase = OriginalMatrixMap.IndexBase64();
285 throw "Amesos_Klu::CreateLocalMatrixAndExporters: Unknown Global Indices";
354 Epetra_CrsMatrix *CrsMatrix =
dynamic_cast<Epetra_CrsMatrix *
>(
SerialMatrix_);
355 bool StorageOptimized = ( CrsMatrix != 0 && CrsMatrix->StorageOptimized() );
357 if (
AddToDiag_ != 0.0 ) StorageOptimized = false ;
361 if ( ! StorageOptimized ) {
371 int NumEntriesThisRow;
373 if( StorageOptimized ) {
378 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
385 Ap[MyRow+1] =
Ap[MyRow] + NumEntriesThisRow ;
394 if ( firsttime && CrsMatrix == 0 ) {
400 if ( CrsMatrix != 0 ) {
402 ExtractMyRowView( MyRow, NumEntriesThisRow, RowValues,
408 ExtractMyRowCopy( MyRow, MaxNumEntries_,
418 Ap[MyRow] = Ai_index ;
419 for (
int j = 0; j < NumEntriesThisRow; j++ ) {
420 VecAi[Ai_index] = ColIndices[j] ;
422 VecAval[Ai_index] = RowValues[j] ;
423 if (ColIndices[j] == MyRow) {
429 for (
int j = 0; j < NumEntriesThisRow; j++ ) {
430 VecAval[Ai_index] = RowValues[j] ;
431 if (ColIndices[j] == MyRow) {
438 Ap[MyRow] = Ai_index ;
501 Comm().Broadcast(&symbolic_ok, 1, 0);
503 if ( !symbolic_ok ) {
532 bool factor_with_pivoting = true ;
559 factor_with_pivoting = false ;
564 if ( factor_with_pivoting ) {
586 Comm().Broadcast(&numeric_ok, 1, 0);
588 if ( ! numeric_ok ) {
609 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints64() !=
610 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints64() ) {
700 Epetra_CrsMatrix *CrsMatrixA =
dynamic_cast<Epetra_CrsMatrix *
>(
RowMatrixA_);
701 if ( CrsMatrixA == 0 )
709 if ( CrsMatrixA == 0 ) {
739 Epetra_MultiVector* vecX = 0 ;
740 Epetra_MultiVector* vecB = 0 ;
742 #ifdef HAVE_AMESOS_EPETRAEXT 750 lose_this_ = (
int *) malloc( 300 ) ;
757 if ( lose_this_[0] == 12834 ) {
758 std::cout << __FILE__ <<
"::" << __LINE__ <<
" very unlikely to happen " << std::endl ;
768 Epetra_MultiVector* OrigVecX ;
769 Epetra_MultiVector* OrigVecB ;
783 #ifdef HAVE_AMESOS_EPETRAEXT 797 if ((vecX == 0) || (vecB == 0))
807 #ifdef HAVE_AMESOS_EPETRAEXT 808 if(vecX_rcp==Teuchos::null)
813 if(vecB_rcp==Teuchos::null)
835 Epetra_Import *UseImport;
895 Epetra_Import *UseImport;
899 vecX->Export( *
SerialX_, *UseImport, Insert ) ;
934 <<
" and " <<
numentries_ <<
" nonzeros" << std::endl;
935 std::cout <<
"Amesos_Klu : Nonzero elements per row = " 937 std::cout <<
"Amesos_Klu : Percentage of nonzero elements = " 939 std::cout <<
"Amesos_Klu : Use transpose = " <<
UseTranspose_ << std::endl;
970 std::string p =
"Amesos_Klu : ";
973 std::cout << p <<
"Time to convert matrix to Klu format = " 974 << ConTime <<
" (s)" << std::endl;
975 std::cout << p <<
"Time to redistribute matrix = " 976 << MatTime <<
" (s)" << std::endl;
977 std::cout << p <<
"Time to redistribute vectors = " 978 << VecTime <<
" (s)" << std::endl;
979 std::cout << p <<
"Number of symbolic factorizations = " 981 std::cout << p <<
"Time for sym fact = " 982 << SymTime *
NumSymbolicFact_ <<
" (s), avg = " << SymTime <<
" (s)" << std::endl;
983 std::cout << p <<
"Number of numeric factorizations = " 985 std::cout << p <<
"Time for num fact = " 986 << NumTime *
NumNumericFact_ <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
987 std::cout << p <<
"Number of solve phases = " 989 std::cout << p <<
"Time for solve = " 990 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
995 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
996 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
997 std::cout << p <<
"(the above time does not include KLU time)" << std::endl;
998 std::cout << p <<
"Amesos interface time / total time = " << OveTime / tt << std::endl;
Teuchos::RCP< Epetra_Import > ImportRangeToSerial_
int NumSymbolicFact_
Number of symbolic factorization phases.
int PerformNumericFactorization()
long long numentries_
Number of non-zero entries in Problem_->GetOperator()
int amesos_klu_refactor(int Ap [], int Ai [], double Ax [], klu_symbolic *Symbolic, klu_numeric *Numeric, klu_common *Common)
double GetTime(const std::string what) const
Gets the cumulative time using the string.
klu_symbolic * amesos_klu_analyze(int n, int Ap [], int Ai [], klu_common *Common)
bool TrustMe_
If true, no checks are made and the matrix is assume to be distributed.
std::vector< int > Ap
Ap, Ai, Aval form the compressed row storage used by Klu Ai and Aval can point directly into a matrix...
Teuchos::RCP< Amesos_StandardIndex > StdIndexRange_
Epetra_RowMatrix * RowMatrixA_
Operator converted to a RowMatrix.
long long NumGlobalElements_
Number of rows and columns in the Problem_->GetOperator()
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
T & get(ParameterList &l, const std::string &name)
int Solve()
Solves A X = B (or AT x = B)
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if UseDataInPlace_ == 1 )
struct klu_common_struct klu_common
klu_numeric * amesos_klu_factor(int Ap [], int Ai [], double Ax [], klu_symbolic *Symbolic, klu_common *Common)
double rcond_threshold_
If error is greater than this value, perform symbolic and numeric factorization with full partial piv...
bool isSublist(const std::string &name) const
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Points to a Serial Copy of A (unused if UseDataInPlace_==1)
int amesos_klu_free_numeric(klu_numeric **Numeric, klu_common *Common)
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
RCP< ParameterList > sublist(const RCP< ParameterList > ¶mList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< Epetra_Import > ImportDomainToSerial_
Epetra_CrsMatrix * CrsMatrixA_
Operator converted to a CrsMatrix.
Teuchos::RCP< klu_numeric > Numeric_
int MtxRedistTime_
Quick access ids for the individual timings.
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Teuchos::RCP< Epetra_MultiVector > SerialX_
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RCP< Amesos_StandardIndex > StdIndexDomain_
Epetra_RowMatrix * SerialMatrix_
Points to a Serial Copy of A.
int NumNumericFact_
Number of numeric factorization phases.
double * SerialXBvalues_
Pointer to the actual values in the serial version of X and B.
std::vector< double > RowValuesV_
Only used for RowMatrices to extract copies.
int CreateLocalMatrixAndExporters()
int NumSolve_
Number of solves.
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
int amesos_klu_solve(klu_symbolic *Symbolic, klu_numeric *Numeric, int ldim, int nrhs, double B [], klu_common *Common)
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
std::vector< int > ColIndicesV_
Only used for RowMatrices to extract copies.
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer to process 0.
int ConvertToKluCRS(bool firsttime)
int PerformSymbolicFactorization()
void ComputeTrueResidual(const Epetra_RowMatrix &Matrix, const Epetra_MultiVector &X, const Epetra_MultiVector &B, const bool UseTranspose, const std::string prefix) const
Computes the true residual, B - Matrix * X, and prints the results.
bool PrintTiming_
If true, prints timing information in the destructor.
Epetra_RowMatrix * StdIndexMatrix_
Points to a Contiguous Copy of A.
void PrintStatus() const
Prints information about the factorization and solution phases.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool UseTranspose_
If true, the transpose of A is used.
~Amesos_Klu(void)
Amesos_Klu Destructor.
bool PrintStatus_
If true, print additional information in the destructor.
void PrintTiming() const
Prints timing information.
Teuchos::RCP< klu_common > common_
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
bool UseTranspose() const
Returns the current UseTranspose setting.
int amesos_klu_rcond(klu_symbolic *Symbolic, klu_numeric *Numeric, klu_common *Common)
int NumVectors_
Number of vectors in RHS and LHS.
bool Reindex_
If true, the Amesos class should reindex the matrix to standard indexing (i.e.
Teuchos::RCP< klu_symbolic > Symbolic_
Teuchos::RCP< Amesos_Klu_Pimpl > PrivateKluData_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
int amesos_klu_defaults(klu_common *Common)
int UseDataInPlace_
1 if Problem_->GetOperator() is stored entirely on process 0
int amesos_klu_free_symbolic(klu_symbolic **Symbolic, klu_common *Common)
int amesos_klu_tsolve(klu_symbolic *Symbolic, klu_numeric *Numeric, int ldim, int nrhs, double B [], klu_common *Common)
Teuchos::RCP< Epetra_MultiVector > SerialBextract_
const int NumericallySingularMatrixError
void ComputeVectorNorms(const Epetra_MultiVector &X, const Epetra_MultiVector &B, const std::string prefix) const
Computes the norms of X and B and print the results.
const int StructurallySingularMatrixError
bool isParameter(const std::string &name) const
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Teuchos::RCP< Epetra_MultiVector > SerialB_
Serial versions of the LHS and RHS (may point to the original vector if serial)
Amesos_Klu(const Epetra_LinearProblem &LinearProblem)
Amesos_Klu Constructor.
bool MatrixShapeOK() const
Returns true if KLU can handle this matrix shape.
int verbose_
Toggles the output level.
std::vector< double > VecAval
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
Teuchos::RCP< Epetra_MultiVector > SerialXextract_
Serial versions of the LHS and RHS (if necessary)
Teuchos::RCP< Amesos_StandardIndex > StdIndex_
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void PrintLine() const
Prints line on std::cout.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
double AddToDiag_
Add this value to the diagonal.