eigh.hpp

Go to the documentation of this file.

// Copyright (c) 2019--present, The Simons Foundation
// This file is part of TRIQS/nda and is licensed under the Apache License, Version 2.0.
// SPDX-License-Identifier: Apache-2.0
// See LICENSE in the root of this distribution for details.

 
#pragma once
 
#include "../basic_array.hpp"
#include "../blas/tools.hpp"
#include "../concepts.hpp"
#include "../declarations.hpp"
#include "../exceptions.hpp"
#include "../lapack/syev.hpp"
#include "../lapack/sygv.hpp"
#include "../lapack/heev.hpp"
#include "../lapack/hegv.hpp"
#include "../layout/policies.hpp"
#include "../macros.hpp"
#include "../matrix_functions.hpp"
#include "../mem/address_space.hpp"
#include "../traits.hpp"
 
#include <complex>
#include <type_traits>
#include <utility>
 
namespace nda::linalg {

 
  namespace detail {
 
    // Perform the call to the LAPACK routines syev/heev for eigh_in_place and eigvalsh_in_place.
    template <typename A>
    auto eigh_impl(A &&a, char jobz) { // NOLINT (temporary views are allowed here)
      // early return if the matrix is empty
      if (a.empty()) return array<double, 1>{};
 
      // make the call to syev/heev
      auto lambda = array<double, 1>(a.extent(0));
      int info    = 0;
      if constexpr (is_complex_v<get_value_t<A>>) {
        info = nda::lapack::heev(a, lambda, jobz);
      } else {
        info = nda::lapack::syev(a, lambda, jobz);
      }
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::detail::eigh_impl: syev/heev routine failed: info = " << info;
 
      return lambda;
    }
 
    // Perform the call to the LAPACK routines sygv/hegv for eigh_in_place and eigvalsh_in_place.
    template <typename A, typename B>
    auto eigh_impl(A &&a, B &&b, char jobz, int itype) { // NOLINT (temporary views are allowed here)
      // early return if the matrix is empty
      if (a.empty()) return array<double, 1>{};
 
      // make the call to sygv/hegv
      auto lambda = array<double, 1>(a.extent(0));
      int info    = 0;
      if constexpr (is_complex_v<get_value_t<A>>) {
        info = nda::lapack::hegv(a, b, lambda, jobz, itype);
      } else {
        info = nda::lapack::sygv(a, b, lambda, jobz, itype);
      }
      if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::detail::eigh_impl: sygv/hegv routine failed: info = " << info;
 
      return lambda;
    }
 
  } // namespace detail

  template <MemoryMatrix A>
    requires(nda::mem::have_host_compatible_addr_space<A> and is_blas_lapack_v<get_value_t<A>> and nda::blas::has_F_layout<A>)
  auto eigh_in_place(A &&a) {
    return detail::eigh_impl(std::forward<A>(a), 'V');
  }

  template <MemoryMatrix A, MemoryMatrix B>
    requires(nda::mem::have_host_compatible_addr_space<A, B> and is_blas_lapack_v<get_value_t<A>> and have_same_value_type_v<A, B>
             and nda::blas::has_F_layout<A> and nda::blas::has_F_layout<B>)
  auto eigh_in_place(A &&a, B &&b, int itype = 1) {
    return detail::eigh_impl(std::forward<A>(a), std::forward<B>(b), 'V', itype);
  }

  template <Matrix A>
    requires(Scalar<get_value_t<A>>)
  auto eigh(A const &a) {
    using value_t = std::conditional_t<is_complex_v<get_value_t<A>>, std::complex<double>, double>;
    auto a_copy   = matrix<value_t, F_layout>{a};
    auto lambda   = eigh_in_place(a_copy);
    return std::make_pair(lambda, a_copy);
  }

  template <Matrix A, Matrix B>
    requires(Scalar<get_value_t<A>> and Scalar<get_value_t<B>>)
  auto eigh(A const &a, B const &b, int itype = 1) {
    using value_t = std::conditional_t<is_complex_v<get_value_t<A>> or is_complex_v<get_value_t<B>>, std::complex<double>, double>;
    auto a_copy   = matrix<value_t, F_layout>{a};
    auto b_copy   = matrix<value_t, F_layout>{b};
    auto lambda   = eigh_in_place(a_copy, b_copy, itype);
    return std::make_pair(lambda, a_copy);
  }

  template <MemoryMatrix A>
    requires(nda::mem::have_host_compatible_addr_space<A> and is_blas_lapack_v<get_value_t<A>> and nda::blas::has_F_layout<A>)
  auto eigvalsh_in_place(A &&a) {
    return detail::eigh_impl(std::forward<A>(a), 'N');
  }

  template <MemoryMatrix A, MemoryMatrix B>
    requires(nda::mem::have_host_compatible_addr_space<A, B> and is_blas_lapack_v<get_value_t<A>> and have_same_value_type_v<A, B>
             and nda::blas::has_F_layout<A> and nda::blas::has_F_layout<B>)
  auto eigvalsh_in_place(A &&a, B &&b, int itype = 1) {
    return detail::eigh_impl(std::forward<A>(a), std::forward<B>(b), 'N', itype);
  }

  template <Matrix A>
    requires(Scalar<get_value_t<A>>)
  auto eigvalsh(A const &a) {
    using value_t = std::conditional_t<is_complex_v<get_value_t<A>>, std::complex<double>, double>;
    auto a_copy   = matrix<value_t, F_layout>{a};
    return eigvalsh_in_place(a_copy);
  }

  template <Matrix A, Matrix B>
    requires(Scalar<get_value_t<A>> and Scalar<get_value_t<B>>)
  auto eigvalsh(A const &a, B const &b, int itype = 1) {
    using value_t = std::conditional_t<is_complex_v<get_value_t<A>> or is_complex_v<get_value_t<B>>, std::complex<double>, double>;
    auto a_copy   = matrix<value_t, F_layout>{a};
    auto b_copy   = matrix<value_t, F_layout>{b};
    return eigvalsh_in_place(a_copy, b_copy, itype);
  }

 
} // namespace nda::linalg