44 #define ICNTL(I) icntl[(I)-1] 45 #define CNTL(I) cntl[(I)-1] 46 #define INFOG(I) infog[(I)-1] 47 #define INFO(I) info[(I)-1] 48 #define RINFOG(I) rinfog[(I)-1] 52 IsComputeSchurComplementOK_(false),
65 NumSchurComplementRows_(-1),
82 #if defined(MUMPS_4_9) || defined(MUMPS_5_0) 154 MPI_Comm_free( &MUMPSComm_ );
175 if (
Comm().NumProc() == 1)
183 #ifdef EXTRA_DEBUG_INFO 187 std::cout <<
" Matrix = " << std::endl ;
188 Eptr->
Print( std::cout ) ;
199 std::vector<int> Indices;
200 std::vector<double> Values;
201 Indices.resize(MaxNumEntries);
202 Values.resize(MaxNumEntries);
206 for (
int i = 0; i < ptr->
NumMyRows() ; ++i) {
213 NumEntries, &Values[0],
217 for (
int j = 0 ; j < NumEntries ; ++j) {
218 if (OnlyValues ==
false) {
219 Row[count] = GlobalRow + 1;
227 Val[count] = Values[j];
234 assert (count <= ptr->NumMyNonzeros());
256 std::map<int,int>::iterator i_iter;
257 for (i_iter =
ICNTL.begin() ; i_iter !=
ICNTL.end() ; ++i_iter)
259 int pos = i_iter->first;
260 int val = i_iter->second;
261 if (pos < 1 || pos > 40)
continue;
262 MDS.ICNTL(pos) = val;
265 std::map<int,double>::iterator d_iter;
266 for (d_iter =
CNTL.begin() ; d_iter !=
CNTL.end() ; ++d_iter)
268 int pos = d_iter->first;
269 double val = d_iter->second;
270 if (pos < 1 || pos > 5)
continue;
276 if (
Comm().NumProc() != 1)
MDS.ICNTL(18)= 3;
277 else MDS.ICNTL(18)= 0;
280 else MDS.ICNTL(9) = 1;
286 else MDS.ICNTL(19) = 0;
326 for (
int i = 1 ; i <= 40 ; ++i)
329 sprintf(what,
"ICNTL(%d)", i);
335 for (
int i = 1 ; i <= 5 ; ++i)
338 sprintf(what,
"CNTL(%d)", i);
347 PermIn = MumpsParams.
get<
int*>(
"PermIn");
354 ColSca = MumpsParams.
get<
double *>(
"ColScaling");
361 RowSca = MumpsParams.
get<
double *>(
"RowScaling");
374 #ifndef HAVE_AMESOS_MPI_C2F 380 int NumGlobalNonzeros, NumRows;
387 int OptNumProcs1 = 1 +
EPETRA_MAX(NumRows/10000, NumGlobalNonzeros/100000);
392 int OptNumProcs2 = (int)sqrt(1.0 *
Comm().
NumProc());
393 if (OptNumProcs2 < 1) OptNumProcs2 = 1;
430 #if defined(HAVE_MPI) && defined(HAVE_AMESOS_MPI_C2F) 434 MPI_Comm_free(&MUMPSComm_);
436 std::vector<int> ProcsInGroup(
MaxProcs_);
440 MPI_Group OrigGroup, MumpsGroup;
441 MPI_Comm_group(MPI_COMM_WORLD, &OrigGroup);
442 MPI_Group_incl(OrigGroup,
MaxProcs_, &ProcsInGroup[0], &MumpsGroup);
443 MPI_Comm_create(MPI_COMM_WORLD, MumpsGroup, &MUMPSComm_);
446 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
449 #ifndef HAVE_AMESOS_OLD_MUMPS 450 MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
452 MDS.comm_fortran = (F_INT) MPI_Comm_c2f( MUMPSComm_);
461 assert (MpiComm != 0);
463 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->
GetMpiComm());
466 #ifndef HAVE_AMESOS_OLD_MUMPS 467 MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f(MpiComm->
GetMpiComm());
469 MDS.comm_fortran = (F_INT) MPI_Comm_c2f(MpiComm->
GetMpiComm());
479 assert (MpiComm != 0);
480 MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->
GetMpiComm());
510 if (
Comm().NumProc() != 1)
522 if (
Comm().MyPID() == 0)
559 bool Wrong = AnyWrong > 0 ;
584 if (
Comm().NumProc() != 1)
606 bool Wrong = AnyWrong > 0 ;
629 if ((vecX == 0) || (vecB == 0))
638 for (
int j = 0 ; j < NumVectors; j++)
645 (*vecX)[j][i] = (*vecB)[j][i];
646 MDS.rhs = (*vecX)[j];
661 for (
int j = 0 ; j < NumVectors; j++)
663 if (
Comm().MyPID() == 0)
664 MDS.rhs = SerialVector[j];
706 if(
Comm().MyPID() == 0 )
729 if(
Comm().MyPID() != 0 )
return 0;
737 (*DenseSchurComplement_)(i,j) =
MDS.schur[pos];
749 int Amesos_Mumps::ComputeSchurComplement(
bool flag,
int NumSchurComplementRows,
750 int * SchurComplementRows)
757 if(
Comm().MyPID() == 0 )
770 if (
Comm().MyPID() != 0 )
return;
777 std::cout <<
"Amesos_Mumps : Nonzero elements per row = " 779 std::cout <<
"Amesos_Mumps : Percentage of nonzero elements = " 781 std::cout <<
"Amesos_Mumps : Use transpose = " <<
UseTranspose_ << std::endl;
783 if (
MatrixProperty_ == 0) std::cout <<
"Amesos_Mumps : Matrix is general unsymmetric" << std::endl;
784 if (
MatrixProperty_ == 2) std::cout <<
"Amesos_Mumps : Matrix is general symmetric" << std::endl;
785 if (
MatrixProperty_ == 1) std::cout <<
"Amesos_Mumps : Matrix is SPD" << std::endl;
786 std::cout <<
"Amesos_Mumps : Available process(es) = " <<
Comm().
NumProc() << std::endl;
787 std::cout <<
"Amesos_Mumps : Using " <<
MaxProcs_ <<
" process(es)" << std::endl;
789 std::cout <<
"Amesos_Mumps : Estimated FLOPS for elimination = " 790 <<
MDS.RINFOG(1) << std::endl;
791 std::cout <<
"Amesos_Mumps : Total FLOPS for assembly = " 792 <<
MDS.RINFOG(2) << std::endl;
793 std::cout <<
"Amesos_Mumps : Total FLOPS for elimination = " 794 <<
MDS.RINFOG(3) << std::endl;
796 std::cout <<
"Amesos_Mumps : Total real space to store the LU factors = " 797 <<
MDS.INFOG(9) << std::endl;
798 std::cout <<
"Amesos_Mumps : Total integer space to store the LU factors = " 799 <<
MDS.INFOG(10) << std::endl;
800 std::cout <<
"Amesos_Mumps : Total number of iterative steps refinement = " 801 <<
MDS.INFOG(15) << std::endl;
802 std::cout <<
"Amesos_Mumps : Estimated size of MUMPS internal data\n" 803 <<
"Amesos_Mumps : for running factorization = " 804 <<
MDS.INFOG(16) <<
" Mbytes" << std::endl;
805 std::cout <<
"Amesos_Mumps : for running factorization = " 806 <<
MDS.INFOG(17) <<
" Mbytes" << std::endl;
807 std::cout <<
"Amesos_Mumps : Allocated during factorization = " 808 <<
MDS.INFOG(19) <<
" Mbytes" << std::endl;
816 bool Wrong = ((
MDS.INFOG(1) != 0) || (
MDS.INFO(1) != 0))
823 std::cerr <<
"Amesos_Mumps : ERROR" << std::endl;
824 std::cerr <<
"Amesos_Mumps : INFOG(1) = " <<
MDS.INFOG(1) << std::endl;
825 std::cerr <<
"Amesos_Mumps : INFOG(2) = " <<
MDS.INFOG(2) << std::endl;
828 if (
MDS.INFO(1) != 0 && Wrong)
830 std::cerr <<
"Amesos_Mumps : On process " <<
Comm().
MyPID()
831 <<
", INFO(1) = " <<
MDS.INFO(1) << std::endl;
832 std::cerr <<
"Amesos_Mumps : On process " <<
Comm().
MyPID()
833 <<
", INFO(2) = " <<
MDS.INFO(2) << std::endl;
859 std::string p =
"Amesos_Mumps : ";
862 std::cout << p <<
"Time to convert matrix to MUMPS format = " 863 << ConTime <<
" (s)" << std::endl;
864 std::cout << p <<
"Time to redistribute matrix = " 865 << MatTime <<
" (s)" << std::endl;
866 std::cout << p <<
"Time to redistribute vectors = " 867 << VecTime <<
" (s)" << std::endl;
868 std::cout << p <<
"Number of symbolic factorizations = " 870 std::cout << p <<
"Time for sym fact = " 871 << SymTime <<
" (s), avg = " << SymTime <<
" (s)" << std::endl;
872 std::cout << p <<
"Number of numeric factorizations = " 874 std::cout << p <<
"Time for num fact = " 875 << NumTime <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
876 std::cout << p <<
"Number of solve phases = " 878 std::cout << p <<
"Time for solve = " 879 << SolTime <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
903 assert (
Comm().NumProc() != 1);
907 if (
Comm().MyPID() == 0)
921 assert (
Comm().NumProc() != 1);
957 if (
Comm().MyPID()) i = 0;
RCP< Epetra_Import > SerialImporter_
Importer from Matrix.OperatorDomainMap() to SerialMap_.
int NumSymbolicFact_
Number of symbolic factorization phases.
~Amesos_Mumps(void)
Amesos_Mumps Destructor.
double GetTime(const std::string what) const
Gets the cumulative time using the string.
int Solve()
Solves A X = B (or AT x = B)
virtual const Epetra_Map & RowMatrixRowMap() const=0
int NumericFactorization()
Performs NumericFactorization on the matrix A.
std::map< int, int > ICNTL
RCP< Epetra_CrsMatrix > CrsSchurComplement_
Pointer to the Schur complement, as CrsMatrix.
RCP< Epetra_Map > SerialMap_
Map with all elements on process 0 (for solution and rhs).
int * GetINFO()
Gets the pointer to the INFO array (defined on all processes).
int NumSchurComplementRows_
Number of rows in the Schur complement (if required)
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
T & get(ParameterList &l, const std::string &name)
std::map< int, double > CNTL
int * PermIn_
PermIn for MUMPS.
double * GetRINFO()
Gets the pointer to the RINFO array (defined on all processes).
int MtxConvTime_
Quick access pointers to the internal timers.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
bool isSublist(const std::string &name) const
void SetCNTL(int pos, double value)
Set CNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
RCP< ParameterList > sublist(const RCP< ParameterList > ¶mList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int * SchurComplementRows_
Rows for the Schur complement (if required)
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
RCP< Epetra_Map > RedistrMap_
Redistributed matrix.
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to be solved.
std::vector< int > Col
column indices of nonzero elements
int NumNumericFact_
Number of numeric factorization phases.
virtual int NumGlobalNonzeros() const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
Epetra_Import & RedistrImporter()
Returns a reference for the redistributed importer.
int NumSolve_
Number of solves.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
Epetra_Import & SerialImporter()
Returns a reference to the importer for SerialMap().
virtual int MaxNumEntries() 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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int MatrixProperty_
Set the matrix property.
int SetColScaling(double *ColSca)
Set column scaling.
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
Epetra_RowMatrix & RedistrMatrix(const bool ImportMatrix=false)
Returns a reference to the redistributed matrix, imports it is ImportMatrix is true.
int ConvertToTriplet(const bool OnlyValues)
Converts to MUMPS format (COO format).
void PrintStatus() const
Prints information about the factorization and solution phases.
int MaxProcs_
Maximum number of processors in the MUMPS' communicator.
MPI_Comm GetMpiComm() const
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.
virtual int NumMyRows() const=0
bool PrintStatus_
If true, print additional information in the destructor.
void SetICNTL(int pos, int value)
Set ICNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
virtual void Print(std::ostream &os) const
RCP< Epetra_CrsMatrix > RedistrMatrix_
Redistributed matrix (only if MaxProcs_ > 1).
bool UseTranspose() const
Returns the current UseTranspose setting.
int SetOrdering(int *PermIn)
Sets ordering.
Amesos_Mumps(const Epetra_LinearProblem &LinearProblem)
Amesos_Mumps Constructor.
Epetra_RowMatrix & Matrix()
Returns a reference to the linear system matrix.
Epetra_MultiVector * GetRHS() const
virtual int NumProc() const=0
std::vector< int > Row
row indices of nonzero elements
int CheckError()
Checks for MUMPS error, prints them if any. See MUMPS' manual.
Epetra_Map & SerialMap()
Returns a reference to the map with all elements on process 0.
void PrintTiming() const
Prints timing information.
const int NumericallySingularMatrixError
RCP< Epetra_SerialDenseMatrix > DenseSchurComplement_
Pointer to the Schur complement,as DenseMatrix.
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
bool isParameter(const std::string &name) const
std::vector< double > Val
values of nonzero elements
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
bool IsComputeSchurComplementOK_
true if the Schur complement has been computed (need to free memory)
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
void CheckParameters()
Check input parameters.
Epetra_MultiVector * GetLHS() const
virtual const Epetra_Map & RowMatrixColMap() const=0
int * GetINFOG()
Get the pointer to the INFOG array (defined on host only).
int SetRowScaling(double *RowSca)
Set row scaling.
int verbose_
Toggles the output level.
RCP< Epetra_Import > RedistrImporter_
Redistributed importer (from Matrix().RowMatrixRowMap() to RedistrMatrix().RowMatrixRowMap()).
double * RowSca_
Row and column scaling.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
Epetra_Operator * GetOperator() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
bool UseTranspose_
If true, solve the problem with AT.
virtual int NumGlobalRows() const=0
virtual int NumMyNonzeros() const=0
void PrintLine() const
Prints line on std::cout.
void Destroy()
Destroys all data associated with this object.
double * GetRINFOG()
Gets the pointer to the RINFOG array (defined on host only).
Epetra_Map & RedistrMap()
Returns a reference to the map for redistributed matrix.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
double AddToDiag_
Add this value to the diagonal.