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