MFEM  v4.5.1
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  hypreParMat.HostRead();
59  hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
60  hypreParMat.HypreRead();
61 #if MFEM_HYPRE_VERSION >= 21600
62  hypre_CSRMatrixBigJtoJ(csr_op);
63 #endif
64 
65  m = parcsr_op->global_num_rows;
66  first_row = parcsr_op->first_row_index;
67  nnz_loc = csr_op->num_nonzeros;
68  m_loc = csr_op->num_rows;
69 
70  height = m_loc;
71  width = m_loc;
72 
73  double *csr_nzval = csr_op->data;
74  int *csr_colind = csr_op->j;
75 
76  delete[] csr_rowptr;
77  delete[] reordered_csr_colind;
78  delete[] reordered_csr_nzval;
79  csr_rowptr = new int[m_loc + 1];
80  reordered_csr_colind = new int[nnz_loc];
81  reordered_csr_nzval = new double[nnz_loc];
82 
83  for (int i = 0; i <= m_loc; i++)
84  {
85  csr_rowptr[i] = (csr_op->i)[i];
86  }
87 
88  // CPardiso expects the column indices to be sorted for each row
89  std::vector<int> permutation_idx(nnz_loc);
90  std::iota(permutation_idx.begin(), permutation_idx.end(), 0);
91  for (int i = 0; i < m_loc; i++)
92  {
93  std::sort(permutation_idx.begin() + csr_rowptr[i],
94  permutation_idx.begin() + csr_rowptr[i + 1],
95  [csr_colind](int i1, int i2)
96  {
97  return csr_colind[i1] < csr_colind[i2];
98  });
99  }
100 
101  for (int i = 0; i < nnz_loc; i++)
102  {
103  reordered_csr_colind[i] = csr_colind[permutation_idx[i]];
104  reordered_csr_nzval[i] = csr_nzval[permutation_idx[i]];
105  }
106 
107  hypre_CSRMatrixDestroy(csr_op);
108 
109  // The number of row in global matrix, rhs element and solution vector that
110  // begins the input domain belonging to this MPI process
111  iparm[40] = first_row;
112 
113  // The number of row in global matrix, rhs element and solution vector that
114  // ends the input domain belonging to this MPI process
115  iparm[41] = first_row + m_loc - 1;
116 
117  // Analyze inputs
118  phase = 11;
119  cluster_sparse_solver(pt,
120  &maxfct,
121  &mnum,
122  &mtype,
123  &phase,
124  &m,
125  reordered_csr_nzval,
126  csr_rowptr,
127  reordered_csr_colind,
128  &idum,
129  &nrhs,
130  iparm,
131  &msglvl,
132  &ddum,
133  &ddum,
134  &comm_,
135  &error);
136 
137  MFEM_ASSERT(error == 0, "CPardiso analyze input error");
138 
139  // Numerical factorization
140  phase = 22;
141  cluster_sparse_solver(pt,
142  &maxfct,
143  &mnum,
144  &mtype,
145  &phase,
146  &m,
147  reordered_csr_nzval,
148  csr_rowptr,
149  reordered_csr_colind,
150  &idum,
151  &nrhs,
152  iparm,
153  &msglvl,
154  &ddum,
155  &ddum,
156  &comm_,
157  &error);
158 
159  MFEM_ASSERT(error == 0, "CPardiso factorization input error");
160 }
161 
162 void CPardisoSolver::Mult(const Vector &b, Vector &x) const
163 {
164  // Solve
165  phase = 33;
166  cluster_sparse_solver(pt,
167  &maxfct,
168  &mnum,
169  &mtype,
170  &phase,
171  &m,
172  reordered_csr_nzval,
173  csr_rowptr,
174  reordered_csr_colind,
175  &idum,
176  &nrhs,
177  iparm,
178  &msglvl,
179  b.GetData(),
180  x.GetData(),
181  &comm_,
182  &error);
183 
184  MFEM_ASSERT(error == 0, "Pardiso solve error");
185 }
186 
187 void CPardisoSolver::SetPrintLevel(int print_level)
188 {
189  msglvl = print_level;
190 }
191 
193 {
194  mtype = mat_type;
195 }
196 
198 {
199  // Release all internal memory
200  phase = -1;
201  cluster_sparse_solver(pt,
202  &maxfct,
203  &mnum,
204  &mtype,
205  &phase,
206  &m,
207  reordered_csr_nzval,
208  csr_rowptr,
209  reordered_csr_colind,
210  &idum,
211  &nrhs,
212  iparm,
213  &msglvl,
214  &ddum,
215  &ddum,
216  &comm_,
217  &error);
218 
219  MFEM_ASSERT(error == 0, "CPardiso free error");
220 
221  delete[] csr_rowptr;
222  delete[] reordered_csr_colind;
223  delete[] reordered_csr_nzval;
224 }
225 
226 } // namespace mfem
227 
228 #endif // MFEM_USE_MKL_CPARDISO
229 #endif // MFEM_USE_MPI
void Mult(const Vector &b, Vector &x) const override
Solve.
Definition: cpardiso.cpp:162
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:209
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:187
void SetMatrixType(MatType mat_type)
Set the matrix type.
Definition: cpardiso.cpp:192
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28