eigenelements.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2023 Simons Foundation
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0.txt
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Authors: Olivier Parcollet, Nils Wentzell
 
#pragma once
 
#include "./det_and_inverse.hpp"
#include "../basic_array.hpp"
#include "../declarations.hpp"
#include "../exceptions.hpp"
#include "../lapack/interface/cxx_interface.hpp"
#include "../layout/policies.hpp"
#include "../macros.hpp"
#include "../traits.hpp"
 
#include <type_traits>
#include <utility>
 
namespace nda::linalg {
 
  namespace detail {
 
    // Dispatch the call to the appropriate LAPACK routine based on the value type of the matrix.
    template <typename M>
    auto _eigen_element_impl(M &&m, char compz) { // NOLINT (temporary views are allowed here)
      using value_type = typename std::decay_t<M>::value_type;
 
      // runtime checks
      EXPECTS((not m.empty()));
      EXPECTS(is_matrix_square(m, true));
      EXPECTS(m.indexmap().is_contiguous());
 
      // set up the workspace
      int dim   = m.extent(0);
      int lwork = 64 * dim;
      array<double, 1> ev(dim);
      array<value_type, 1> work(lwork);
      array<double, 1> work2(is_complex_v<value_type> ? lwork : 0);
 
#if defined(__has_feature)
#if __has_feature(memory_sanitizer)
      work2 = 0;
      work  = 0;
      ev    = 0;
#endif
#endif
 
      // call the correct LAPACK routine
      int info = 0;
      if constexpr (not is_complex_v<value_type>) {
        lapack::f77::syev(compz, 'U', dim, m.data(), dim, ev.data(), work.data(), lwork, info);
      } else {
        lapack::f77::heev(compz, 'U', dim, m.data(), dim, ev.data(), work.data(), lwork, work2.data(), info);
      }
      if (info) NDA_RUNTIME_ERROR << "Error in nda::linalg::detail::_eigen_element_impl: Diagonalization error";
      return ev;
    }
 
  } // namespace detail
 
  template <typename M>
  auto eigenelements(M const &m) {
    auto m_copy = matrix<typename M::value_type, F_layout>(m);
    auto ev     = detail::_eigen_element_impl(m_copy, 'V');
    return std::pair<array<double, 1>, typename M::regular_type>{ev, m_copy};
  }
 
  template <typename M>
  auto eigenvalues(M const &m) {
    auto m_copy = matrix<typename M::value_type, F_layout>(m);
    return detail::_eigen_element_impl(m_copy, 'N');
  }
 
  template <typename M>
  auto eigenvalues_in_place(M &m) {
    return detail::_eigen_element_impl(m, 'N');
  }
 
} // namespace nda::linalg