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_IntSerialDenseVector.h"
38#include "Epetra_SerialDenseVector.h"
39#include "Epetra_IntVector.h"
40#include "Epetra_Vector.h"
41#include "Epetra_SerialDenseMatrix.h"
42#include "Epetra_Util.h"
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)
129 Teuchos::ParameterList ParamList;
154 MPI_Comm_free( &MUMPSComm_ );
174 Epetra_RowMatrix* ptr;
175 if (
Comm().NumProc() == 1)
183#ifdef EXTRA_DEBUG_INFO
184 Epetra_CrsMatrix* Eptr =
dynamic_cast<Epetra_CrsMatrix*
>( ptr );
185 if ( ptr->NumGlobalNonzeros() < 300 )
SetICNTL(4,3 );
186 if ( ptr->NumGlobalNonzeros() < 42 && Eptr ) {
187 std::cout <<
" Matrix = " << std::endl ;
188 Eptr->Print( std::cout ) ;
194 Row.resize(ptr->NumMyNonzeros());
195 Col.resize(ptr->NumMyNonzeros());
196 Val.resize(ptr->NumMyNonzeros());
198 int MaxNumEntries = ptr->MaxNumEntries();
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) {
208 int GlobalRow = ptr->RowMatrixRowMap().GID(i);
212 ierr = ptr->ExtractMyRowCopy(i, MaxNumEntries,
213 NumEntries, &Values[0],
217 for (
int j = 0 ; j < NumEntries ; ++j) {
218 if (OnlyValues ==
false) {
219 Row[count] = GlobalRow + 1;
220 Col[count] = ptr->RowMatrixColMap().GID(Indices[j]) + 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;
312 if (ParameterList.isParameter(
"NoDestroy"))
313 NoDestroy_ = ParameterList.get<
bool>(
"NoDestroy");
321 if (ParameterList.isSublist(
"mumps"))
323 Teuchos::ParameterList MumpsParams = ParameterList.sublist(
"mumps") ;
326 for (
int i = 1 ; i <= 40 ; ++i)
329 sprintf(what,
"ICNTL(%d)", i);
330 if (MumpsParams.isParameter(what))
331 SetICNTL(i, MumpsParams.get<
int>(what));
335 for (
int i = 1 ; i <= 5 ; ++i)
338 sprintf(what,
"CNTL(%d)", i);
339 if (MumpsParams.isParameter(what))
340 SetCNTL(i, MumpsParams.get<
double>(what));
344 if (MumpsParams.isParameter(
"PermIn"))
347 PermIn = MumpsParams.get<
int*>(
"PermIn");
351 if (MumpsParams.isParameter(
"ColScaling"))
354 ColSca = MumpsParams.get<
double *>(
"ColScaling");
358 if (MumpsParams.isParameter(
"RowScaling"))
361 RowSca = MumpsParams.get<
double *>(
"RowScaling");
374#ifndef HAVE_AMESOS_MPI_C2F
380 int NumGlobalNonzeros, NumRows;
382 NumGlobalNonzeros =
Matrix().NumGlobalNonzeros();
383 NumRows =
Matrix().NumGlobalRows();
387 int OptNumProcs1 = 1 + EPETRA_MAX(NumRows/10000, NumGlobalNonzeros/100000);
388 OptNumProcs1 = EPETRA_MIN(
Comm().NumProc(),OptNumProcs1);
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_);
460 const Epetra_MpiComm* MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&
Comm());
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());
478 const Epetra_MpiComm* MpiComm =
dynamic_cast<const Epetra_MpiComm*
>(&
Comm());
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)
558 Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ;
559 bool Wrong = AnyWrong > 0 ;
584 if (
Comm().NumProc() != 1)
605 Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ;
606 bool Wrong = AnyWrong > 0 ;
625 Epetra_MultiVector* vecX =
Problem_->GetLHS();
626 Epetra_MultiVector* vecB =
Problem_->GetRHS();
627 int NumVectors = vecX->NumVectors();
629 if ((vecX == 0) || (vecB == 0))
632 if (NumVectors != vecB->NumVectors())
635 if (
Comm().NumProc() == 1)
638 for (
int j = 0 ; j < NumVectors; j++)
644 for (
int i = 0 ; i <
Matrix().NumMyRows() ; ++i)
645 (*vecX)[j][i] = (*vecB)[j][i];
646 MDS.rhs = (*vecX)[j];
655 Epetra_MultiVector SerialVector(
SerialMap(),NumVectors);
661 for (
int j = 0 ; j < NumVectors; j++)
663 if (
Comm().MyPID() == 0)
664 MDS.rhs = SerialVector[j];
695Epetra_CrsMatrix * Amesos_Mumps::GetCrsSchurComplement()
706 if(
Comm().MyPID() == 0 )
724Epetra_SerialDenseMatrix * Amesos_Mumps::GetDenseSchurComplement()
729 if(
Comm().MyPID() != 0 )
return 0;
737 (*DenseSchurComplement_)(i,j) =
MDS.schur[pos];
749int Amesos_Mumps::ComputeSchurComplement(
bool flag,
int NumSchurComplementRows,
750 int * SchurComplementRows)
757 if(
Comm().MyPID() == 0 )
770 if (
Comm().MyPID() != 0 )
return;
775 std::cout <<
"Amesos_Mumps : Matrix has " <<
Matrix().NumGlobalRows() <<
" rows"
776 <<
" and " <<
Matrix().NumGlobalNonzeros() <<
" nonzeros" << std::endl;
777 std::cout <<
"Amesos_Mumps : Nonzero elements per row = "
778 << 1.0*
Matrix().NumGlobalNonzeros()/
Matrix().NumGlobalRows() << std::endl;
779 std::cout <<
"Amesos_Mumps : Percentage of nonzero elements = "
780 << 100.0*
Matrix().NumGlobalNonzeros()/(pow(
Matrix().NumGlobalRows(),2.0)) << std::endl;
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))
821 if (
Comm().MyPID() == 0 && Wrong)
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);
956 int i =
Matrix().NumGlobalRows();
957 if (
Comm().MyPID()) i = 0;
const int NumericallySingularMatrixError
const int StructurallySingularMatrixError
#define AMESOS_CHK_ERR(a)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
double AddToDiag_
Add this value to the diagonal.
int MatrixProperty_
Set the matrix property.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
RCP< Epetra_Map > SerialMap_
Map with all elements on process 0 (for solution and rhs).
double * GetRINFO()
Gets the pointer to the RINFO array (defined on all processes).
bool UseTranspose_
If true, solve the problem with AT.
Epetra_RowMatrix & Matrix()
Returns a reference to the linear system matrix.
void SetCNTL(int pos, double value)
Set CNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
void SetICNTL(int pos, int value)
Set ICNTL[pos] to value. pos is expressed in FORTRAN style (starting from 1).
Epetra_Map & RedistrMap()
Returns a reference to the map for redistributed matrix.
double * RowSca_
Row and column scaling.
int MtxConvTime_
Quick access pointers to the internal timers.
void PrintStatus() const
Prints information about the factorization and solution phases.
std::map< int, double > CNTL
Epetra_Import & RedistrImporter()
Returns a reference for the redistributed importer.
std::map< int, int > ICNTL
int NumericFactorization()
Performs NumericFactorization on the matrix A.
Epetra_RowMatrix & RedistrMatrix(const bool ImportMatrix=false)
Returns a reference to the redistributed matrix, imports it is ImportMatrix is true.
int * GetINFO()
Gets the pointer to the INFO array (defined on all processes).
int * SchurComplementRows_
Rows for the Schur complement (if required)
bool IsComputeSchurComplementOK_
true if the Schur complement has been computed (need to free memory)
RCP< Epetra_CrsMatrix > CrsSchurComplement_
Pointer to the Schur complement, as CrsMatrix.
int SetColScaling(double *ColSca)
Set column scaling.
int NumSchurComplementRows_
Number of rows in the Schur complement (if required)
std::vector< double > Val
values of nonzero elements
int SetOrdering(int *PermIn)
Sets ordering.
Amesos_Mumps(const Epetra_LinearProblem &LinearProblem)
Amesos_Mumps Constructor.
int MaxProcs_
Maximum number of processors in the MUMPS' communicator.
double * GetRINFOG()
Gets the pointer to the RINFOG array (defined on host only).
RCP< Epetra_Map > RedistrMap_
Redistributed matrix.
RCP< Epetra_CrsMatrix > RedistrMatrix_
Redistributed matrix (only if MaxProcs_ > 1).
RCP< Epetra_Import > RedistrImporter_
Redistributed importer (from Matrix().RowMatrixRowMap() to RedistrMatrix().RowMatrixRowMap()).
std::vector< int > Col
column indices of nonzero elements
int * GetINFOG()
Get the pointer to the INFOG array (defined on host only).
std::vector< int > Row
row indices of nonzero elements
void PrintTiming() const
Prints timing information.
bool UseTranspose() const
Returns the current UseTranspose setting.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int * PermIn_
PermIn for MUMPS.
Epetra_Import & SerialImporter()
Returns a reference to the importer for SerialMap().
int CheckError()
Checks for MUMPS error, prints them if any. See MUMPS' manual.
int ConvertToTriplet(const bool OnlyValues)
Converts to MUMPS format (COO format).
const Epetra_LinearProblem * Problem_
Pointer to the linear problem to be solved.
void CheckParameters()
Check input parameters.
RCP< Epetra_Import > SerialImporter_
Importer from Matrix.OperatorDomainMap() to SerialMap_.
RCP< Epetra_SerialDenseMatrix > DenseSchurComplement_
Pointer to the Schur complement,as DenseMatrix.
int SetRowScaling(double *RowSca)
Set row scaling.
Epetra_Map & SerialMap()
Returns a reference to the map with all elements on process 0.
~Amesos_Mumps(void)
Amesos_Mumps Destructor.
void Destroy()
Destroys all data associated with \sl this object.
bool PrintTiming_
If true, prints timing information in the destructor.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int verbose_
Toggles the output level.
int NumSymbolicFact_
Number of symbolic factorization phases.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
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)
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
int NumNumericFact_
Number of numeric factorization phases.
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 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.