30 #include "Epetra_Map.h" 31 #include "Epetra_Import.h" 32 #include "Epetra_RowMatrix.h" 33 #include "Epetra_CrsMatrix.h" 34 #include "Epetra_Vector.h" 35 #include "Epetra_Util.h" 36 #include "Epetra_Time.h" 37 #include "Epetra_Comm.h" 38 #include "Epetra_LinearProblem.h" 42 #ifdef AMESOS_SUPERLU_PRE_JULY2005 43 # include "dsp_defs.h" 45 # include "slu_ddefs.h" 61 A.Store =
B.Store =
X.Store =
L.Store =
U.Store =
NULL;
76 NumGlobalNonzeros_(-1),
78 FactorizationOK_(false),
79 FactorizationDone_(false),
88 SerialCrsMatrixA_(
Teuchos::null),
100 dCreate_Dense_Matrix( &(
data_->
X),
105 SLU::SLU_DN, SLU::SLU_D, SLU::SLU_GE);
107 dCreate_Dense_Matrix( &(
data_->
B),
112 SLU::SLU_DN, SLU::SLU_D, SLU::SLU_GE);
121 Destroy_SuperMatrix_Store(&
data_->
B);
122 Destroy_SuperMatrix_Store(&
data_->
X);
126 Destroy_SuperMatrix_Store(&
data_->
A);
127 Destroy_SuperNode_Matrix(&
data_->
L);
128 Destroy_CompCol_Matrix(&
data_->
U);
155 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints64() !=
156 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints64())
174 const Epetra_Map &OriginalMap =
RowMatrixA_->RowMatrixRowMap() ;
183 int NumMyElements_ = 0;
192 if (
Comm().NumProc() == 1)
198 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) 199 if(
RowMatrixA_->RowMatrixRowMap().GlobalIndicesInt())
203 #if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 204 if(
RowMatrixA_->RowMatrixRowMap().GlobalIndicesLongLong())
208 throw "Amesos_Superlu::ConvertToSerial: Global indices unknown.";
229 std::cout << __FILE__ <<
"::" << __LINE__
230 <<
" this_res = " << this_res
296 Ap_[MyRow] = Ai_index;
298 ierr =
SerialMatrix_->ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
301 Ai_index += NzThisRow;
309 Destroy_SuperMatrix_Store(&
data_->
A);
310 Destroy_SuperNode_Matrix(&
data_->
L);
311 Destroy_CompCol_Matrix(&
data_->
U);
317 &
Ai_[0], &
Ap_[0], SLU::SLU_NR, SLU::SLU_D, SLU::SLU_GE );
358 std::vector<int> ColIndicesV_;
359 std::vector<double> RowValuesV_;
362 Epetra_CrsMatrix *SuperluCrs =
dynamic_cast<Epetra_CrsMatrix *
>(
SerialMatrix_);
363 if ( SuperluCrs == 0 ) {
364 ColIndicesV_.resize(MaxNumEntries_);
365 RowValuesV_.resize(MaxNumEntries_);
368 if ( SuperluCrs != 0 ) {
369 AMESOS_CHK_ERR(SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues, ColIndices));
374 RowValues = &RowValuesV_[0];
375 ColIndices = &ColIndicesV_[0];
378 if (
Ap_[MyRow] != Ai_index)
381 for (
int j = 0; j < NzThisRow; j++) {
382 assert(
Ai_[Ai_index] == ColIndices[j]) ;
383 Aval_[Ai_index] = RowValues[j] ;
393 Destroy_SuperMatrix_Store(&
data_->
A);
394 Destroy_SuperNode_Matrix(&
data_->
L);
395 Destroy_CompCol_Matrix(&
data_->
U);
399 &
Ai_[0], &
Ap_[0], SLU::SLU_NR, SLU::SLU_D, SLU::SLU_GE );
430 set_default_options( &SLUopt ) ;
437 SLUopt.Fact = SLU::DOFACT;
447 assert( SLUopt.Trans == NOTRANS ) ;
449 SLUopt.Trans = TRANS ;
454 std::cout <<
" SLUopt.ColPerm = " << SLUopt.ColPerm << std::endl ;
455 std::cout <<
" SLUopt.Equil = " << SLUopt.Equil << std::endl ;
456 std::cout <<
" SLUopt.Fact = " << SLUopt.Fact << std::endl ;
457 std::cout <<
" SLUopt.IterRefine = " << SLUopt.IterRefine << std::endl ;
458 std::cout <<
" data_->A.Stype = " <<
data_->
A.Stype
459 <<
" SLU_NC = " << SLU_NC
460 <<
" SLU_NR = " << SLU_NR
462 std::cout <<
" SLUopt.ColPerm = " << SLUopt.ColPerm << std::endl ;
467 SLU::DNformat* Bstore = (SLU::DNformat *) (
data_->
B.Store) ;
472 SLU::DNformat* Xstore = (SLU::DNformat *) (
data_->
X.Store) ;
476 SLU::SuperLUStat_t SLU_stat ;
477 SLU::StatInit( &SLU_stat ) ;
478 assert( SLUopt.Fact == SLU::DOFACT);
479 dgssvx( &(SLUopt), &(
data_->
A),
485 SLU::StatFree( &SLU_stat ) ;
507 Epetra_MultiVector* vecX =
Problem_->GetLHS();
508 Epetra_MultiVector* vecB =
Problem_->GetRHS();
511 if (vecX == 0 || vecB == 0)
514 int nrhs = vecX->NumVectors();
515 if (nrhs != vecB->NumVectors())
521 Epetra_MultiVector* SerialB = 0;
522 Epetra_MultiVector* SerialX = 0;
524 double *SerialXvalues ;
525 double *SerialBvalues ;
527 if (
Comm().NumProc() == 1)
536 SerialX =
new Epetra_MultiVector(
SerialMap(), nrhs);
537 SerialB =
new Epetra_MultiVector(
SerialMap(), nrhs);
558 ierr = SerialX->ExtractView(&SerialXvalues, &SerialXlda);
562 ierr = SerialB->ExtractView(&SerialBvalues, &SerialBlda);
566 SLU::SuperMatrix& dataX = (
data_->
X) ;
570 SLU::DNformat* Xstore = (SLU::DNformat *) (
data_->
X.Store) ;
571 Xstore->lda = SerialXlda;
572 Xstore->nzval = &SerialXvalues[0];
574 SLU::SuperMatrix& dataB = (
data_->
B) ;
578 SLU::DNformat* Bstore = (SLU::DNformat *) (
data_->
B.Store) ;
579 Bstore->lda = SerialBlda;
580 Bstore->nzval = &SerialBvalues[0];
584 char fact, trans, refact;
599 SLU::SuperLUStat_t SLU_stat ;
600 SLU::StatInit( &SLU_stat ) ;
604 SLUopt.Trans = SLU::TRANS;
606 SLUopt.Trans = SLU::NOTRANS;
610 dgssvx( &(SLUopt), &(
data_->
A),
617 StatFree( &SLU_stat ) ;
627 if (
Comm().NumProc() != 1)
650 if (
Comm().NumProc() != 1)
651 Comm().Broadcast(&Ierr, 1, 0);
660 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
663 std::string p =
"Amesos_Superlu : ";
666 long long n =
GetProblem()->GetMatrix()->NumGlobalRows64();
667 long long nnz =
GetProblem()->GetMatrix()->NumGlobalNonzeros64();
669 std::cout << p <<
"Matrix has " << n <<
" rows" 670 <<
" and " <<
nnz <<
" nonzeros" << std::endl;
671 std::cout << p <<
"Nonzero elements per row = " 672 << 1.0 *
nnz / n << std::endl;
673 std::cout << p <<
"Percentage of nonzero elements = " 674 << 100.0 *
nnz /(pow(
double(n),
double(2.0))) << std::endl;
675 std::cout << p <<
"Use transpose = " <<
UseTranspose_ << std::endl;
685 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
701 std::string p =
"Amesos_Superlu : ";
704 std::cout << p <<
"Time to convert matrix to SuperLU format = " << ConTime <<
" (s)" << std::endl;
705 std::cout << p <<
"Time to redistribute matrix = " << MatTime <<
" (s)" << std::endl;
706 std::cout << p <<
"Time to redistribute vectors = " << VecTime <<
" (s)" << std::endl;
707 std::cout << p <<
"Number of numeric factorizations = " <<
NumNumericFact_ << std::endl;
708 std::cout << p <<
"Time for num fact = " 710 << NumTime <<
" (s)" << std::endl;
711 std::cout << p <<
"Number of solve phases = " <<
NumSolve_ << std::endl;
712 std::cout << p <<
"Time for solve = " 713 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
717 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
718 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
719 std::cout << p <<
"(the above time does not include SuperLU time)" << std::endl;
720 std::cout << p <<
"Amesos interface time / total time = " << OveTime / tt << std::endl;
double GetTime(const std::string what) const
Gets the cumulative time using the string.
std::vector< int > etree_
Teuchos::RCP< Epetra_Map > SerialMap_
Contains a map with all elements assigned to processor 0.
std::vector< double > Aval_
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
UF_long CHOLMOD() nnz(cholmod_sparse *A, cholmod_common *Common)
Epetra_RowMatrix * SerialMatrix_
For parallel runs, stores the matrix defined on SerialMap_.
SLU::fact_t refactor_option
SLU::superlu_options_t SLU_options
int Solve()
Solves A X = B (or AT x = B)
void PrintTiming() const
Prints timing information.
bool MatrixShapeOK() const
Returns true if the solver can handle this matrix shape.
std::vector< int > perm_c_
std::vector< int > Ap_
stores the matrix in SuperLU format.
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
int Factor()
Factors the matrix, no previous factorization available.
int MtxConvTime_
Quick access pointer to internal timing data.
bool UseTranspose_
If true, solve the linear system with the transpose of the matrix.
long long NumGlobalNonzeros_
Global number of nonzeros in the matrix.
int NumNumericFact_
Number of numeric factorization phases.
SLU::mem_usage_t mem_usage
int ReFactor()
Re-factors the matrix.
SLU::factor_param_t iparam
int NumSolve_
Number of solves.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
Amesos_Superlu(const Epetra_LinearProblem &LinearProblem)
Amesos_Superlu Constructor.
const Epetra_Import & ImportToSerial() const
Returns a reference to the importer.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
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.
std::vector< double > ferr_
bool UseTranspose() const
Returns the current UseTranspose setting.
bool PrintStatus_
If true, print additional information in the destructor.
#define AMESOS_RETURN(amesos_err)
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
void PrintStatus() const
Prints status information.
Epetra_RowMatrix * RowMatrixA_
Pointer to the linear system matrix.
double * DummyArray
stores the matrix in SuperLU format.
int ConvertToSerial()
Sets up the matrix on processor 0.
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrixA_
Contains a matrix with all rows assigned to processor 0.
Teuchos::RCP< Epetra_Import > ImportToSerial_
Importer from distributed to SerialMap_.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
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.
SLUData * data_
Main structure for SuperLU.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
~Amesos_Superlu()
Amesos_Superlu Destructor.
std::vector< double > berr_
std::vector< int > perm_r_
bool FactorizationOK_
If true, the factorization has been successfully computed.
int iam_
Process number (i.e. Comm().MyPID()
const Epetra_Map & SerialMap() const
Returns a reference to the serial map.
long long NumGlobalRows_
Global size of the matrix.
void PrintLine() const
Prints line on std::cout.
std::vector< int > Ai_
stores the matrix in SuperLU format.
const Epetra_LinearProblem * Problem_
Pointer to the user's defined linear problem.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().