Line data    Source code

       1             : // SPDX-FileCopyrightText: 2025 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/utils/tensor.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisPair.hpp"
       7             : #include "pairinteraction/utils/eigen_assertion.hpp"
       8             : #include "pairinteraction/utils/traits.hpp"
       9             : 
      10             : #include <Eigen/SparseCore>
      11             : #include <algorithm>
      12             : #include <limits>
      13             : #include <memory>
      14             : #include <oneapi/tbb.h>
      15             : 
      16             : namespace pairinteraction::utils {
      17             : 
      18             : template <typename Scalar>
      19             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
      20         233 : calculate_tensor_product(const std::shared_ptr<const BasisPair<Scalar>> &basis_initial,
      21             :                          const std::shared_ptr<const BasisPair<Scalar>> &basis_final,
      22             :                          const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix1,
      23             :                          const Eigen::SparseMatrix<Scalar, Eigen::RowMajor> &matrix2) {
      24             :     using real_t = typename traits::NumTraits<Scalar>::real_t;
      25         233 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
      26             : 
      27         233 :     oneapi::tbb::concurrent_vector<Eigen::Triplet<Scalar>> triplets;
      28             : 
      29             :     // Loop over the rows of the first matrix in parallel (outer index == row)
      30         233 :     oneapi::tbb::parallel_for(
      31       29907 :         oneapi::tbb::blocked_range<Eigen::Index>(0, matrix1.outerSize()), [&](const auto &range) {
      32       60974 :             for (Eigen::Index row1 = range.begin(); row1 != range.end(); ++row1) {
      33             : 
      34       31300 :                 const auto &range_row2 = basis_final->get_index_range(row1);
      35             : 
      36             :                 // Loop over the rows of the second matrix that are energetically allowed
      37      342368 :                 for (auto row2 = static_cast<Eigen::Index>(range_row2.min());
      38      342368 :                      row2 < static_cast<Eigen::Index>(range_row2.max()); ++row2) {
      39             : 
      40      311068 :                     Eigen::Index row = basis_final->get_ket_index_from_tuple(row1, row2);
      41      311068 :                     if (row < 0) {
      42      239565 :                         continue;
      43             :                     }
      44             : 
      45             :                     // Loop over the non-zero column elements of the first matrix
      46     1215155 :                     for (typename Eigen::SparseMatrix<Scalar, Eigen::RowMajor>::InnerIterator it1(
      47             :                              matrix1, row1);
      48     1215155 :                          it1; ++it1) {
      49             : 
      50     1143652 :                         Eigen::Index col1 = it1.col();
      51     1143652 :                         Scalar value1 = it1.value();
      52             : 
      53             :                         // Calculate the minimum and maximum values of the index pointer of the
      54             :                         // second matrix
      55     1143652 :                         Eigen::Index begin_idxptr2 = matrix2.outerIndexPtr()[row2];
      56     1143652 :                         Eigen::Index end_idxptr2 = matrix2.outerIndexPtr()[row2 + 1];
      57             : 
      58             :                         // The minimum value is chosen such that we start with an energetically
      59             :                         // allowed column
      60     1143652 :                         const auto &range_col2 = basis_initial->get_index_range(it1.index());
      61     1143652 :                         begin_idxptr2 +=
      62     1143652 :                             std::distance(matrix2.innerIndexPtr() + begin_idxptr2,
      63     1143652 :                                           std::lower_bound(matrix2.innerIndexPtr() + begin_idxptr2,
      64     1143652 :                                                            matrix2.innerIndexPtr() + end_idxptr2,
      65     1143652 :                                                            range_col2.min()));
      66             : 
      67             :                         // Loop over the non-zero column elements of the second matrix that are
      68             :                         // energetically allowed (we break the loop if the index pointer corresponds
      69             :                         // to a column that is not energetically allowed)
      70     2373926 :                         for (Eigen::Index idxptr2 = begin_idxptr2; idxptr2 < end_idxptr2;
      71             :                              ++idxptr2) {
      72             : 
      73     2166958 :                             Eigen::Index col2 = matrix2.innerIndexPtr()[idxptr2];
      74     2166958 :                             if (col2 >= static_cast<Eigen::Index>(range_col2.max())) {
      75      936684 :                                 break;
      76             :                             }
      77             : 
      78     1230274 :                             Eigen::Index col = basis_initial->get_ket_index_from_tuple(col1, col2);
      79     1230274 :                             if (col < 0) {
      80       38250 :                                 continue;
      81             :                             }
      82             : 
      83     1192024 :                             Scalar value2 = matrix2.valuePtr()[idxptr2];
      84             : 
      85             :                             // Store the entry
      86           0 :                             Scalar value = value1 * value2;
      87     1192024 :                             if (std::abs(value) > numerical_precision) {
      88     1191148 :                                 triplets.emplace_back(row, col, value);
      89             :                             }
      90             :                         }
      91             :                     }
      92             :                 }
      93             :             }
      94             :         });
      95             : 
      96             :     // Construct the combined matrix from the triplets
      97         233 :     Eigen::SparseMatrix<Scalar, Eigen::RowMajor> matrix(basis_final->get_number_of_kets(),
      98         233 :                                                         basis_initial->get_number_of_kets());
      99         233 :     matrix.setFromTriplets(triplets.begin(), triplets.end());
     100         233 :     matrix.makeCompressed();
     101             : 
     102         466 :     return matrix;
     103         233 : }
     104             : 
     105             : // Explicit instantiations
     106             : template Eigen::SparseMatrix<double, Eigen::RowMajor>
     107             : calculate_tensor_product(const std::shared_ptr<const BasisPair<double>> &,
     108             :                          const std::shared_ptr<const BasisPair<double>> &,
     109             :                          const Eigen::SparseMatrix<double, Eigen::RowMajor> &,
     110             :                          const Eigen::SparseMatrix<double, Eigen::RowMajor> &);
     111             : template Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor>
     112             : calculate_tensor_product(const std::shared_ptr<const BasisPair<std::complex<double>>> &,
     113             :                          const std::shared_ptr<const BasisPair<std::complex<double>>> &,
     114             :                          const Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> &,
     115             :                          const Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor> &);
     116             : 
     117             : } // namespace pairinteraction::utils