39 #include "trilinos_klu_decl.h" 47 template<
class T,
class DeleteFunctor>
48 class DeallocFunctorDeleteWithCommon
51 DeallocFunctorDeleteWithCommon(
52 const RCP<trilinos_klu_common> &common
53 ,DeleteFunctor deleteFunctor
55 : common_(common), deleteFunctor_(deleteFunctor)
59 if(ptr) deleteFunctor_(&ptr,&*common_);
62 Teuchos::RCP<trilinos_klu_common> common_;
63 DeleteFunctor deleteFunctor_;
64 DeallocFunctorDeleteWithCommon();
67 template<
class T,
class DeleteFunctor>
68 DeallocFunctorDeleteWithCommon<T,DeleteFunctor>
69 deallocFunctorDeleteWithCommon(
70 const RCP<trilinos_klu_common> &common
71 ,DeleteFunctor deleteFunctor
74 return DeallocFunctorDeleteWithCommon<T,DeleteFunctor>(common,deleteFunctor);
85 Teuchos::RCP<trilinos_klu_common>
common_;
105 Teuchos::ParameterList ParamList ;
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 ;
219 int NumMyElements_ = 0 ;
243 #ifdef HAVE_AMESOS_EPETRAEXT 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 273 int indexBase = OriginalMatrixMap.
IndexBase();
278 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 280 long long indexBase = OriginalMatrixMap.
IndexBase64();
285 throw "Amesos_Klu::CreateLocalMatrixAndExporters: Unknown Global Indices";
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 ;
461 if (ParameterList.isParameter(
"TrustMe") )
462 TrustMe_ = ParameterList.get<
bool>(
"TrustMe" );
468 if (ParameterList.isSublist(
"Klu") ) {
469 Teuchos::ParameterList KluParams = ParameterList.sublist(
"Klu") ;
491 ,deallocFunctorDeleteWithCommon<trilinos_klu_symbolic>(
PrivateKluData_->common_,trilinos_klu_free_symbolic)
503 if ( !symbolic_ok ) {
532 bool factor_with_pivoting = true ;
544 int result = trilinos_klu_refactor (&
Ap[0],
Ai,
Aval,
548 const bool refactor_ok = result == 1 &&
PrivateKluData_->common_->status == TRILINOS_KLU_OK ;
559 factor_with_pivoting = false ;
564 if ( factor_with_pivoting ) {
573 rcpWithDealloc( trilinos_klu_factor(&
Ap[0],
Ai,
Aval,
575 deallocFunctorDeleteWithCommon<trilinos_klu_numeric>(
PrivateKluData_->common_,trilinos_klu_free_numeric)
588 if ( ! numeric_ok ) {
609 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints64() !=
610 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints64() ) {
701 if ( CrsMatrixA == 0 )
709 if ( CrsMatrixA == 0 ) {
742 #ifdef HAVE_AMESOS_EPETRAEXT 743 Teuchos::RCP<Epetra_MultiVector> vecX_rcp;
744 Teuchos::RCP<Epetra_MultiVector> vecB_rcp;
750 lose_this_ = (
int *) malloc( 300 ) ;
757 if ( lose_this_[0] == 12834 ) {
758 std::cout << __FILE__ <<
"::" << __LINE__ <<
" very unlikely to happen " << std::endl ;
783 #ifdef HAVE_AMESOS_EPETRAEXT 797 if ((vecX == 0) || (vecB == 0))
807 #ifdef HAVE_AMESOS_EPETRAEXT 808 if(vecX_rcp==Teuchos::null)
809 SerialX_ = Teuchos::rcp(vecX,
false);
813 if(vecB_rcp==Teuchos::null)
814 SerialB_ = Teuchos::rcp(vecB,
false);
818 SerialB_ = Teuchos::rcp(vecB,
false);
819 SerialX_ = Teuchos::rcp(vecX,
false);
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;
const Epetra_Map & RowMap() const
Teuchos::RCP< Epetra_Import > ImportRangeToSerial_
bool StorageOptimized() const
int NumSymbolicFact_
Number of symbolic factorization phases.
int PerformNumericFactorization()
long long numentries_
Number of non-zero entries in Problem_->GetOperator()
double GetTime(const std::string what) const
Gets the cumulative time using the string.
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...
virtual const Epetra_Map & RowMatrixRowMap() const=0
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.
bool SameAs(const Epetra_BlockMap &Map) const
int Solve()
Solves A X = B (or AT x = B)
Teuchos::RCP< Epetra_Map > SerialMap_
Points to a Serial Map (unused if UseDataInPlace_ == 1 )
double rcond_threshold_
If error is greater than this value, perform symbolic and numeric factorization with full partial piv...
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Points to a Serial Copy of A (unused if UseDataInPlace_==1)
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
Teuchos::RCP< Epetra_Import > ImportDomainToSerial_
virtual const Epetra_Map & OperatorDomainMap() const=0
Epetra_CrsMatrix * CrsMatrixA_
Operator converted to a CrsMatrix.
int MtxRedistTime_
Quick access ids for the individual timings.
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Teuchos::RCP< trilinos_klu_common > common_
Teuchos::RCP< Epetra_MultiVector > SerialX_
long long NumGlobalElements64() const
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.
virtual int MyPID() const=0
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.
virtual int MaxNumEntries() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
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.
bool GlobalIndicesInt() const
int NumMyElements() const
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.
Teuchos::RCP< trilinos_klu_symbolic > Symbolic_
void PrintTiming() const
Prints timing information.
virtual long long NumGlobalCols64() const=0
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.
virtual int Broadcast(double *MyVals, int Count, int Root) const=0
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.
virtual long long NumGlobalNonzeros64() const=0
Teuchos::RCP< Amesos_Klu_Pimpl > PrivateKluData_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
long long IndexBase64() const
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
int UseDataInPlace_
1 if Problem_->GetOperator() is stored entirely on process 0
Epetra_MultiVector * GetRHS() const
virtual int NumProc() const=0
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
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.
Epetra_MultiVector * GetLHS() const
Epetra_RowMatrix * GetMatrix() const
bool MatrixShapeOK() const
Returns true if KLU can handle this matrix shape.
int verbose_
Toggles the output level.
std::vector< double > VecAval
Teuchos::RCP< trilinos_klu_numeric > Numeric_
bool GlobalIndicesLongLong() const
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_
Epetra_Operator * GetOperator() const
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
void PrintLine() const
Prints line on std::cout.
virtual long long NumGlobalRows64() const=0
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
double AddToDiag_
Add this value to the diagonal.