36 #ifdef HAVE_AMESOS_PARDISO_MKL 38 #include "mkl_pardiso.h" 39 #define F77_PARDISO PARDISO 43 #define F77_PARDISOINIT F77_FUNC(pardisoinit, PARDISOINIT) 44 #define F77_PARDISO F77_FUNC(pardiso, PARDISO) 48 (
void *,
int *,
int *,
int *,
double *,
int *);
51 (
void *,
int *,
int *,
int *,
int *,
int *,
52 double *,
int *,
int *,
int *,
int *,
int *,
53 int *,
double *,
double *,
int *,
double *);
56 #define IPARM(i) iparm_[i-1] 74 pardiso_initialized_(false)
76 for (
int i = 0; i < 64; i++) {
105 for (
int i = 0; i < 64; i++) {
120 #ifdef HAVE_AMESOS_PARDISO_MKL 146 if (
Comm().MyPID() == 0)
147 NumMyRows = NumGlobalRows;
177 if (
Comm().MyPID() == 0)
184 std::vector<int> Indices(MaxNumEntries);
185 std::vector<double> Values(MaxNumEntries);
193 int ierr, NumEntriesThisRow;
196 &Values[0], &Indices[0]);
200 ia_[i + 1] =
ia_[i] + NumEntriesThisRow;
202 for (
int j = 0 ; j < NumEntriesThisRow ; ++j)
207 ja_[count] = Indices[j] + 1;
208 aa_[count] = Values[j];
244 iparm_[7] =
param_.
get<
int>(
"Max num of iterative refinement steps", 0);
245 iparm_[9] =
param_.
get<
int>(
"Perturbation for pivot elements 10^-k", 13);
246 iparm_[10] =
param_.
get<
int>(
"Use (non-)symmetric scaling vectors", 1);
249 iparm_[17] =
param_.
get<
int>(
"Number of non-zeros in LU; -1 to compute", -1);
262 if (
Comm().MyPID() == 0)
275 #ifndef HAVE_AMESOS_PARDISO_MKL 281 char* var = getenv(
"OMP_NUM_THREADS");
283 sscanf( var,
"%d", &num_procs );
284 IPARM(3) = num_procs;
295 #ifdef HAVE_AMESOS_PARDISO_MKL 311 if (
Comm().MyPID() == 0) {
327 if (
Comm().MyPID() == 0)
334 #ifdef HAVE_AMESOS_PARDISO_MKL 350 if (
Comm().MyPID() == 0) {
364 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
365 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() )
444 if ((X == 0) || (
B == 0))
448 if (NumVectors !=
B->NumVectors())
457 if (
Comm().NumProc() == 1)
476 if (
Comm().MyPID() == 0)
478 double* SerialXValues;
479 double* SerialBValues;
491 for (
int i = 0 ; i < NumVectors ; ++i)
492 #ifdef HAVE_AMESOS_PARDISO_MKL
496 SerialBValues + i *
n,
497 SerialXValues + i *
n,
503 SerialBValues + i *
n,
504 SerialXValues + i *
n,
514 if (
Comm().MyPID() == 0) {
524 if (
Comm().NumProc() != 1)
550 std::string p =
"Amesos_Pardiso : ";
556 std::cout << p <<
"Matrix has " <<
n <<
" rows" 557 <<
" and " << nnz <<
" nonzeros" << std::endl;
558 std::cout << p <<
"Nonzero elements per row = " 559 << 1.0 * nnz /
n << std::endl;
560 std::cout << p <<
"Percentage of nonzero elements = " 561 << 100.0 * nnz /(pow(
n,2.0)) << std::endl;
562 std::cout << p <<
"Use transpose = " <<
UseTranspose_ << std::endl;
563 std::cout << p <<
"Number of performed iterative ref. steps = " <<
IPARM(9) << std::endl;
564 std::cout << p <<
"Peak memory symbolic factorization = " <<
IPARM(15) << std::endl;
565 std::cout << p <<
"Permanent memory symbolic factorization = " <<
IPARM(16) << std::endl;
566 std::cout << p <<
"Memory numerical fact. and solution = " <<
IPARM(17) << std::endl;
567 std::cout << p <<
"Number of nonzeros in factors = " <<
IPARM(18) << std::endl;
568 std::cout << p <<
"MFlops of factorization = " <<
IPARM(19) << std::endl;
569 std::cout << p <<
"CG/CGS diagnostic = " <<
IPARM(20) << std::endl;
570 std::cout << p <<
"Inertia: Number of positive eigenvalues = " <<
IPARM(22) << std::endl;
571 std::cout << p <<
"Inertia: Number of negative eigenvalues = " <<
IPARM(23) << std::endl;
600 std::string p =
"Amesos_Pardiso : ";
603 std::cout << p <<
"Time to convert matrix to Pardiso format = " 604 << ConTime <<
" (s)" << std::endl;
605 std::cout << p <<
"Time to redistribute matrix = " 606 << MatTime <<
" (s)" << std::endl;
607 std::cout << p <<
"Time to redistribute vectors = " 608 << VecTime <<
" (s)" << std::endl;
609 std::cout << p <<
"Number of symbolic factorizations = " 611 std::cout << p <<
"Time for sym fact = " 612 << SymTime <<
" (s), avg = " << SymTime <<
" (s)" << std::endl;
613 std::cout << p <<
"Number of numeric factorizations = " 615 std::cout << p <<
"Time for num fact = " 616 << NumTime <<
" (s), avg = " << NumTime <<
" (s)" << std::endl;
617 std::cout << p <<
"Number of solve phases = " 619 std::cout << p <<
"Time for solve = " 620 << SolTime <<
" (s), avg = " << SolTime <<
" (s)" << std::endl;
633 std::cerr <<
"Amesos: PARDISO returned error code " <<
error << std::endl;
634 std::cerr <<
"Amesos: Related message from manual is:" << std::endl;
639 std::cerr <<
"Input inconsistent" << std::endl;
642 std::cerr <<
"Not enough memory" << std::endl;
645 std::cerr <<
"Reordering problems" << std::endl;
648 std::cerr <<
"Zero pivot, numerical fact. or iterative refinement problem. " << std::endl;
651 std::cerr <<
"Unclassified (internal) error" << std::endl;
654 std::cerr <<
"Preordering failed (matrix types 11, 13 only)" << std::endl;
657 std::cerr <<
"Diagonal matrix problem." << std::endl;
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
int NumSymbolicFact_
Number of symbolic factorization phases.
Teuchos::RCP< Epetra_RowMatrix > SerialMatrix_
double GetTime(const std::string what) const
Gets the cumulative time using the string.
Epetra_Import & Importer()
virtual const Epetra_Map & RowMatrixRowMap() const=0
int PerformNumericFactorization()
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
Amesos_Pardiso(const Epetra_LinearProblem &LinearProblem)
Constructor.
T & get(ParameterList &l, const std::string &name)
bool isSublist(const std::string &name) const
Epetra_CrsMatrix & SerialCrsMatrix()
RCP< ParameterList > sublist(const RCP< ParameterList > ¶mList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrix_
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
const Epetra_Map & Map() const
int CheckError(const int error) const
int NumNumericFact_
Number of numeric factorization phases.
virtual int NumGlobalNonzeros() const=0
virtual int MyPID() const=0
int msglvl_
Actual matrix for solution phase (always 1)
void PrintStatus() const
Prints information about the factorization and solution phases.
int NumSolve_
Number of solves.
void PrintTiming() const
Prints timing information.
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().
std::vector< double > aa_
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
#define AMESOS_CHK_ERR(a)
int SetParameters(Teuchos::ParameterList &ParameterList)
Set parameters from the input parameters list, returns 0 if successful.
const Epetra_RowMatrix * Matrix_
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.
int MtxConvTime_
Quick access pointers to the internal timing data.
virtual int NumMyRows() const=0
bool PrintStatus_
If true, print additional information in the destructor.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
bool UseTranspose() const
Returns the current UseTranspose setting.
#define AMESOS_RETURN(amesos_err)
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
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 int Broadcast(double *MyVals, int Count, int Root) const=0
int PerformSymbolicFactorization()
~Amesos_Pardiso()
Destructor.
Epetra_MultiVector * GetRHS() const
virtual int NumProc() const=0
Epetra_RowMatrix & SerialMatrix()
bool pardiso_initialized_
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 Epetra_RowMatrix & Matrix() const
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Teuchos::RCP< Epetra_Map > SerialMap_
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
Epetra_MultiVector * GetLHS() const
int ExtractView(double **A, int *MyLDA) const
int verbose_
Toggles the output level.
int debug_
Sets the level of debug_ output.
Teuchos::RCP< Epetra_Import > Importer_
int Solve()
Solves A X = B (or AT X = B)
#define AMESOS_CHK_ERRV(amesos_err)
bool UseTranspose_
If true, the transpose of A is used.
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)
virtual int NumGlobalRows() const=0
bool MatrixShapeOK() const
Returns true if PARDISO can handle this matrix shape.
void PrintLine() const
Prints line on std::cout.
Teuchos::ParameterList param_
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
double AddToDiag_
Add this value to the diagonal.