Line data Source code
1 : // SPDX-FileCopyrightText: 2025 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvr.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 <fmt/core.h>
13 : #include <limits>
14 : #include <spdlog/spdlog.h>
15 :
16 : #ifdef WITH_MKL
17 : #include <mkl.h>
18 : #elif WITH_LAPACKE
19 : #include <lapacke.h>
20 : #endif
21 :
22 : namespace pairinteraction {
23 : #if defined(WITH_MKL) || defined(WITH_LAPACKE)
24 5 : lapack_int evr(int matrix_layout, char jobz, char range, char uplo, lapack_int n, float *a,
25 : lapack_int lda, float vl, float vu, lapack_int il, lapack_int iu, float abstol,
26 : lapack_int *m, float *w, float *z, lapack_int ldz, lapack_int *isuppz) {
27 5 : return LAPACKE_ssyevr(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w,
28 5 : z, ldz, isuppz);
29 : };
30 :
31 8 : lapack_int evr(int matrix_layout, char jobz, char range, char uplo, lapack_int n, double *a,
32 : lapack_int lda, double vl, double vu, lapack_int il, lapack_int iu, double abstol,
33 : lapack_int *m, double *w, double *z, lapack_int ldz, lapack_int *isuppz) {
34 8 : return LAPACKE_dsyevr(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w,
35 8 : z, ldz, isuppz);
36 : };
37 :
38 2 : lapack_int evr(int matrix_layout, char jobz, char range, char uplo, lapack_int n,
39 : lapack_complex_float *a, lapack_int lda, float vl, float vu, lapack_int il,
40 : lapack_int iu, float abstol, lapack_int *m, float *w, lapack_complex_float *z,
41 : lapack_int ldz, lapack_int *isuppz) {
42 2 : return LAPACKE_cheevr(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w,
43 2 : z, ldz, isuppz);
44 : };
45 :
46 3 : lapack_int evr(int matrix_layout, char jobz, char range, char uplo, lapack_int n,
47 : lapack_complex_double *a, lapack_int lda, double vl, double vu, lapack_int il,
48 : lapack_int iu, double abstol, lapack_int *m, double *w, lapack_complex_double *z,
49 : lapack_int ldz, lapack_int *isuppz) {
50 3 : return LAPACKE_zheevr(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w,
51 3 : z, ldz, isuppz);
52 : };
53 :
54 : template <typename Scalar>
55 : template <typename ScalarLim>
56 18 : EigenSystemH<Scalar> DiagonalizerLapackeEvr<Scalar>::dispatch_eigh(
57 : const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
58 : std::optional<real_t> min_eigenvalue, std::optional<real_t> max_eigenvalue, double rtol) const {
59 : using real_lim_t = typename traits::NumTraits<ScalarLim>::real_t;
60 18 : int dim = matrix.rows();
61 18 : bool has_range = min_eigenvalue.has_value() || max_eigenvalue.has_value();
62 :
63 : // Subtract the mean of the diagonal elements from the diagonal
64 18 : real_t shift{};
65 18 : Eigen::MatrixX<ScalarLim> shifted_matrix =
66 : this->template subtract_mean<ScalarLim>(matrix, shift, rtol);
67 18 : real_lim_t scaling = shifted_matrix.norm();
68 :
69 18 : if (scaling == 0) {
70 2 : scaling = 1; // Avoid division by zero if the matrix is zero
71 : }
72 :
73 18 : shifted_matrix /= scaling; // This seems to increase the numerical stability of lapacke_evr
74 :
75 : // Diagonalize the shifted matrix
76 18 : lapack_int m = 0; // Number of eigenvalues found
77 18 : char jobz = 'V'; // Compute eigenvectors
78 18 : char range_char = has_range ? 'V' : 'A'; // Compute all eigenvalues if no range is specified
79 18 : char uplo = 'U'; // Matrix is stored in upper-triangular part
80 18 : real_lim_t abstol = 0.1 * rtol; // 0.1 is a safety factor, abstol ~ rtol because ||H||=1
81 18 : lapack_int il = 0; // Lower index bound if 'I'
82 18 : lapack_int iu = 0; // Upper index bound if 'I'
83 18 : real_lim_t vl = (min_eigenvalue.value_or(-1) - shift) / scaling; // Lower eval bounds if 'V'
84 18 : real_lim_t vu = (max_eigenvalue.value_or(1) - shift) / scaling; // Upper eval bounds if 'V'
85 :
86 18 : Eigen::VectorX<real_lim_t> evals(dim); // Eigenvalues
87 18 : Eigen::MatrixX<ScalarLim> evecs(dim, dim); // Eigenvectors
88 18 : std::vector<lapack_int> isuppz(static_cast<size_t>(2 * dim)); // Workspace
89 18 : lapack_int info =
90 18 : evr(LAPACK_COL_MAJOR, jobz, range_char, uplo, dim, shifted_matrix.data(), dim, vl, vu, il,
91 : iu, abstol, &m, evals.data(), evecs.data(), dim, isuppz.data());
92 :
93 18 : if (info != 0) {
94 0 : if (info < 0) {
95 0 : throw std::invalid_argument(
96 : fmt::format("Diagonalization error: Argument {} to the "
97 : "lapacke_evr routine had an illegal value (the counting of the "
98 : "arguments starts with one). For a documentation of lapacke_evr, see "
99 : "https://www.intel.com/content/www/us/en/docs/onemkl/"
100 : "developer-reference-c/2025-1/syevr.html.",
101 : -info));
102 : }
103 0 : throw std::runtime_error(
104 : fmt::format("Diagonalization error: The lapacke_evr routine failed with error code {}. "
105 : "Try to set a lower 'rtol'.",
106 : info));
107 : }
108 :
109 : // Restrict to the first m eigenvectors and eigenvalues because the rest is not calculated
110 18 : evals.conservativeResize(m);
111 18 : evals *= scaling;
112 :
113 18 : return {evecs.leftCols(m).sparseView(1, 0.5 * rtol / std::sqrt(dim)).template cast<Scalar>(),
114 18 : this->add_mean(evals, shift)};
115 36 : }
116 :
117 : template <typename Scalar>
118 6 : DiagonalizerLapackeEvr<Scalar>::DiagonalizerLapackeEvr(FloatType float_type)
119 6 : : DiagonalizerInterface<Scalar>(float_type) {}
120 :
121 : template <typename Scalar>
122 : EigenSystemH<Scalar>
123 5 : DiagonalizerLapackeEvr<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
124 : double rtol) const {
125 5 : return this->eigh(matrix, std::nullopt, std::nullopt, rtol);
126 : }
127 :
128 : template <typename Scalar>
129 : EigenSystemH<Scalar>
130 18 : DiagonalizerLapackeEvr<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
131 : std::optional<real_t> min_eigenvalue,
132 : std::optional<real_t> max_eigenvalue, double rtol) const {
133 18 : switch (this->float_type) {
134 7 : case FloatType::FLOAT32:
135 : return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT32>>(
136 7 : matrix, min_eigenvalue, max_eigenvalue, rtol);
137 11 : case FloatType::FLOAT64:
138 : return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT64>>(
139 11 : matrix, min_eigenvalue, max_eigenvalue, rtol);
140 0 : default:
141 0 : throw std::invalid_argument("Unsupported floating point precision.");
142 : }
143 : }
144 :
145 : #else
146 :
147 : template <typename Scalar>
148 : DiagonalizerLapackeEvr<Scalar>::DiagonalizerLapackeEvr(FloatType float_type)
149 : : DiagonalizerInterface<Scalar>(float_type) {
150 : throw std::runtime_error("The lapacke_evr routine is not available in this build. Please use a "
151 : "different diagonalizer.");
152 : }
153 :
154 : template <typename Scalar>
155 : EigenSystemH<Scalar> DiagonalizerLapackeEvr<Scalar>::eigh(
156 : const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/, double /*rtol*/) const {
157 : std::abort(); // can't happen because the constructor throws
158 : }
159 :
160 : template <typename Scalar>
161 : EigenSystemH<Scalar> DiagonalizerLapackeEvr<Scalar>::eigh(
162 : const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/,
163 : std::optional<real_t> /*min_eigenvalue*/, std::optional<real_t> /*max_eigenvalue*/,
164 : double /*rtol*/) const {
165 : std::abort(); // can't happen because the constructor throws
166 : }
167 :
168 : #endif // WITH_MKL || WITH_LAPACKE
169 :
170 : // Explicit instantiations
171 : template class DiagonalizerLapackeEvr<double>;
172 : template class DiagonalizerLapackeEvr<std::complex<double>>;
173 : } // namespace pairinteraction
|
