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