Line data Source code
1 : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/diagonalize/DiagonalizerFeast.hpp"
5 :
6 : #include "pairinteraction/enums/FloatType.hpp"
7 : #include "pairinteraction/utils/eigen_assertion.hpp"
8 : #include "pairinteraction/utils/eigen_compat.hpp"
9 :
10 : #include <Eigen/Dense>
11 : #include <cmath>
12 : #include <optional>
13 : #include <spdlog/spdlog.h>
14 :
15 : #ifdef WITH_MKL
16 : #include <fmt/core.h>
17 : #include <fmt/format.h>
18 : #include <mkl.h>
19 : #endif // WITH_MKL
20 :
21 : namespace pairinteraction {
22 : #ifdef WITH_MKL
23 0 : void feast(const char *uplo, const MKL_INT *n, const float *a, const MKL_INT *lda, MKL_INT *fpm,
24 : float *epsout, MKL_INT *loop, const float *emin, const float *emax, MKL_INT *m0,
25 : float *e, float *x, MKL_INT *m, float *res, MKL_INT *info) {
26 0 : sfeast_syev(uplo, n, a, lda, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
27 0 : }
28 :
29 1 : void feast(const char *uplo, const MKL_INT *n, const double *a, const MKL_INT *lda, MKL_INT *fpm,
30 : double *epsout, MKL_INT *loop, const double *emin, const double *emax, MKL_INT *m0,
31 : double *e, double *x, MKL_INT *m, double *res, MKL_INT *info) {
32 1 : dfeast_syev(uplo, n, a, lda, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
33 1 : }
34 :
35 2 : void feast(const char *uplo, const MKL_INT *n, const MKL_Complex8 *a, const MKL_INT *lda,
36 : MKL_INT *fpm, float *epsout, MKL_INT *loop, const float *emin, const float *emax,
37 : MKL_INT *m0, float *e, MKL_Complex8 *x, MKL_INT *m, float *res, MKL_INT *info) {
38 2 : cfeast_heev(uplo, n, a, lda, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
39 2 : }
40 :
41 3 : void feast(const char *uplo, const MKL_INT *n, const MKL_Complex16 *a, const MKL_INT *lda,
42 : MKL_INT *fpm, double *epsout, MKL_INT *loop, const double *emin, const double *emax,
43 : MKL_INT *m0, double *e, MKL_Complex16 *x, MKL_INT *m, double *res, MKL_INT *info) {
44 3 : zfeast_heev(uplo, n, a, lda, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
45 3 : }
46 :
47 : template <typename Scalar>
48 : template <typename ScalarLim>
49 : EigenSystemH<Scalar>
50 6 : DiagonalizerFeast<Scalar>::dispatch_eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
51 : real_t min_eigenvalue, real_t max_eigenvalue,
52 : double rtol) const {
53 : using real_lim_t = typename traits::NumTraits<ScalarLim>::real_t;
54 6 : int dim = matrix.rows();
55 :
56 : // Subtract the mean of the diagonal elements from the diagonal
57 6 : real_t shift{};
58 6 : Eigen::MatrixX<ScalarLim> hamiltonian =
59 : this->template subtract_mean<ScalarLim>(matrix, shift, rtol);
60 :
61 : // Diagonalize the shifted matrix
62 6 : int m0 = std::min(dim, this->m0);
63 :
64 6 : Eigen::VectorX<real_lim_t> evals(dim);
65 6 : Eigen::MatrixX<ScalarLim> evecs(dim, m0); // the first m columns will contain the eigenvectors
66 :
67 6 : double targeted_trace_relative_error = rtol;
68 6 : int precision_feast = static_cast<int>(std::ceil(-std::log10(targeted_trace_relative_error)));
69 :
70 6 : std::vector<MKL_INT> fpm(128);
71 6 : feastinit(fpm.data());
72 6 : fpm[0] = 0; // disable terminal output
73 6 : fpm[1] = 8; // number of contour points
74 6 : fpm[4] = 0; // do not use initial subspace
75 6 : fpm[26] = 0; // disables matrix checker
76 6 : fpm[2] = precision_feast; // single precision stopping criteria
77 6 : fpm[6] = precision_feast; // double precision stopping criteria
78 6 : MKL_INT m{}; // will contain the number of eigenvalues
79 6 : std::vector<real_lim_t> e(m0); // will contain the first m eigenvalues
80 6 : char uplo = 'F'; // full matrix is stored
81 6 : MKL_INT info{}; // will contain return codes
82 6 : real_lim_t epsout{}; // will contain relative error
83 6 : MKL_INT loop{}; // will contain number of used refinement
84 6 : std::vector<real_lim_t> res(m0); // will contain the residual errors
85 6 : real_lim_t min_eigenvalue_lim = min_eigenvalue - shift;
86 6 : real_lim_t max_eigenvalue_lim = max_eigenvalue - shift;
87 :
88 6 : feast(&uplo, &dim, hamiltonian.data(), &dim, fpm.data(), &epsout, &loop, &min_eigenvalue_lim,
89 : &max_eigenvalue_lim, &m0, evals.data(), evecs.data(), &m, res.data(), &info);
90 :
91 : // https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-1/extended-eigensolver-output-details.html
92 6 : if (info != 0) {
93 0 : if (info == 202) {
94 0 : throw std::invalid_argument(
95 : "Diagonalization error: Problem with size of the system n (n≤0).");
96 : }
97 0 : if (info == 201) {
98 0 : throw std::invalid_argument(
99 : "Diagonalization error: Problem with size of initial subspace m0 (m0≤0 or m0>n).");
100 : }
101 0 : if (info == 200) {
102 0 : throw std::invalid_argument(
103 : "Diagonalization error: Problem with emin,emax (emin≥emax).");
104 : }
105 0 : if (info == 3) {
106 0 : throw std::invalid_argument(
107 : "Diagonalization error: Size of the subspace m0 is too small (m0<m).");
108 : }
109 0 : if (info == 2) {
110 0 : throw std::runtime_error(
111 : "Diagonalization error: No convergence (number of iteration loops >fpm[3]).");
112 : }
113 0 : if (info == 1) {
114 0 : throw std::runtime_error(
115 : "Diagonalization error: No eigenvalue found in the search interval.");
116 : }
117 0 : if (info == -1) {
118 0 : throw std::runtime_error(
119 : "Diagonalization error: Internal error for allocation memory.");
120 : }
121 0 : if (info == -2) {
122 0 : throw std::runtime_error(
123 : "Diagonalization error: Internal error of the inner system solver. Possible "
124 : "reasons: not enough memory for inner linear system solver or inconsistent input.");
125 : }
126 0 : if (info == -3) {
127 0 : throw std::runtime_error(
128 : "Diagonalization error: Internal error of the reduced eigenvalue solver. Possible "
129 : "cause: matrix may not be positive definite.");
130 : }
131 0 : if (info == -4) {
132 0 : throw std::invalid_argument("Diagonalization error: Matrix is not positive definite.");
133 : }
134 0 : if (info <= 100) {
135 0 : throw std::invalid_argument(
136 : fmt::format("Diagonalization error: Argument {} to the FEAST *interface* "
137 : "had an illegal value (the counting of the arguments starts with one). "
138 : "For a documentation of the feast interface, see "
139 : "https://www.intel.com/content/www/us/en/docs/onemkl/"
140 : "developer-reference-c/2025-1/feast-syev-feast-heev.html.",
141 0 : -info - 100));
142 : }
143 0 : if (info >= 100) {
144 0 : throw std::invalid_argument(fmt::format(
145 : "Diagonalization error: Argument {} to the FEAST "
146 : "*initialization* had an illegal value (the counting of the arguments starts with "
147 : "one). For a documentation of the feast initialization, see "
148 : "https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-1/"
149 : "extended-eigensolver-input-parameters.html.",
150 0 : info - 100));
151 : }
152 0 : throw std::runtime_error(fmt::format(
153 : "Diagonalization error: The FEAST routine failed with error code {}.", info));
154 : }
155 :
156 : // Restrict to the first m eigenvectors and eigenvalues because the rest is not calculated
157 6 : evecs.conservativeResize(dim, m);
158 6 : evals.conservativeResize(m);
159 :
160 6 : return {evecs.sparseView(1, 0.5 * rtol / std::sqrt(dim)).template cast<Scalar>(),
161 6 : this->add_mean(evals, shift)};
162 12 : }
163 :
164 : template <typename Scalar>
165 3 : DiagonalizerFeast<Scalar>::DiagonalizerFeast(int m0, FloatType float_type)
166 3 : : DiagonalizerInterface<Scalar>(float_type), m0(m0) {
167 3 : if (m0 <= 0) {
168 0 : throw std::invalid_argument("The size of the initial subspace m0 (i.e., the number of "
169 : "maximally obtainable eigenvalues) must be positive.");
170 : }
171 3 : }
172 :
173 : template <typename Scalar>
174 : EigenSystemH<Scalar>
175 0 : DiagonalizerFeast<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/,
176 : double /*rtol*/) const {
177 0 : throw std::invalid_argument("The FEAST routine requires a search interval.");
178 : }
179 :
180 : template <typename Scalar>
181 : EigenSystemH<Scalar>
182 6 : DiagonalizerFeast<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
183 : std::optional<real_t> min_eigenvalue,
184 : std::optional<real_t> max_eigenvalue, double rtol) const {
185 6 : if (!min_eigenvalue.has_value() || !max_eigenvalue.has_value()) {
186 0 : throw std::invalid_argument("The FEAST routine requires a search interval.");
187 : }
188 6 : switch (this->float_type) {
189 2 : case FloatType::FLOAT32:
190 : return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT32>>(
191 2 : matrix, min_eigenvalue.value(), max_eigenvalue.value(), rtol);
192 4 : case FloatType::FLOAT64:
193 : return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT64>>(
194 4 : matrix, min_eigenvalue.value(), max_eigenvalue.value(), rtol);
195 0 : default:
196 0 : throw std::invalid_argument("Unsupported floating point precision.");
197 : }
198 : }
199 :
200 : #else
201 :
202 : template <typename Scalar>
203 : DiagonalizerFeast<Scalar>::DiagonalizerFeast(int m0, FloatType float_type)
204 : : DiagonalizerInterface<Scalar>(float_type), m0(m0) {
205 : throw std::runtime_error(
206 : "The FEAST routine is not available in this build. Please use a different diagonalizer.");
207 : }
208 :
209 : template <typename Scalar>
210 : EigenSystemH<Scalar>
211 : DiagonalizerFeast<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/,
212 : double /*rtol*/) const {
213 : std::abort(); // can't happen because the constructor throws
214 : }
215 :
216 : template <typename Scalar>
217 : EigenSystemH<Scalar>
218 : DiagonalizerFeast<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/,
219 : std::optional<real_t> /*min_eigenvalue*/,
220 : std::optional<real_t> /*max_eigenvalue*/, double /*rtol*/) const {
221 : std::abort(); // can't happen because the constructor throws
222 : }
223 :
224 : #endif // WITH_MKL
225 :
226 : // Explicit instantiations
227 : template class DiagonalizerFeast<double>;
228 : template class DiagonalizerFeast<std::complex<double>>;
229 : } // namespace pairinteraction
|
