zgesvd_tallskinny

Source: tests/dense/gesvd.cpp

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

// Containers.
#include <myramath/utility/Number.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/gesvd.h>
#include <myramath/dense/dimm.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Construct random input.
  auto A = Matrix<Number>::random(I,J);
  int K = I < J ? I : J;
  // Perform SVD and grab references to components.
  auto USV = gesvd(A);
  const auto& U = USV.U;
  const auto& S = USV.S;
  const auto& V = USV.V;
  // Check U*S*V-A
  Precision decomposition_error = frobenius(U*S*V-A);
  myra::out() << "  |U*S*V-A| = " << decomposition_error << std::endl;
  REQUIRE(decomposition_error < tolerance);
  // Check orthogonality of U
  Precision unitary_u = frobenius( gemm(U,'H',U) - Matrix<Number>::identity(K) );
  myra::out() << "  |U'U-I| = " << unitary_u << std::endl;
  REQUIRE(unitary_u < tolerance);
  // Check orthogonality of V.
  Precision unitary_v = frobenius( gemm(V,V,'H') - Matrix<Number>::identity(K) );
  myra::out() << "  |VV'-I| = " << unitary_v << std::endl;
  REQUIRE(unitary_u < tolerance);
  }

} // namespace

// Tall-skinny tests.
ADD_TEST("sgesvd_tallskinny","[dense][lapack]")
  { test<NumberS>(47,23,1.0e-4f); }

ADD_TEST("dgesvd_tallskinny","[dense][lapack]")
  { test<NumberD>(47,23,1.0e-12); }

ADD_TEST("cgesvd_tallskinny","[dense][lapack]")
  { test<NumberC>(47,23,1.0e-4f); }

ADD_TEST("zgesvd_tallskinny","[dense][lapack]")
  { test<NumberZ>(47,23,1.0e-12); }

// Short-fat tests.
ADD_TEST("sgesvd_shortfat","[dense][lapack]")
  { test<NumberS>(13,33,1.0e-4f); }

ADD_TEST("dgesvd_shortfat","[dense][lapack]")
  { test<NumberD>(13,33,1.0e-12); }

ADD_TEST("cgesvd_shortfat","[dense][lapack]")
  { test<NumberC>(13,33,1.0e-4f); }

ADD_TEST("zgesvd_shortfat","[dense][lapack]")
  { test<NumberZ>(13,33,1.0e-12); }

// Square tests.
ADD_TEST("sgesvd_square","[dense][lapack]")
  { test<NumberS>(25,25,1.0e-4f); }

ADD_TEST("dgesvd_square","[dense][lapack]")
  { test<NumberD>(25,25,1.0e-12); }

ADD_TEST("cgesvd_square","[dense][lapack]")
  { test<NumberC>(25,25,1.0e-4f); }

ADD_TEST("zgesvd_square","[dense][lapack]")
  { test<NumberZ>(25,25,1.0e-12); }

// Stress tests.
ADD_TEST("sgesvd_stress","[dense][lapack][.]")
  {
  int T = 100;
  for (int t = 0; t < T; ++t)
    test<NumberS>( random_int(1,50), random_int(1,50), 1.0e-4f );
  }

ADD_TEST("dgesvd_stress","[dense][lapack][.]")
  {
  int T = 100;
  for (int t = 0; t < T; ++t)
    test<NumberD>( random_int(1,50), random_int(1,50), 1.0e-12 );
  }

ADD_TEST("cgesvd_stress","[dense][lapack][.]")
  {
  int T = 100;
  for (int t = 0; t < T; ++t)
    test<NumberC>( random_int(1,50), random_int(1,50), 1.0e-4f );
  }

ADD_TEST("zgesvd_stress","[dense][lapack][.]")
  {
  int T = 100;
  for (int t = 0; t < T; ++t)
    test<NumberZ>( random_int(1,50), random_int(1,50), 1.0e-12 );
  }

Results: [PASS]

std::complex<double>
  |U*S*V-A| = 1.32584e-13
  |U'U-I| = 7.33825e-15
  |VV'-I| = 7.26709e-15

Go back to Summary of /test programs.