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 0 : 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 0 : return LAPACKE_ssyevr(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w,
28 0 : 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 13 : 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 13 : int dim = matrix.rows();
61 13 : bool has_range = min_eigenvalue.has_value() || max_eigenvalue.has_value();
62 :
63 : // Subtract the mean of the diagonal elements from the diagonal
64 13 : real_t shift{};
65 13 : Eigen::MatrixX<ScalarLim> shifted_matrix =
66 : this->template subtract_mean<ScalarLim>(matrix, shift, rtol);
67 13 : real_lim_t scaling = shifted_matrix.norm();
68 :
69 13 : if (scaling < std::numeric_limits<real_lim_t>::epsilon()) {
70 2 : scaling = 1; // Avoid division by zero if the matrix is zero
71 : }
72 :
73 13 : shifted_matrix /= scaling; // This seems to increase the numerical stability of lapacke_evr
74 :
75 : // Diagonalize the shifted matrix
76 13 : lapack_int m = 0; // Number of eigenvalues found
77 13 : char jobz = 'V'; // Compute eigenvectors
78 13 : char range_char = has_range ? 'V' : 'A'; // Compute all eigenvalues if no range is specified
79 13 : char uplo = 'U'; // Matrix is stored in upper-triangular part
80 13 : real_lim_t abstol = 0.1 * rtol; // 0.1 is a safety factor, abstol ~ rtol because ||H||=1
81 13 : lapack_int il = 0; // Lower index bound if 'I'
82 13 : lapack_int iu = 0; // Upper index bound if 'I'
83 13 : real_lim_t vl = (min_eigenvalue.value_or(-1) - shift) / scaling; // Lower eval bounds if 'V'
84 13 : real_lim_t vu = (max_eigenvalue.value_or(1) - shift) / scaling; // Upper eval bounds if 'V'
85 :
86 13 : Eigen::VectorX<real_lim_t> evals(dim); // Eigenvalues
87 13 : Eigen::MatrixX<ScalarLim> evecs(dim, dim); // Eigenvectors
88 13 : std::vector<lapack_int> isuppz(static_cast<size_t>(2 * dim)); // Workspace
89 13 : lapack_int info =
90 13 : 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 13 : 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 13 : evals.conservativeResize(m);
111 13 : evals *= scaling;
112 :
113 13 : return {evecs.leftCols(m).sparseView(1, 0.5 * rtol / std::sqrt(dim)).template cast<Scalar>(),
114 13 : this->add_mean(evals, shift)};
115 26 : }
116 :
117 : template <typename Scalar>
118 4 : DiagonalizerLapackeEvr<Scalar>::DiagonalizerLapackeEvr(FloatType float_type)
119 4 : : DiagonalizerInterface<Scalar>(float_type) {}
120 :
121 : template <typename Scalar>
122 : EigenSystemH<Scalar>
123 0 : DiagonalizerLapackeEvr<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
124 : double rtol) const {
125 0 : return this->eigh(matrix, std::nullopt, std::nullopt, rtol);
126 : }
127 :
128 : template <typename Scalar>
129 : EigenSystemH<Scalar>
130 13 : 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 13 : switch (this->float_type) {
134 2 : case FloatType::FLOAT32:
135 : return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT32>>(
136 2 : 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
|
