svd.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 "../blas/tools.hpp"
#include "../concepts.hpp"
#include "../declarations.hpp"
#include "../exceptions.hpp"
#include "../lapack/gesvd.hpp"
#include "../layout/policies.hpp"
#include "../macros.hpp"
#include "../mem/address_space.hpp"
#include "../mem/policies.hpp"
#include "../traits.hpp"
 
#include <algorithm>
#include <tuple>
#include <type_traits>
#include <utility>
 
namespace nda::linalg {


  template <blas_lapack::BlasArray<2> A>
  auto svd_in_place(A &&a) { // NOLINT (temporary views are allowed here)
    using layout_policy       = nda::detail::layout_to_policy<typename std::remove_cvref_t<A>::layout_t>::type;
    constexpr auto addr_space = mem::get_addr_space<A>;
 
    // vector s and matrices U and V^H
    auto const [m, n] = a.shape();
    auto s            = vector<get_fp_t<A>, heap<addr_space>>(std::min(m, n));
    auto U            = matrix<get_value_t<A>, layout_policy, heap<addr_space>>(m, m);
    auto VH           = matrix<get_value_t<A>, layout_policy, heap<addr_space>>(n, n);
 
    // call lapack gesvd
    int info = lapack::gesvd(a, s, U, VH);
    if (info != 0) NDA_RUNTIME_ERROR << "Error in nda::linalg::svd_in_place: gesvd returned a non-zero value: info = " << info;
 
    return std::make_tuple(U, s, VH);
  }

  template <Matrix A>
    requires(is_blas_lapack_v<get_value_t<A>>)
  auto svd(A const &a) { // NOLINT (temporary views are allowed here)
    return svd_in_place(basic_array{a});
  }

 
} // namespace nda::linalg