Add svd and svd_in_place functions to compute SVDs

Author: Thoemi09

c++/nda/linalg.hpp

@@ -31,3 +31,4 @@
 #include "./linalg/matmul.hpp"
 #include "./linalg/norm.hpp"
 #include "./linalg/solve.hpp"
+#include "./linalg/svd.hpp"

c++/nda/linalg/svd.hpp

@@ -0,0 +1,112 @@
+// 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
+
+/**
+ * @file
+ * @brief Provides functions to compute the singular value decoomposition of a matrix.
+ */
+
+#pragma once
+
+#include "../basic_array.hpp"
+#include "../blas/tools.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>
+
+namespace nda {
+
+  /**
+   * @addtogroup linalg_tools
+   * @{
+   */
+
+  /**
+   * @brief Compute the singular value decomposition (SVD) of a matrix in place.
+   *
+   * @details The function computes the SVD of a given m-by-n matrix \f$ \mathbf{A} \f$:
+   * \f[
+   *   \mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^H \; ,
+   * \f]
+   * where \f$ \mathbf{U} \f$ is a unitary m-by-m matrix, \f$ \mathbf{V} \f$ is a unitary n-by-n matrix and \f$
+   * \mathbf{S} \f$ is an m-by-n matrix with non-negative real numbers on the diagonal.
+   *
+   * It first constructs the output vector \f$ \mathbf{s} \f$, which contains the singular values, and the output
+   * matrices \f$ \mathbf{U} \f$ and \f$ \mathbf{V}^H \f$. It then calls nda::lapack::gesvd to compute the SVD.
+   *
+   * @note If the input matrix \f$ \mathbf{A} \f$ is in Fortran layout, the output matrices \f$ \mathbf{U} \f$ and
+   * \f$ \mathbf{V}^H \f$ are also in Fortran layout. Otherwise, they are in C layout.
+   *
+   * @tparam A nda::MemoryMatrix type.
+   * @param a Input/output matrix. On entry, the m-by-n matrix \f$ \mathbf{A} \f$. On exit, the contents of \f$
+   * \mathbf{A} \f$ are destroyed.
+   * @return `std::tuple` containing \f$ \mathbf{U} \f$, \f$ \mathbf{s} \f$ and \f$ \mathbf{V}^H \f$.
+   */
+  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       = 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 s  = vector<double, heap<addr_space>>(std::min(a.extent(0), a.extent(1)));
+    auto U  = matrix<get_value_t<A>, layout_policy, heap<addr_space>>(a.extent(0), a.extent(0));
+    auto VH = matrix<get_value_t<A>, layout_policy, heap<addr_space>>(a.extent(1), a.extent(1));
+
+    // 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);
+  }
+
+  /**
+   * @brief Compute the singular value decomposition (SVD) of a matrix.
+   *
+   * @details The function computes the SVD of a given m-by-n matrix \f$ \mathbf{A} \f$:
+   * \f[
+   *   \mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^H \; ,
+   * \f]
+   * where \f$ \mathbf{U} \f$ is a unitary m-by-m matrix, \f$ \mathbf{V} \f$ is a unitary n-by-n matrix and \f$
+   * \mathbf{S} \f$ is an m-by-n matrix with non-negative real numbers on the diagonal.
+   *
+   * It first makes a copy of the input matrix \f$ \mathbf{A} \f$ and then calls nda::svd_in_place with the copy.
+   *
+   * @tparam A nda::MemoryMatrix type.
+   * @param a Input matrix \f$ \mathbf{A} \f$.
+   * @return `std::tuple` containing \f$ \mathbf{U} \f$, \f$ \mathbf{s} \f$ and \f$ \mathbf{V}^H \f$.
+   */
+  template <Matrix A>
+    requires(is_blas_lapack_v<get_value_t<A>>)
+  auto svd(A const &a) { // NOLINT (temporary views are allowed here)
+    using layout_policy       = detail::layout_to_policy<typename A::layout_t>::type;
+    constexpr auto addr_space = mem::get_addr_space<A>;
+    auto a_copy               = matrix<get_value_t<A>, layout_policy, heap<addr_space>>(a);
+    return svd_in_place(a_copy);
+  }
+
+  /** @} */
+
+} // namespace nda

test/c++/nda_linear_algebra.cpp

@@ -404,7 +404,7 @@ void test_solve() {
 
   // solve A * X = B using the exact matrix inverse
   auto Ainv = matrix_t{{-24, 18, 5}, {20, -15, -4}, {-5, 4, 1}};
-  auto X   = matrix_t{Ainv * B};
+  auto X    = matrix_t{Ainv * B};
   EXPECT_ARRAY_NEAR(matrix_t{A * X}, B);
 
   // solve A * X = B using solve_in_place
@@ -415,7 +415,7 @@ void test_solve() {
   EXPECT_ARRAY_NEAR(X, Bcopy);
 
   // solve A * x = b using solve_in_place
-  Acopy = A;
+  Acopy  = A;
   auto b = vector_t{B(nda::range::all, 0)};
   nda::solve_in_place(Acopy, b);
   EXPECT_ARRAY_NEAR(A * b, B(nda::range::all, 0));
@@ -438,3 +438,34 @@ TEST(NDA, LinearAlgebraSolve) {
   test_solve<std::complex<double>, nda::C_layout>();
   test_solve<std::complex<double>, nda::F_layout>();
 }
+
+// Check the SVD of a matrix.
+template <typename value_t, typename Layout>
+void test_svd() {
+  using matrix_t = nda::matrix<value_t, Layout>;
+
+  auto A = matrix_t{{2, -2, 1}, {-4, -8, -8}};
+  auto s = nda::vector<double>{12, 3};
+
+  // compute the SVD of A
+  auto [U_1, s_1, VH_1] = nda::svd(A);
+  auto S_1              = matrix_t::zeros(A.shape());
+  for (auto i : nda::range(2)) S_1(i, i) = s_1(i);
+  EXPECT_ARRAY_NEAR(s_1, s, 1e-14);
+  EXPECT_ARRAY_NEAR(A, U_1 * S_1 * VH_1, 1e-14);
+
+  // compute the SVD of A in place
+  auto A_copy           = A;
+  auto [U_2, s_2, VH_2] = nda::svd_in_place(A_copy);
+  auto S_2              = matrix_t::zeros(A.shape());
+  for (auto i : nda::range(2)) S_2(i, i) = s_2(i);
+  EXPECT_ARRAY_NEAR(s, s_2, 1e-14);
+  EXPECT_ARRAY_NEAR(A, U_2 * S_2 * VH_2, 1e-14);
+}
+
+TEST(NDA, LinearAlgebraSVD) {
+  test_svd<double, nda::C_layout>();
+  test_svd<double, nda::F_layout>();
+  test_svd<std::complex<double>, nda::C_layout>();
+  test_svd<std::complex<double>, nda::F_layout>();
+}