Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
       5             : 
       6             : #include "pairinteraction/enums/FloatType.hpp"
       7             : #include "pairinteraction/utils/eigen_assertion.hpp"
       8             : #include "pairinteraction/utils/eigen_compat.hpp"
       9             : #include "pairinteraction/utils/traits.hpp"
      10             : 
      11             : #include <Eigen/Dense>
      12             : #include <Eigen/Eigenvalues>
      13             : #include <cmath>
      14             : 
      15             : namespace pairinteraction {
      16             : 
      17             : template <typename Scalar>
      18          15 : DiagonalizerEigen<Scalar>::DiagonalizerEigen(FloatType float_type)
      19          15 :     : DiagonalizerInterface<Scalar>(float_type) {}
      20             : 
      21             : template <typename Scalar>
      22             : template <typename ScalarLim>
      23             : EigenSystemH<Scalar>
      24         121 : DiagonalizerEigen<Scalar>::dispatch_eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
      25             :                                          double rtol) const {
      26             :     using real_t = typename traits::NumTraits<Scalar>::real_t;
      27         121 :     int dim = matrix.rows();
      28             : 
      29             :     // Subtract the mean of the diagonal elements from the diagonal
      30         121 :     real_t shift{};
      31         121 :     Eigen::MatrixX<ScalarLim> shifted_matrix =
      32             :         this->template subtract_mean<ScalarLim>(matrix, shift, rtol);
      33             : 
      34             :     // Diagonalize the shifted matrix
      35         120 :     Eigen::SelfAdjointEigenSolver<Eigen::MatrixX<ScalarLim>> eigensolver;
      36         119 :     eigensolver.compute(shifted_matrix);
      37             : 
      38         242 :     return {eigensolver.eigenvectors()
      39         121 :                 .sparseView(1, 0.5 * rtol / std::sqrt(dim))
      40          53 :                 .template cast<Scalar>(),
      41         121 :             this->add_mean(eigensolver.eigenvalues(), shift)};
      42         239 : }
      43             : 
      44             : template <typename Scalar>
      45             : EigenSystemH<Scalar>
      46         121 : DiagonalizerEigen<Scalar>::eigh(const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix,
      47             :                                 double rtol) const {
      48         121 :     switch (this->float_type) {
      49          67 :     case FloatType::FLOAT32:
      50          67 :         return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT32>>(matrix, rtol);
      51          54 :     case FloatType::FLOAT64:
      52          54 :         return dispatch_eigh<traits::restricted_t<Scalar, FloatType::FLOAT64>>(matrix, rtol);
      53           0 :     default:
      54           0 :         throw std::invalid_argument("Unsupported floating point precision.");
      55             :     }
      56             : }
      57             : 
      58             : // Explicit instantiations
      59             : template class DiagonalizerEigen<double>;
      60             : template class DiagonalizerEigen<std::complex<double>>;
      61             : } // namespace pairinteraction