utils.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 "../basic_functions.hpp"
#include "../concepts.hpp"
#include "../declarations.hpp"
#include "../exceptions.hpp"
#include "../layout/policies.hpp"
#include "../macros.hpp"
#include "../mem/address_space.hpp"
#include "../traits.hpp"
 
#include <utility>
 
namespace nda::linalg {


  auto get_permutation_vector(Vector auto const &ipiv, int m) {
    static_assert(nda::mem::have_host_compatible_addr_space<decltype(ipiv)>);
    EXPECTS(m >= ipiv.size());
    auto sigma = vector<int>{arange<int>(m)};
    for (int i = 0; i < ipiv.size(); ++i) std::swap(sigma(i), sigma(ipiv(i) - 1));
    return sigma;
  }

  template <Scalar T, typename LP = F_layout>
  auto get_permutation_matrix(Vector auto const &sigma, bool column_permutations = false) {
    static_assert(nda::mem::have_host_compatible_addr_space<decltype(sigma)>);
    int m  = sigma.size();
    auto P = matrix<T, LP>::zeros(m, m);
    for (int i = 0; i < m; ++i) (column_permutations ? P(sigma(i), i) : P(i, sigma(i))) = T{1};
    return P;
  }

  template <Scalar T, typename LP = F_layout>
  auto get_permutation_matrix(Vector auto const &ipiv, int m) {
    return get_permutation_matrix<T, LP>(get_permutation_vector(ipiv, m));
  }

 
} // namespace nda::linalg