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