MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
cpardiso.cpp
Go to the documentation of this file.
1 #include "cpardiso.hpp"
2 #include "hypre.hpp"
3 #include <algorithm>
4 #include <vector>
5 #include <numeric>
6 
7 #ifdef MFEM_USE_MPI
8 #ifdef MFEM_USE_MKL_CPARDISO
9 
10 namespace mfem
11 {
13 {
14  comm_ = MPI_Comm_c2f(comm);
15 
16  // Indicate that default parameters are changed
17  iparm[0] = 1;
18  // Use METIS for fill-in reordering
19  iparm[1] = 2;
20  // Do not write the solution into the x vector data
21  iparm[5] = 0;
22  // Maximum number of iterative refinement steps
23  iparm[7] = 2;
24  // Perturb the pivot elements with 1E-13
25  iparm[9] = 13;
26  // Use nonsymmetric permutation
27  iparm[10] = 1;
28  // Perform a check on the input data
29  iparm[26] = 1;
30  // 0-based indexing in CSR data structure
31  iparm[34] = 1;
32  // All inputs are distributed between MPI processes
33  iparm[39] = 2;
34  // Maximum number of numerical factorizations
35  maxfct = 1;
36  // Which factorization to use. This parameter is ignored and always assumed
37  // to be equal to 1. See MKL documentation.
38  mnum = 1;
39  // Print statistical information in file
40  msglvl = 0;
41  // Initialize error flag
42  error = 0;
43  // Real nonsymmetric matrix
44  mtype = MatType::REAL_NONSYMMETRIC;
45  // Number of right hand sides
46  nrhs = 1;
47 };
48 
50 {
51  auto hypreParMat = dynamic_cast<const HypreParMatrix &>(op);
52 
53  MFEM_ASSERT(hypreParMat, "Must pass HypreParMatrix as Operator");
54 
55  auto parcsr_op = static_cast<hypre_ParCSRMatrix *>(
56  const_cast<HypreParMatrix &>(hypreParMat));
57 
58  hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
59 #if MFEM_HYPRE_VERSION >= 21600
60  hypre_CSRMatrixBigJtoJ(csr_op);
61 #endif
62 
63  m = parcsr_op->global_num_rows;
64  first_row = parcsr_op->first_row_index;
65  nnz_loc = csr_op->num_nonzeros;
66  m_loc = csr_op->num_rows;
67 
68  height = m_loc;
69  width = m_loc;
70 
71  double *csr_nzval = csr_op->data;
72  int *csr_colind = csr_op->j;
73 
74  delete[] csr_rowptr;
75  delete[] reordered_csr_colind;
76  delete[] reordered_csr_nzval;
77  csr_rowptr = new int[m_loc + 1];
78  reordered_csr_colind = new int[nnz_loc];
79  reordered_csr_nzval = new double[nnz_loc];
80 
81  for (int i = 0; i <= m_loc; i++)
82  {
83  csr_rowptr[i] = (csr_op->i)[i];
84  }
85 
86  // CPardiso expects the column indices to be sorted for each row
87  std::vector<int> permutation_idx(nnz_loc);
88  std::iota(permutation_idx.begin(), permutation_idx.end(), 0);
89  for (int i = 0; i < m_loc; i++)
90  {
91  std::sort(permutation_idx.begin() + csr_rowptr[i],
92  permutation_idx.begin() + csr_rowptr[i + 1],
93  [csr_colind](int i1, int i2)
94  {
95  return csr_colind[i1] < csr_colind[i2];
96  });
97  }
98 
99  for (int i = 0; i < nnz_loc; i++)
100  {
101  reordered_csr_colind[i] = csr_colind[permutation_idx[i]];
102  reordered_csr_nzval[i] = csr_nzval[permutation_idx[i]];
103  }
104 
105  hypre_CSRMatrixDestroy(csr_op);
106 
107  // The number of row in global matrix, rhs element and solution vector that
108  // begins the input domain belonging to this MPI process
109  iparm[40] = first_row;
110 
111  // The number of row in global matrix, rhs element and solution vector that
112  // ends the input domain belonging to this MPI process
113  iparm[41] = first_row + m_loc - 1;
114 
115  // Analyze inputs
116  phase = 11;
117  cluster_sparse_solver(pt,
118  &maxfct,
119  &mnum,
120  &mtype,
121  &phase,
122  &m,
123  reordered_csr_nzval,
124  csr_rowptr,
125  reordered_csr_colind,
126  &idum,
127  &nrhs,
128  iparm,
129  &msglvl,
130  &ddum,
131  &ddum,
132  &comm_,
133  &error);
134 
135  MFEM_ASSERT(error == 0, "CPardiso analyze input error");
136 
137  // Numerical factorization
138  phase = 22;
139  cluster_sparse_solver(pt,
140  &maxfct,
141  &mnum,
142  &mtype,
143  &phase,
144  &m,
145  reordered_csr_nzval,
146  csr_rowptr,
147  reordered_csr_colind,
148  &idum,
149  &nrhs,
150  iparm,
151  &msglvl,
152  &ddum,
153  &ddum,
154  &comm_,
155  &error);
156 
157  MFEM_ASSERT(error == 0, "CPardiso factorization input error");
158 }
159 
160 void CPardisoSolver::Mult(const Vector &b, Vector &x) const
161 {
162  // Solve
163  phase = 33;
164  cluster_sparse_solver(pt,
165  &maxfct,
166  &mnum,
167  &mtype,
168  &phase,
169  &m,
170  reordered_csr_nzval,
171  csr_rowptr,
172  reordered_csr_colind,
173  &idum,
174  &nrhs,
175  iparm,
176  &msglvl,
177  b.GetData(),
178  x.GetData(),
179  &comm_,
180  &error);
181 
182  MFEM_ASSERT(error == 0, "Pardiso solve error");
183 }
184 
185 void CPardisoSolver::SetPrintLevel(int print_level)
186 {
187  msglvl = print_level;
188 }
189 
191 {
192  mtype = mat_type;
193 }
194 
196 {
197  // Release all internal memory
198  phase = -1;
199  cluster_sparse_solver(pt,
200  &maxfct,
201  &mnum,
202  &mtype,
203  &phase,
204  &m,
205  reordered_csr_nzval,
206  csr_rowptr,
207  reordered_csr_colind,
208  &idum,
209  &nrhs,
210  iparm,
211  &msglvl,
212  &ddum,
213  &ddum,
214  &comm_,
215  &error);
216 
217  MFEM_ASSERT(error == 0, "CPardiso free error");
218 
219  delete[] csr_rowptr;
220  delete[] reordered_csr_colind;
221  delete[] reordered_csr_nzval;
222 }
223 
224 } // namespace mfem
225 
226 #endif // MFEM_USE_MKL_CPARDISO
227 #endif // MFEM_USE_MPI
void Mult(const Vector &b, Vector &x) const override
Solve.
Definition: cpardiso.cpp:160
CPardisoSolver(MPI_Comm comm)
Construct a new CPardisoSolver object.
Definition: cpardiso.cpp:12
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
double b
Definition: lissajous.cpp:42
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void SetOperator(const Operator &op) override
Set the Operator object and perform factorization.
Definition: cpardiso.cpp:49
Vector data type.
Definition: vector.hpp:60
void SetPrintLevel(int print_lvl)
Set the print level for MKL CPardiso.
Definition: cpardiso.cpp:185
void SetMatrixType(MatType mat_type)
Set the matrix type.
Definition: cpardiso.cpp:190
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28