inverse_LowerMatrix

Source: tests/dense/inverse.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/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/inverse.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

// Generates DiagonalMatrix A, computes inv(A), checks A*inv(A)=I
template<class Number> void test_diagonal(int N, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Construct random matrix and compute its inverse.
  auto A = DiagonalMatrix<Number>::random(N);
  DiagonalMatrix<Number> inv_A = inverse(A);
  // Check A*inv_A = I
  Precision error = frobenius(A*inv_A-DiagonalMatrix<Number>::identity(N));
  myra::out() << "  |A*inv(A)-I| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

// Generates LowerMatrix A, computes inv(A), checks A*inv(A)=I
template<class Number> void test_lower(int N, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Construct random matrix and compute its inverse.
  auto A = LowerMatrix<Number>::random(N);
  // Make A well-conditioned by putting 1's on diagonal.
  for (int n = 0; n < N; ++n)
    A(n,n) = Number(1);
  LowerMatrix<Number> inv_A = inverse(A);
  // Check A*inv_A = I. Note that there's no LowerMatrix*LowerMatrix kernel, 
  // so have to promote both to explicit square Matrix's to multiply.
  Precision error = frobenius(A.make_Matrix()*inv_A.make_Matrix()-Matrix<Number>::identity(N));
  myra::out() << "  |A*inv(A)-I| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

// Generates square Matrix A, computes inv(A), checks A*inv(A)=I
template<class Number> void test_square(int N, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Construct random matrix and compute its inverse.
  auto A = Matrix<Number>::random(N,N);
  Matrix<Number> inv_A = inverse(A);  
  // Check A*inv_A = I
  Precision error = frobenius(A*inv_A-Matrix<Number>::identity(N));
  myra::out() << "  |A*inv(A)-I| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

} // namespace

ADD_TEST("inverse_DiagonalMatrix","[dense][lapack]")
  {
  myra::out() << "Testing explicit inverse routines (DiagonalMatrix)" << std::endl;
  test_diagonal<float >(15,1.0e-4f);
  test_diagonal<double>(16,1.0e-9);
  test_diagonal<std::complex<float > >(17,1.0e-4f);
  test_diagonal<std::complex<double> >(18,1.0e-9);
  myra::out() << std::endl;
  }

ADD_TEST("inverse_LowerMatrix","[dense][lapack]")
  {
  myra::out() << "Testing explicit inverse routines (LowerMatrix)" << std::endl;
  test_lower<float >(15,1.0e-4f);
  test_lower<double>(16,1.0e-9);
  test_lower<std::complex<float > >(17,1.0e-4f);
  test_lower<std::complex<double> >(18,1.0e-9);
  myra::out() << std::endl;
  }

ADD_TEST("inverse_Matrix","[dense][lapack]")
  {
  myra::out() << "Testing explicit inverse routines (square Matrix)" << std::endl;
  test_square<float >(15,1.0e-4f);
  test_square<double>(16,1.0e-9);
  test_square<std::complex<float > >(17,1.0e-4f);
  test_square<std::complex<double> >(18,1.0e-9);
  myra::out() << std::endl;  
  }

Results: [PASS]

Testing explicit inverse routines (LowerMatrix)
float
  |A*inv(A)-I| = 1.05797e-06
double
  |A*inv(A)-I| = 1.20657e-15
std::complex<float>
  |A*inv(A)-I| = 5.95333e-06
std::complex<double>
  |A*inv(A)-I| = 1.3808e-14

Go back to Summary of /test programs.