33 #include "Epetra_Map.h" 34 #include "Epetra_Import.h" 35 #include "Epetra_CrsMatrix.h" 36 #include "Epetra_Util.h" 38 #include "superlu_ddefs.h" 39 #include "supermatrix.h" 43 #if SUPERLU_DIST_MAJOR_VERSION > 6 || (SUPERLU_DIST_MAJOR_VERSION == 6 && SUPERLU_DIST_MINOR_VERSION > 2) 44 #define ScalePermstruct_t dScalePermstruct_t 45 #define LUstruct_t dLUstruct_t 46 #define SOLVEstruct_t dSOLVEstruct_t 47 #define ScalePermstructInit dScalePermstructInit 48 #define ScalePermstructFree dScalePermstructFree 49 #define Destroy_LU dDestroy_LU 50 #define LUstructFree dLUstructFree 51 #define LUstructInit dLUstructInit 67 #if SUPERLU_DIST_MAJOR_VERSION > 4 85 for ( i = 1; i*i <= NumProcs; i++ )
88 for ( NumProcRows = i-1 ; done == false ; ) {
89 int NumCols = NumProcs / NumProcRows ;
90 if ( NumProcRows * NumCols == NumProcs )
103 if (nprow < 1 ) nprow = 1;
104 npcol = MaxProcesses / nprow;
106 if( nprow <=0 || npcol <= 0 || MaxProcesses <=0 ) {
107 std::cerr <<
"Amesos_Superludist: wrong value for MaxProcs (" 108 << MaxProcesses <<
"), or nprow (" << nprow
109 <<
") or npcol (" << npcol <<
")" << std::endl;
121 FactorizationDone_(0),
122 FactorizationOK_(false),
126 PrintNonzeros_(false),
131 IterRefine_(
"NOT SET"),
132 ReplaceTinyPivot_(true),
149 #if (SUPERLU_DIST_MAJOR_VERSION > 5) || ( SUPERLU_DIST_MAJOR_VERSION == 5 && SUPERLU_DIST_MINOR_VERSION > 3) 165 Teuchos::ParameterList ParamList;
198 if (ParameterList.isParameter(
"Redistribute"))
203 if (ParameterList.isSublist(
"Superludist") )
205 const Teuchos::ParameterList& SuperludistParams =
206 ParameterList.sublist(
"Superludist") ;
208 if( SuperludistParams.isParameter(
"ReuseSymbolic") )
210 std::string FactOption =
"NotSet";
212 if( SuperludistParams.isParameter(
"Fact") )
213 FactOption = SuperludistParams.get<std::string>(
"Fact");
215 if( FactOption ==
"SamePattern_SameRowPerm" )
PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm;
217 else if ( FactOption !=
"NotSet" )
220 if (SuperludistParams.isParameter(
"Equil"))
221 Equil_ = SuperludistParams.get<
bool>(
"Equil");
223 if (SuperludistParams.isParameter(
"ColPerm"))
224 ColPerm_ = SuperludistParams.get<std::string>(
"ColPerm");
227 if( SuperludistParams.isParameter(
"perm_c"))
228 perm_c_ = SuperludistParams.get<
int*>(
"perm_c");
230 if (SuperludistParams.isParameter(
"RowPerm"))
231 RowPerm_ = SuperludistParams.get<std::string>(
"RowPerm");
233 if (SuperludistParams.isParameter(
"perm_r"))
234 perm_r_ = SuperludistParams.get<
int*>(
"perm_r");
237 if (SuperludistParams.isParameter(
"IterRefine"))
238 IterRefine_ = SuperludistParams.get<std::string>(
"IterRefine");
240 if (SuperludistParams.isParameter(
"ReplaceTinyPivot"))
243 if (SuperludistParams.isParameter(
"PrintNonzeros"))
265 int iam =
Comm().MyPID();
272 int MyFirstElement = iam * m_per_p + EPETRA_MIN( iam, remainder );
273 int MyFirstNonElement = (iam+1) * m_per_p + EPETRA_MIN( iam+1, remainder );
274 int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
277 NumExpectedElements = 0;
282 const Epetra_Map &OriginalMap =
RowMatrixA_->RowMatrixRowMap();
300 for (
int i = 0 ; i <
UniformMap_->NumGlobalElements(); i++ )
333 if (
Comm().NumProc() == 1)
343 if (
Comm().NumProc() == 1)
362 int MyActualFirstElement =
UniformMatrix().RowMatrixRowMap().MinMyGID() ;
365 Ap_.resize( NumMyElements+1 );
366 Ai_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
367 Aval_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
376 std::vector<double> RowValuesV_(MaxNumEntries_);
377 std::vector<int> ColIndicesV_(MaxNumEntries_);
381 const Epetra_CrsMatrix *SuperluCrs =
dynamic_cast<const Epetra_CrsMatrix *
>(&
UniformMatrix());
385 for (MyRow = 0; MyRow < NumMyElements ; MyRow++)
389 ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
394 ierr =
UniformMatrix().ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
395 &RowValuesV_[0], &ColIndicesV_[0]);
396 RowValues = &RowValuesV_[0];
397 ColIndices = &ColIndicesV_[0];
404 for (
int i = 0 ; i < NzThisRow ; ++i) {
405 if (ColIndices[i] == MyRow) {
412 Ap_[MyRow] = Ai_index ;
413 for (
int j = 0; j < NzThisRow; j++ ) {
415 Aval_[Ai_index] = RowValues[j] ;
419 assert( NumMyElements == MyRow );
420 Ap_[ NumMyElements ] = Ai_index ;
427 const Epetra_MpiComm & comm1 =
dynamic_cast<const Epetra_MpiComm &
> (
Comm());
458 nnz_loc, NumMyElements, MyActualFirstElement,
460 SLU_NR_loc, SLU_D, SLU_GE );
466 #if SUPERLU_DIST_MAJOR_VERSION > 6 || (SUPERLU_DIST_MAJOR_VERSION == 6 && SUPERLU_DIST_MINOR_VERSION > 2) 469 #ifdef HAVE_SUPERLUDIST_LUSTRUCTINIT_2ARG 495 #if (SUPERLU_DIST_MAJOR_VERSION > 5) || ( SUPERLU_DIST_MAJOR_VERSION == 5 && SUPERLU_DIST_MINOR_VERSION > 3) 511 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE 520 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE 605 const Epetra_CrsMatrix *SuperluCrs =
dynamic_cast<const Epetra_CrsMatrix *
>(&
UniformMatrix());
611 std::vector<int> ColIndicesV_(MaxNumEntries_);
612 std::vector<double> RowValuesV_(MaxNumEntries_);
619 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
620 if ( SuperluCrs != 0 ) {
621 ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
625 ierr =
UniformMatrix().ExtractMyRowCopy( MyRow, MaxNumEntries_,
626 NzThisRow, &RowValuesV_[0],
628 RowValues = &RowValuesV_[0];
629 ColIndices = &ColIndicesV_[0];
635 for (
int j = 0; j < NzThisRow; j++ ) {
638 Aval_[Ai_index] = RowValues[j] ;
680 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
681 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints())
738 Epetra_MultiVector* vecX =
Problem_->GetLHS();
739 Epetra_MultiVector* vecB =
Problem_->GetRHS();
741 if (vecX == 0 || vecB == 0)
744 int nrhs = vecX->NumVectors() ;
745 if (vecB->NumVectors() != nrhs)
751 RCP<Epetra_MultiVector> vec_uni;
755 vec_uni = Teuchos::rcp(
new Epetra_MultiVector(*
UniformMap_, nrhs));
757 vec_uni->Import(*vecB,
Importer(), Insert);
762 vecX->Update(1.0, *vecB, 0.0);
763 vec_uni = Teuchos::rcp(vecX,
false);
777 std::vector<double>berr(nrhs);
801 vecX->Export(*vec_uni,
Importer(), Insert);
807 "Amesos_Superludist");
821 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
824 std::string p =
"Amesos_Superludist : ";
830 <<
" and " << NNZ <<
" nonzeros" << std::endl;
831 std::cout << p <<
"Nonzero elements per row = " 833 std::cout << p <<
"Percentage of nonzero elements = " 835 std::cout << p <<
"Use transpose = " <<
UseTranspose() << std::endl;
836 std::cout << p <<
"Redistribute = " <<
Redistribute_ << std::endl;
837 std::cout << p <<
"# available processes = " <<
Comm().NumProc() << std::endl;
838 std::cout << p <<
"# processes used in computation = " <<
nprow_ *
npcol_ 839 <<
" ( = " <<
nprow_ <<
"x" <<
npcol_ <<
")" << std::endl;
840 std::cout << p <<
"Equil = " <<
Equil_ << std::endl;
841 std::cout << p <<
"ColPerm = " <<
ColPerm_ << std::endl;
842 std::cout << p <<
"RowPerm = " <<
RowPerm_ << std::endl;
843 std::cout << p <<
"IterRefine = " <<
IterRefine_ << std::endl;
845 std::cout << p <<
"AddZeroToDiag = " <<
AddZeroToDiag_ << std::endl;
846 std::cout << p <<
"Redistribute = " <<
Redistribute_ << std::endl;
856 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
872 std::string p =
"Amesos_Superludist : ";
875 std::cout << p <<
"Time to convert matrix to Superludist format = " 876 << ConTime <<
" (s)" << std::endl;
877 std::cout << p <<
"Time to redistribute matrix = " 878 << MatTime <<
" (s)" << std::endl;
879 std::cout << p <<
"Time to redistribute vectors = " 880 << VecTime <<
" (s)" << std::endl;
881 std::cout << p <<
"Number of symbolic factorizations = " 883 std::cout << p <<
"Time for sym fact = 0.0 (s), avg = 0.0 (s)" << std::endl;
884 std::cout << p <<
"Number of numeric factorizations = " 886 std::cout << p <<
"Time for num fact = " 887 << NumTime *
NumNumericFact_ <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
888 std::cout << p <<
"Number of solve phases = " 890 std::cout << p <<
"Time for solve = " 891 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
896 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
897 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
898 std::cout << p <<
"(the above time does not include SuperLU_DIST time)" << std::endl;
899 std::cout << p <<
"Amesos interface time / total time = " << OveTime / tt << std::endl;
void PrintTiming() const
Print various timig.
int NumSymbolicFact_
Number of symbolic factorization phases.
double GetTime(const std::string what) const
Gets the cumulative time using the string.
int Solve()
Solves A X = B (or AT x = B)
bool ReuseSymbolic_
Allows FactOption to be used on subsequent calls to pdgssvx from NumericFactorization.
void PrintStatus() const
Print various information about the parameters used by Superludist.
bool FactorizationOK_
true if NumericFactorization() has been successfully called.
bool UseTranspose() const
Always returns false. Superludist doesn't support transpose solve.
~Amesos_Superludist(void)
Amesos_Superludist Destructor.
bool Redistribute_
redistribute the input matrix prior to calling Superludist
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
ScalePermstruct_t ScalePermstruct_
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void SetMaxProcesses(int &MaxProcesses, const Epetra_RowMatrix &A)
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
int NumNumericFact_
Number of numeric factorization phases.
int GridCreated_
true if the SuperLU_DIST's grid has been created (and has to be free'd)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int NumSolve_
Number of solves.
Epetra_CrsMatrix & CrsUniformMatrix()
Epetra_RowMatrix * RowMatrixA_
std::vector< double > Aval_
const Epetra_Import & Importer() const
int NumericFactorization()
Performs NumericFactorization on the matrix A.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
RCP< Epetra_Map > UniformMap_
#define AMESOS_CHK_ERR(a)
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.
gridinfo_t grid_
SuperLU_DIST's grid information.
int * Global_Columns_
Contains the global ID of local columns.
bool PrintStatus_
If true, print additional information in the destructor.
Teuchos::RCP< Amesos_Superlu_Pimpl > PrivateSuperluData_
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
RCP< Epetra_RowMatrix > UniformMatrix_
int NumGlobalRows_
Global dimension of the matrix.
int SetNPRowAndCol(const int MaxProcesses, int &nprow, int &npcol)
const Epetra_RowMatrix & UniformMatrix() const
Amesos_Superludist(const Epetra_LinearProblem &LinearProblem)
Amesos_Superludist Constructor.
Teuchos::RCP< Epetra_Import > Importer_
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.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
SOLVEstruct_t SOLVEstruct_
int Superludist_NumProcRows(int NumProcs)
const Epetra_LinearProblem * Problem_
bool MatrixShapeOK() const
Returns true if SUPERLUDIST can handle this matrix shape.
RCP< Epetra_CrsMatrix > CrsUniformMatrix_
void PrintLine() const
Prints line on std::cout.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
superlu_options_t options_
Vector of options.
double AddToDiag_
Add this value to the diagonal.