svd.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 "../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 <MemoryMatrix A>
    requires(is_blas_lapack_v<get_value_t<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 = nda::mem::get_addr_space<A>;
 
    // vector s and matrices U and V^H
    auto const [m, n] = a.shape();
    auto s            = vector<double, 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::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