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