Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvd.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 <spdlog/spdlog.h>
      14             : 
      15             : #ifdef WITH_MKL
      16             : #include <mkl.h>
      17             : #elif WITH_LAPACKE
      18             : #include <lapacke.h>
      19             : #endif
      20             : 
      21             : namespace pairinteraction {
      22             : #if defined(WITH_MKL) || defined(WITH_LAPACKE)
      23           0 : lapack_int evd(int matrix_layout, char jobz, char uplo, lapack_int n, float *a, lapack_int lda,
      24             :                float *w) {
      25           0 :     return LAPACKE_ssyevd(matrix_layout, jobz, uplo, n, a, lda, w);
      26             : };
      27             : 
      28           2 : lapack_int evd(int matrix_layout, char jobz, char uplo, lapack_int n, double *a, lapack_int lda,
      29             :                double *w) {
      30           2 :     return LAPACKE_dsyevd(matrix_layout, jobz, uplo, n, a, lda, w);
      31             : };
      32             : 
      33           2 : lapack_int evd(int matrix_layout, char jobz, char uplo, lapack_int n, lapack_complex_float *a,
      34             :                lapack_int lda, float *w) {
      35           2 :     return LAPACKE_cheevd(matrix_layout, jobz, uplo, n, a, lda, w);
      36             : };
      37             : 
      38           3 : lapack_int evd(int matrix_layout, char jobz, char uplo, lapack_int n, lapack_complex_double *a,
      39             :                lapack_int lda, double *w) {
      40           3 :     return LAPACKE_zheevd(matrix_layout, jobz, uplo, n, a, lda, w);
      41             : };
      42             : 
      43             : template <typename Scalar>
      44             : template <typename ScalarLim>
      45           7 : EigenSystemH<Scalar> DiagonalizerLapackeEvd<Scalar>::dispatch_eigh(
      46             :     const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix, double rtol) const {
      47             :     using real_lim_t = typename traits::NumTraits<ScalarLim>::real_t;
      48           7 :     int dim = matrix.rows();
      49             : 
      50             :     // Subtract the mean of the diagonal elements from the diagonal
      51           7 :     real_t shift{};
      52           7 :     Eigen::MatrixX<ScalarLim> evecs = this->template subtract_mean<ScalarLim>(matrix, shift, rtol);
      53             : 
      54             :     // Diagonalize the shifted matrix
      55           7 :     char jobz = 'V'; // eigenvalues and eigenvectors are computed
      56           7 :     char uplo = 'U'; // full matrix is stored, upper is used
      57             : 
      58           7 :     Eigen::VectorX<real_lim_t> evals(dim);
      59           7 :     lapack_int info = evd(LAPACK_COL_MAJOR, jobz, uplo, dim, evecs.data(), dim, evals.data());
      60             : 
      61           7 :     if (info != 0) {
      62           0 :         if (info < 0) {
      63           0 :             throw std::invalid_argument(
      64             :                 fmt::format("Diagonalization error: Argument {} to the "
      65             :                             "lapacke_evd routine had an illegal value (the counting of the "
      66             :                             "arguments starts with one). For a documentation of lapacke_evd, see "
      67             :                             "https://www.intel.com/content/www/us/en/docs/onemkl/"
      68             :                             "developer-reference-c/2025-1/syevd.html.",
      69             :                             -info));
      70             :         }
      71           0 :         throw std::runtime_error(fmt::format(
      72             :             "Diagonalization error: The lapacke_evd routine failed with error code {}.", info));
      73             :     }
      74             : 
      75           7 :     return {evecs.sparseView(1, 0.5 * rtol / std::sqrt(dim)).template cast<Scalar>(),
      76           7 :             this->add_mean(evals, shift)};
      77          14 : }
      78             : 
      79             : template <typename Scalar>
      80           4 : DiagonalizerLapackeEvd<Scalar>::DiagonalizerLapackeEvd(FloatType float_type)
      81           4 :     : DiagonalizerInterface<Scalar>(float_type) {}
      82             : 
      83             : template <typename Scalar>
      84             : EigenSystemH<Scalar>
      85           7 : DiagonalizerLapackeEvd<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
      86             :                                      double rtol) const {
      87           7 :     switch (this->float_type) {
      88           2 :     case FloatType::FLOAT32:
      89           2 :         return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT32>>(matrix, rtol);
      90           5 :     case FloatType::FLOAT64:
      91           5 :         return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT64>>(matrix, rtol);
      92           0 :     default:
      93           0 :         throw std::invalid_argument("Unsupported floating point precision.");
      94             :     }
      95             : }
      96             : 
      97             : #else
      98             : 
      99             : template <typename Scalar>
     100             : DiagonalizerLapackeEvd<Scalar>::DiagonalizerLapackeEvd(FloatType float_type)
     101             :     : DiagonalizerInterface<Scalar>(float_type) {
     102             :     throw std::runtime_error("The lapacke_evd routine is not available in this build. Please use a "
     103             :                              "different diagonalizer.");
     104             : }
     105             : 
     106             : template <typename Scalar>
     107             : EigenSystemH<Scalar> DiagonalizerLapackeEvd<Scalar>::eigh(
     108             :     const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> & /*matrix*/, double /*rtol*/) const {
     109             :     std::abort(); // can't happen because the constructor throws
     110             : }
     111             : 
     112             : #endif // WITH_MKL || WITH_LAPACKE
     113             : 
     114             : // Explicit instantiations
     115             : template class DiagonalizerLapackeEvd<double>;
     116             : template class DiagonalizerLapackeEvd<std::complex<double>>;
     117             : } // namespace pairinteraction