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;
#define AMESOS_CHK_ERR(a)
int SetNPRowAndCol(const int MaxProcesses, int &nprow, int &npcol)
int Superludist_NumProcRows(int NumProcs)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
bool AddZeroToDiag_
Adds zero to diagonal of redistributed matrix (some solvers choke on a matrix with a partly empty dia...
double AddToDiag_
Add this value to the diagonal.
bool PrintTiming_
If true, prints timing information in the destructor.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int NumSymbolicFact_
Number of symbolic factorization phases.
int NumSolve_
Number of solves.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
bool PrintStatus_
If true, print additional information in the destructor.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
int NumNumericFact_
Number of numeric factorization phases.
superlu_options_t options_
Vector of options.
gridinfo_t grid_
SuperLU_DIST's grid information.
SOLVEstruct_t SOLVEstruct_
ScalePermstruct_t ScalePermstruct_
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
Teuchos::RCP< Amesos_Superlu_Pimpl > PrivateSuperluData_
void PrintTiming() const
Print various timig.
Epetra_RowMatrix * RowMatrixA_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
int Solve()
Solves A X = B (or AT x = B)
int NumGlobalRows_
Global dimension of the matrix.
int * Global_Columns_
Contains the global ID of local columns.
Epetra_CrsMatrix & CrsUniformMatrix()
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
void PrintStatus() const
Print various information about the parameters used by Superludist.
bool UseTranspose() const
Always returns false. Superludist doesn't support transpose solve.
const Epetra_Import & Importer() const
RCP< Epetra_CrsMatrix > CrsUniformMatrix_
RCP< Epetra_Map > UniformMap_
Amesos_Superludist(const Epetra_LinearProblem &LinearProblem)
Amesos_Superludist Constructor.
Teuchos::RCP< Epetra_Import > Importer_
std::vector< double > Aval_
~Amesos_Superludist(void)
Amesos_Superludist Destructor.
bool Redistribute_
redistribute the input matrix prior to calling Superludist
int GridCreated_
true if the SuperLU_DIST's grid has been created (and has to be free'd)
bool FactorizationOK_
true if NumericFactorization() has been successfully called.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
const Epetra_LinearProblem * Problem_
bool ReuseSymbolic_
Allows FactOption to be used on subsequent calls to pdgssvx from NumericFactorization.
const Epetra_RowMatrix & UniformMatrix() const
bool MatrixShapeOK() const
Returns true if SUPERLUDIST can handle this matrix shape.
RCP< Epetra_RowMatrix > UniformMatrix_
double GetTime(const std::string what) const
Gets the cumulative time using the string.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
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 CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
void SetMaxProcesses(int &MaxProcesses, const Epetra_RowMatrix &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.
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 PrintLine() const
Prints line on std::cout.