30#include "Epetra_Map.h"
31#include "Epetra_Import.h"
32#include "Epetra_CrsMatrix.h"
33#include "Epetra_Vector.h"
34#include "Epetra_Util.h"
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
142 int NumGlobalRows =
Matrix_->NumGlobalRows();
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;
194 ierr =
SerialMatrix().ExtractMyRowCopy(i, MaxNumEntries,
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];
232 if (ParameterList.isSublist(
"Pardiso"))
234 param_ = ParameterList.sublist(
"Pardiso");
238 iparm_[0] =
param_.get<
int>(
"No default parameters", 1);
239 iparm_[1] =
param_.get<
int>(
"Use METIS reordering" , 2);
241 iparm_[3] =
param_.get<
int>(
"Do preconditioned CGS iterations", 0);
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);
248 iparm_[12] =
param_.get<
int>(
"Use (non-)symmetric matchings", 0);
249 iparm_[17] =
param_.get<
int>(
"Number of non-zeros in LU; -1 to compute", -1);
250 iparm_[18] =
param_.get<
int>(
"Mflops for LU fact; -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
309 Comm().Broadcast( &error, 1, 0 );
311 if (
Comm().MyPID() == 0) {
327 if (
Comm().MyPID() == 0)
334#ifdef HAVE_AMESOS_PARDISO_MKL
348 Comm().Broadcast( &error, 1, 0 );
350 if (
Comm().MyPID() == 0) {
364 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
365 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() )
392 if (
Comm().NumProc() != 1)
441 Epetra_MultiVector* X =
Problem_->GetLHS();
442 Epetra_MultiVector* B =
Problem_->GetRHS();
444 if ((X == 0) || (B == 0))
447 int NumVectors = X->NumVectors();
448 if (NumVectors != B->NumVectors())
452 Epetra_MultiVector* SerialB;
453 Epetra_MultiVector* SerialX;
457 if (
Comm().NumProc() == 1)
464 SerialX =
new Epetra_MultiVector(
SerialMap(),NumVectors);
465 SerialB =
new Epetra_MultiVector(
SerialMap(),NumVectors);
467 SerialB->Import(*B,
Importer(),Insert);
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,
512 Comm().Broadcast( &error, 1, 0 );
514 if (
Comm().MyPID() == 0) {
524 if (
Comm().NumProc() != 1)
526 X->Export(*SerialX,
Importer(), Insert);
547 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
550 std::string p =
"Amesos_Pardiso : ";
553 int n =
Matrix().NumGlobalRows();
554 int nnz =
Matrix().NumGlobalNonzeros();
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;
581 if (
Problem_->GetOperator() == 0 ||
Comm().MyPID() != 0)
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;
#define AMESOS_CHK_ERRV(amesos_err)
#define AMESOS_CHK_ERR(a)
#define AMESOS_RETURN(amesos_err)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
double AddToDiag_
Add this value to the diagonal.
Teuchos::RCP< Epetra_RowMatrix > SerialMatrix_
~Amesos_Pardiso()
Destructor.
bool pardiso_initialized_
bool MatrixShapeOK() const
Returns true if PARDISO can handle this matrix shape.
Epetra_Import & Importer()
Amesos_Pardiso(const Epetra_LinearProblem &LinearProblem)
Constructor.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
int SetParameters(Teuchos::ParameterList &ParameterList)
Set parameters from the input parameters list, returns 0 if successful.
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrix_
Teuchos::RCP< Epetra_Map > SerialMap_
const Epetra_RowMatrix * Matrix_
Epetra_RowMatrix & SerialMatrix()
const Epetra_Map & Map() const
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int MtxConvTime_
Quick access pointers to the internal timing data.
bool UseTranspose_
If true, the transpose of A is used.
const Epetra_RowMatrix & Matrix() const
int msglvl_
Actual matrix for solution phase (always 1)
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
bool UseTranspose() const
Returns the current UseTranspose setting.
int CheckError(const int error) const
int Solve()
Solves A X = B (or AT X = B)
Teuchos::RCP< Epetra_Import > Importer_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
std::vector< double > aa_
void PrintStatus() const
Prints information about the factorization and solution phases.
void PrintTiming() const
Prints timing information.
int PerformNumericFactorization()
Teuchos::ParameterList param_
Epetra_CrsMatrix & SerialCrsMatrix()
int PerformSymbolicFactorization()
int debug_
Sets the level of debug_ output.
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.