MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
pardiso.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
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 "pardiso.hpp"
13
14#ifdef MFEM_USE_MKL_PARDISO
15
16#include "sparsemat.hpp"
17
18namespace mfem
19{
20
22{
23 // Indicate that default parameters are changed
24 iparm[0] = 1;
25 // Use METIS for fill-in reordering
26 iparm[1] = 2;
27 // Write the solution into the x vector data
28 iparm[5] = 0;
29 // Maximum number of iterative refinement steps
30 iparm[7] = 2;
31 // Perturb the pivot elements with 1E-13
32 iparm[9] = 13;
33 // Use nonsymmetric permutation
34 iparm[10] = 1;
35 // Perform a check on the input data
36 iparm[26] = 1;
37#ifdef MFEM_USE_SINGLE
38 // Single precision
39 iparm[27] = 1;
40#endif
41 // 0-based indexing in CSR data structure
42 iparm[34] = 1;
43 // Maximum number of numerical factorizations
44 maxfct = 1;
45 // Which factorization to use. This parameter is ignored and always assumed
46 // to be equal to 1. See MKL documentation.
47 mnum = 1;
48 // Print statistical information in file
49 msglvl = 0;
50 // Initialize error flag
51 error = 0;
52 // Real nonsymmetric matrix
54 // Number of right hand sides
55 nrhs = 1;
56}
57
59{
60 auto mat = const_cast<SparseMatrix *>(dynamic_cast<const SparseMatrix *>(&op));
61
62 MFEM_ASSERT(mat, "Must pass SparseMatrix as Operator");
63
64 height = op.Height();
65
66 width = op.Width();
67
68 m = mat->Size();
69
70 nnz = mat->NumNonZeroElems();
71
72 const int *Ap = mat->HostReadI();
73 const int *Ai = mat->HostReadJ();
74 const real_t *Ax = mat->HostReadData();
75
76 csr_rowptr = new int[m + 1];
77 reordered_csr_colind = new int[nnz];
78 reordered_csr_nzval = new real_t[nnz];
79
80 for (int i = 0; i <= m; i++)
81 {
82 csr_rowptr[i] = Ap[i];
83 }
84
85 // Pardiso expects the column indices to be sorted for each row
86 mat->SortColumnIndices();
87
88 for (int i = 0; i < nnz; i++)
89 {
90 reordered_csr_colind[i] = Ai[i];
91 reordered_csr_nzval[i] = Ax[i];
92 }
93
94 // Analyze inputs
95 phase = 11;
96 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &m, reordered_csr_nzval, csr_rowptr,
97 reordered_csr_colind, &idum, &nrhs,
98 iparm, &msglvl, &ddum, &ddum, &error);
99
100 MFEM_ASSERT(error == 0, "Pardiso symbolic factorization error");
101
102 // Numerical factorization
103 phase = 22;
104 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &m, reordered_csr_nzval, csr_rowptr,
105 reordered_csr_colind, &idum, &nrhs,
106 iparm, &msglvl, &ddum, &ddum, &error);
107
108 MFEM_ASSERT(error == 0, "Pardiso numerical factorization error");
109}
110
111void PardisoSolver::Mult(const Vector &b, Vector &x) const
112{
113 // Solve
114 phase = 33;
115 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &m, reordered_csr_nzval, csr_rowptr,
116 reordered_csr_colind, &idum, &nrhs,
117 iparm, &msglvl, b.GetData(), x.GetData(), &error);
118
119 MFEM_ASSERT(error == 0, "Pardiso solve error");
120}
121
122void PardisoSolver::SetPrintLevel(int print_level)
123{
124 msglvl = print_level;
125}
126
128{
129 mtype = mat_type;
130}
131
133{
134 // Release all internal memory
135 phase = -1;
136 PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &m, reordered_csr_nzval, csr_rowptr,
137 reordered_csr_colind, &idum, &nrhs,
138 iparm, &msglvl, &ddum, &ddum, &error);
139
140 MFEM_ASSERT(error == 0, "Pardiso free error");
141
142 delete[] csr_rowptr;
143 delete[] reordered_csr_colind;
144 delete[] reordered_csr_nzval;
145}
146
147} // namespace mfem
148
149#endif // MFEM_USE_MKL_PARDISO
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
PardisoSolver()
Construct a new PardisoSolver object.
Definition pardiso.cpp:21
void SetMatrixType(MatType mat_type)
Set the matrix type.
Definition pardiso.cpp:127
void SetOperator(const Operator &op) override
Set the Operator object and perform factorization.
Definition pardiso.cpp:58
void SetPrintLevel(int print_lvl)
Set the print level for MKL Pardiso.
Definition pardiso.cpp:122
void Mult(const Vector &b, Vector &x) const override
Solve.
Definition pardiso.cpp:111
Data type sparse matrix.
Definition sparsemat.hpp:51
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