tools.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 "../concepts.hpp"
#include "../map.hpp"
#include "../mapped_functions.hpp"
 
#include <complex>
#include <tuple>
#include <type_traits>
#include <utility>
 
namespace nda {

  using dcomplex = std::complex<double>;
 
} // namespace nda
 
namespace nda::blas {


  template <typename A>
  static constexpr bool is_conj_array_expr = false;

  template <MemoryArray A>
  static constexpr bool is_conj_array_expr<expr_call<detail::conj_f, A>> = true;
 
  // Specialization of nda::blas::is_conj_array_expr for cvref types.
  template <typename A>
    requires(!std::is_same_v<A, std::remove_cvref_t<A>>)
  static constexpr bool is_conj_array_expr<A> = is_conj_array_expr<std::remove_cvref_t<A>>;

  template <Array A>
    requires(MemoryArray<A> or is_conj_array_expr<A>)
  MemoryArray decltype(auto) get_array(A &&a) {
    if constexpr (is_conj_array_expr<A>) {
      return std::get<0>(std::forward<A>(a).a);
    } else {
      return std::forward<A>(a);
    }
  }

  template <Array A>
    requires(MemoryArray<A> or is_conj_array_expr<A>)
  static constexpr bool has_F_layout = []() {
    if constexpr (is_conj_array_expr<A>)
      return has_F_layout<decltype(std::get<0>(std::declval<A>().a))>;
    else
      return std::remove_cvref_t<A>::is_stride_order_Fortran();
  }();

  template <Array A>
    requires(MemoryArray<A> or is_conj_array_expr<A>)
  static constexpr bool has_C_layout = []() {
    if constexpr (is_conj_array_expr<A>)
      return has_C_layout<decltype(std::get<0>(std::declval<A>().a))>;
    else
      return std::remove_cvref_t<A>::is_stride_order_C();
  }();

  template <Array A>
  static constexpr char get_op = []() {
    auto constexpr conj      = is_conj_array_expr<A>;
    auto constexpr transpose = has_C_layout<A>;
    static_assert(!(conj and not transpose), "Error in nda::blas::get_op: Cannot use conjugate operation alone in BLAS operations");
    if constexpr (conj and transpose)
      return 'C';
    else if constexpr (transpose)
      return 'T';
    else // !conj and !transpose
      return 'N';
  }();

  template <MemoryArray A>
    requires(get_rank<A> == 1 or get_rank<A> == 2)
  int get_ld(A const &a) {
    if constexpr (get_rank<A> == 1) {
      return a.size();
    } else {
      return a.indexmap().strides()[has_F_layout<A> ? 1 : 0];
    }
  }

  template <MemoryArray A>
    requires(get_rank<A> == 1 or get_rank<A> == 2)
  int get_ncols(A const &a) {
    if constexpr (get_rank<A> == 1) {
      return 1;
    } else {
      return a.shape()[has_F_layout<A> ? 1 : 0];
    }
  }

 
} // namespace nda::blas