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" 56 #if SUPERLU_DIST_MAJOR_VERSION == 5 74 for ( i = 1; i*i <= NumProcs; i++ )
77 for ( NumProcRows = i-1 ; done == false ; ) {
78 int NumCols = NumProcs / NumProcRows ;
79 if ( NumProcRows * NumCols == NumProcs )
92 if (nprow < 1 ) nprow = 1;
93 npcol = MaxProcesses / nprow;
95 if( nprow <=0 || npcol <= 0 || MaxProcesses <=0 ) {
96 std::cerr <<
"Amesos_Superludist: wrong value for MaxProcs (" 97 << MaxProcesses <<
"), or nprow (" << nprow
98 <<
") or npcol (" << npcol <<
")" << std::endl;
110 FactorizationDone_(0),
111 FactorizationOK_(false),
115 PrintNonzeros_(false),
120 IterRefine_(
"NOT SET"),
121 ReplaceTinyPivot_(true),
150 Teuchos::ParameterList ParamList;
183 if (ParameterList.isParameter(
"Redistribute"))
188 if (ParameterList.isSublist(
"Superludist") )
190 const Teuchos::ParameterList& SuperludistParams =
191 ParameterList.sublist(
"Superludist") ;
193 if( SuperludistParams.isParameter(
"ReuseSymbolic") )
195 std::string FactOption =
"NotSet";
197 if( SuperludistParams.isParameter(
"Fact") )
198 FactOption = SuperludistParams.get<std::string>(
"Fact");
200 if( FactOption ==
"SamePattern_SameRowPerm" )
PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm;
202 else if ( FactOption !=
"NotSet" )
205 if (SuperludistParams.isParameter(
"Equil"))
206 Equil_ = SuperludistParams.get<
bool>(
"Equil");
208 if (SuperludistParams.isParameter(
"ColPerm"))
209 ColPerm_ = SuperludistParams.get<std::string>(
"ColPerm");
212 if( SuperludistParams.isParameter(
"perm_c"))
213 perm_c_ = SuperludistParams.get<
int*>(
"perm_c");
215 if (SuperludistParams.isParameter(
"RowPerm"))
216 RowPerm_ = SuperludistParams.get<std::string>(
"RowPerm");
218 if (SuperludistParams.isParameter(
"perm_r"))
219 perm_r_ = SuperludistParams.get<
int*>(
"perm_r");
222 if (SuperludistParams.isParameter(
"IterRefine"))
223 IterRefine_ = SuperludistParams.get<std::string>(
"IterRefine");
225 if (SuperludistParams.isParameter(
"ReplaceTinyPivot"))
228 if (SuperludistParams.isParameter(
"PrintNonzeros"))
250 int iam =
Comm().MyPID();
257 int MyFirstElement = iam * m_per_p + EPETRA_MIN( iam, remainder );
258 int MyFirstNonElement = (iam+1) * m_per_p + EPETRA_MIN( iam+1, remainder );
259 int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
262 NumExpectedElements = 0;
267 const Epetra_Map &OriginalMap =
RowMatrixA_->RowMatrixRowMap();
285 for (
int i = 0 ; i <
UniformMap_->NumGlobalElements(); i++ )
318 if (
Comm().NumProc() == 1)
328 if (
Comm().NumProc() == 1)
347 int MyActualFirstElement =
UniformMatrix().RowMatrixRowMap().MinMyGID() ;
350 Ap_.resize( NumMyElements+1 );
351 Ai_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
352 Aval_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ;
361 std::vector<double> RowValuesV_(MaxNumEntries_);
362 std::vector<int> ColIndicesV_(MaxNumEntries_);
366 const Epetra_CrsMatrix *SuperluCrs =
dynamic_cast<const Epetra_CrsMatrix *
>(&
UniformMatrix());
370 for (MyRow = 0; MyRow < NumMyElements ; MyRow++)
374 ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
379 ierr =
UniformMatrix().ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
380 &RowValuesV_[0], &ColIndicesV_[0]);
381 RowValues = &RowValuesV_[0];
382 ColIndices = &ColIndicesV_[0];
389 for (
int i = 0 ; i < NzThisRow ; ++i) {
390 if (ColIndices[i] == MyRow) {
397 Ap_[MyRow] = Ai_index ;
398 for (
int j = 0; j < NzThisRow; j++ ) {
400 Aval_[Ai_index] = RowValues[j] ;
404 assert( NumMyElements == MyRow );
405 Ap_[ NumMyElements ] = Ai_index ;
412 const Epetra_MpiComm & comm1 =
dynamic_cast<const Epetra_MpiComm &
> (
Comm());
443 nnz_loc, NumMyElements, MyActualFirstElement,
445 SLU_NR_loc, SLU_D, SLU_GE );
450 #ifdef HAVE_SUPERLUDIST_LUSTRUCTINIT_2ARG 487 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE 496 #ifdef HAVE_SUPERLUDIST_ENUM_NAMESPACE 581 const Epetra_CrsMatrix *SuperluCrs =
dynamic_cast<const Epetra_CrsMatrix *
>(&
UniformMatrix());
587 std::vector<int> ColIndicesV_(MaxNumEntries_);
588 std::vector<double> RowValuesV_(MaxNumEntries_);
595 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
596 if ( SuperluCrs != 0 ) {
597 ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues,
601 ierr =
UniformMatrix().ExtractMyRowCopy( MyRow, MaxNumEntries_,
602 NzThisRow, &RowValuesV_[0],
604 RowValues = &RowValuesV_[0];
605 ColIndices = &ColIndicesV_[0];
611 for (
int j = 0; j < NzThisRow; j++ ) {
614 Aval_[Ai_index] = RowValues[j] ;
656 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
657 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints())
714 Epetra_MultiVector* vecX =
Problem_->GetLHS();
715 Epetra_MultiVector* vecB =
Problem_->GetRHS();
717 if (vecX == 0 || vecB == 0)
720 int nrhs = vecX->NumVectors() ;
721 if (vecB->NumVectors() != nrhs)
727 RCP<Epetra_MultiVector> vec_uni;
731 vec_uni = Teuchos::rcp(
new Epetra_MultiVector(*
UniformMap_, nrhs));
733 vec_uni->Import(*vecB,
Importer(), Insert);
738 vecX->Update(1.0, *vecB, 0.0);
739 vec_uni = Teuchos::rcp(vecX,
false);
753 std::vector<double>berr(nrhs);
777 vecX->Export(*vec_uni,
Importer(), Insert);
783 "Amesos_Superludist");
797 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
800 std::string p =
"Amesos_Superludist : ";
806 <<
" and " << NNZ <<
" nonzeros" << std::endl;
807 std::cout << p <<
"Nonzero elements per row = " 809 std::cout << p <<
"Percentage of nonzero elements = " 811 std::cout << p <<
"Use transpose = " <<
UseTranspose() << std::endl;
812 std::cout << p <<
"Redistribute = " <<
Redistribute_ << std::endl;
813 std::cout << p <<
"# available processes = " <<
Comm().NumProc() << std::endl;
814 std::cout << p <<
"# processes used in computation = " <<
nprow_ *
npcol_ 815 <<
" ( = " <<
nprow_ <<
"x" <<
npcol_ <<
")" << std::endl;
816 std::cout << p <<
"Equil = " <<
Equil_ << std::endl;
817 std::cout << p <<
"ColPerm = " <<
ColPerm_ << std::endl;
818 std::cout << p <<
"RowPerm = " <<
RowPerm_ << std::endl;
819 std::cout << p <<
"IterRefine = " <<
IterRefine_ << std::endl;
821 std::cout << p <<
"AddZeroToDiag = " <<
AddZeroToDiag_ << std::endl;
822 std::cout << p <<
"Redistribute = " <<
Redistribute_ << std::endl;
832 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
848 std::string p =
"Amesos_Superludist : ";
851 std::cout << p <<
"Time to convert matrix to Superludist format = " 852 << ConTime <<
" (s)" << std::endl;
853 std::cout << p <<
"Time to redistribute matrix = " 854 << MatTime <<
" (s)" << std::endl;
855 std::cout << p <<
"Time to redistribute vectors = " 856 << VecTime <<
" (s)" << std::endl;
857 std::cout << p <<
"Number of symbolic factorizations = " 859 std::cout << p <<
"Time for sym fact = 0.0 (s), avg = 0.0 (s)" << std::endl;
860 std::cout << p <<
"Number of numeric factorizations = " 862 std::cout << p <<
"Time for num fact = " 863 << NumTime *
NumNumericFact_ <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
864 std::cout << p <<
"Number of solve phases = " 866 std::cout << p <<
"Time for solve = " 867 << SolTime *
NumSolve_ <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
872 std::cout << p <<
"Total time spent in Amesos = " << tt <<
" (s) " << std::endl;
873 std::cout << p <<
"Total time spent in the Amesos interface = " << OveTime <<
" (s)" << std::endl;
874 std::cout << p <<
"(the above time does not include SuperLU_DIST time)" << std::endl;
875 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.