matrix_functions.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 "./accessors.hpp"
#include "./concepts.hpp"
#include "./declarations.hpp"
#include "./layout/policies.hpp"
#include "./layout_transforms.hpp"
#include "./macros.hpp"
#include "./mapped_functions.hpp"
#include "./mem/policies.hpp"
#include "./stdutil/array.hpp"
 
#include <algorithm>
#include <concepts>
#include <iostream>
#include <ranges>
#include <type_traits>
 
namespace nda {


  template <Scalar S = double, std::integral Int = long>
  auto eye(Int dim) {
    auto r = matrix<S>(dim, dim);
    r      = S{1};
    return r;
  }

  template <ArrayOfRank<2> M>
  auto trace(M const &m) {
    static_assert(get_rank<M> == 2, "Error in nda::trace: Array/View must have rank 2");
    EXPECTS(m.shape()[0] == m.shape()[1]);
    auto r = get_value_t<M>{};
    auto d = m.shape()[0];
    for (int i = 0; i < d; ++i) r += m(i, i);
    return r;
  }

  template <ArrayOfRank<2> M>
  ArrayOfRank<2> auto dagger(M const &m) {
    if constexpr (is_complex_v<typename M::value_type>)
      return conj(transpose(m));
    else
      return transpose(m);
  }

  template <MemoryArrayOfRank<2> M>
  ArrayOfRank<1> auto diagonal(M &&m) { // NOLINT
    long dim    = std::min(m.shape()[0], m.shape()[1]);
    long stride = stdutil::sum(m.indexmap().strides());
    using vector_view_t =
       basic_array_view<std::remove_reference_t<decltype(*m.data())>, 1, C_stride_layout, 'V', nda::default_accessor, nda::borrowed<>>;
    return vector_view_t{C_stride_layout::mapping<1>{{dim}, {stride}}, m.data()};
  }

  template <typename V>
    requires(std::ranges::contiguous_range<V> or ArrayOfRank<V, 1>)
  ArrayOfRank<2> auto diag(V const &v) {
    if constexpr (std::ranges::contiguous_range<V>) {
      return diag(nda::basic_array_view{v});
    } else {
      auto m      = matrix<std::remove_const_t<typename V::value_type>>::zeros({v.size(), v.size()});
      diagonal(m) = v;
      return m;
    }
  }

  template <ArrayOfRank<2> A, ArrayOfRank<2> B>
    requires(std::same_as<get_value_t<A>, get_value_t<B>>) // NB the get_value_t gets rid of const if any
  matrix<get_value_t<A>> vstack(A const &a, B const &b) {
    static_assert(get_rank<A> == 2, "Error in nda::vstack: Only rank 2 arrays/views are allowed");
    static_assert(get_rank<A> == 2, "Error in nda::vstack: Only rank 2 arrays/views are allowed");
    EXPECTS_WITH_MESSAGE(a.shape()[1] == b.shape()[1], "Error in nda::vstack: The second dimension of the two matrices must be equal");
 
    auto [n, q] = a.shape();
    auto p      = b.shape()[0];
    matrix<get_value_t<A>> res(n + p, q);
    res(range(0, n), range::all)     = a;
    res(range(n, n + p), range::all) = b;
    return res;
  }

  template <Matrix A>
  bool is_matrix_square(A const &a, bool print_error = false) {
    auto const [m, n] = a.shape();
    if (m != n and print_error) std::cerr << "Error in nda::is_matrix_square: Dimensions are: (" << m << "," << n << ")\n" << std::endl;
    return m == n;
  }

  template <Matrix A>
  bool is_matrix_diagonal(A const &a, bool print_error = false) {
    bool const r = is_matrix_square(a) and a == diag(diagonal(a));
    if (not r and print_error) std::cerr << "Error in nda::is_matrix_diagonal: Non-diagonal matrix: " << a << std::endl;
    return r;
  }

 
} // namespace nda