drank

Source: tests/dense/rank.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>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/nullspace.h>
#include <myramath/dense/rangespace.h>

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, int K, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // From IxJ sized matrix with rank K.
  auto U = Matrix<Number>::random(I,K);
  auto V = Matrix<Number>::random(K,J);
  auto A = U*V;
  auto Ah = hermitian(A);
  // Extract nullspace/rangespace of A.
  auto null_A = nullspace(A);
  auto range_A = rangespace(A);
  int rank = range_A.size().second;
  int nullity = null_A.size().second;
  // Extract nullspace of A'
  auto null_Ah = nullspace(Ah);
  auto range_Ah = rangespace(Ah);
  // Verify nullspaces form orthonormal bases.
  Precision error1 = frobenius(gemm(null_A,'H',null_A) - Matrix<Number>::identity(nullity));
  myra::out() << "  |null(A)'null(A)-I| = " << error1 << std::endl;
  Precision error2 = frobenius(gemm(null_Ah,'H',null_Ah) - Matrix<Number>::identity(nullity));  
  myra::out() << "  |null(A')'null(A')-I| = " << error2 << std::endl;
  // Verify A*null(A) and A'*null(A') are vanishingly small.
  Precision error3 = frobenius(gemm(A,null_A));
  myra::out() << "  |A*null(A)| = " << error3 << std::endl;
  Precision error4 = frobenius(gemm(Ah,null_Ah));
  myra::out() << "  |A'*null(A')| = " << error4 << std::endl;  
  // Verify rangespaces form orthonormal basis.
  Precision error5 = frobenius(gemm(range_A,'H',range_A) - Matrix<Number>::identity(rank));
  myra::out() << "  |range(A)'range(A)-I| = " << error5 << std::endl;
  Precision error6 = frobenius(gemm(range_Ah,'H',range_Ah) - Matrix<Number>::identity(rank));  
  myra::out() << "  |range(A')'range(A')-I| = " << error6 << std::endl;
  // Verify range(A')'*null(A) and range(A)*null(A') are vanishingly small.
  Precision error7 = frobenius( gemm(range_Ah,'H',null_A) );
  myra::out() << "  |range(A')*null(A)| = " << error7 << std::endl;
  Precision error8 = frobenius( gemm(range_A,'H',null_Ah) );
  myra::out() << "  |range(A)*null(A')| = " << error8 << std::endl;
  // Assert everything.
  REQUIRE(error1 < tolerance);
  REQUIRE(error2 < tolerance);
  REQUIRE(error3 < tolerance);
  REQUIRE(error4 < tolerance);
  REQUIRE(error5 < tolerance);
  REQUIRE(error6 < tolerance);
  REQUIRE(error7 < tolerance);
  REQUIRE(error8 < tolerance);  
  }

} // namespace

ADD_TEST("srank","[dense][lapack]")
  { test<NumberS>(20,17,12,1.0e-4f); }

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

ADD_TEST("crank","[dense][lapack]")
  { test<NumberC>(20,17,12,1.0e-4f); }

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

Results: [PASS]

double
  |null(A)'null(A)-I| = 1.78696e-15
  |null(A')'null(A')-I| = 1.98733e-15
  |A*null(A)| = 7.68344e-15
  |A'*null(A')| = 4.60185e-14
  |range(A)'range(A)-I| = 2.86933e-15
  |range(A')'range(A')-I| = 3.06559e-15
  |range(A')*null(A)| = 3.11091e-15
  |range(A)*null(A')| = 9.4904e-14

Go back to Summary of /test programs.