outer_product.hpp

Go to the documentation of this file.

// Copyright (c) 2020--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 "../blas/ger.hpp"
#include "../blas/tools.hpp"
#include "../concepts.hpp"
#include "../declarations.hpp"
#include "../macros.hpp"
#include "../mem/address_space.hpp"
#include "../mem/policies.hpp"
#include "../stdutil/array.hpp"
#include "../traits.hpp"
 
#include <array>
#include <type_traits>
 
namespace nda::linalg {

  template <blas_lapack::BlasArray A, blas_lapack::BlasArrayFor<A> B>
    requires(blas_lapack::has_C_layout<A, B> or blas_lapack::has_F_layout<A, B>)
  auto outer_product(A const &a, B const &b) {
    // check the input arrays/views
    EXPECTS(a.is_contiguous());
    EXPECTS(b.is_contiguous());
 
    // get the return type
    auto constexpr rank    = get_rank<A> + get_rank<B>;
    auto constexpr algebra = []() {
      if constexpr (get_algebra<A> == 'V' and get_algebra<B> == 'V') {
        return 'M';
      } else {
        return 'A';
      }
    }();
    using layout_pol = std::conditional_t<get_rank<A> == 1, typename B::layout_policy_t::contiguous_t, typename A::layout_policy_t::contiguous_t>;
    using cont_pol   = heap<mem::common_addr_space<A, B>>;
    using return_t   = basic_array<get_value_t<A>, rank, layout_pol, algebra, cont_pol>;
 
    // use ger to calculate the outer product
    auto res   = return_t::zeros(stdutil::join(a.shape(), b.shape()));
    auto a_vec = reshape(a, std::array{a.size()});
    auto b_vec = reshape(b, std::array{b.size()});
    auto mat   = reshape(res, std::array{a.size(), b.size()});
    blas::ger(1.0, a_vec, b_vec, mat);
 
    return res;
  }
 
} // namespace nda::linalg