utils.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024 Simons Foundation
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0.txt
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Authors: Thomas Hahn, Olivier Parcollet, Nils Wentzell

 
#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