Permutation

Source: tests/sparse/Permutation.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/dense/Vector.h>
#include <myramath/dense/Matrix.h>
#include <myramath/sparse/Permutation.h>
#include <myramath/sparse/SparseMatrix.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/inverse.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/euclidean.h>
#include <myramath/sparse/gemm.h>
#include <myramath/sparse/gemv.h>
#include <myramath/sparse/frobenius.h>

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

#include <myramath/dense/intRange.h>

using namespace myra;

namespace {

// Tests gemm() operations involving Permutation's.
void test_Permutation1(int N)
  {
  // Define a random Permutation and its inverse.
  auto P = Permutation::random(N);
  P.verify();
  auto A = P.make_SparseMatrix<double>();

  // Define random Matrix X.
  auto X = Matrix<double>::random(N,N);

  // Check P*X = A*X
    {
    Matrix<double> PX = X; gemm_inplace(P,PX);
    Matrix<double> AX = gemm(A,X);
    double error = frobenius(PX-AX);
    myra::out() << "  |PX-AX| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  // Check P'*X = A'*X
    {
    Matrix<double> PtX = X; gemm_inplace(P,'T',PtX);
    Matrix<double> AtX = gemm(A,'T',X);
    double error = frobenius(PtX-AtX);
    myra::out() << "  |P'X-A'X| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  // Check X*P = X*A
    {
    Matrix<double> XP = X; gemm_inplace(XP,P);
    Matrix<double> XA = gemm(X,A);
    double error = frobenius(XP-XA);
    myra::out() << "  |XP-XA| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  // Check X*P' = X*A'
    {
    Matrix<double> XPt = X; gemm_inplace(XPt,P,'T');
    Matrix<double> XAt = gemm(X,A,'T');
    double error = frobenius(XPt-XAt);
    myra::out() << "  |XP'-XA'| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }
  
  }

// Tests gemv() operations involving Permutation's.
void test_Permutation2(int N)
  {

  // Define a random Permutation and its inverse.
  auto P = Permutation::random(N);
  P.verify();
  auto A = P.make_SparseMatrix<double>();
  
  // Define random Vector x.
  auto x = Vector<double>::random(N);

  // Check P*x = A*x
    {
    Vector<double> Px = x; gemv_inplace(P,Px);
    Vector<double> Ax = gemv(A,x);
    double error = frobenius(Px-Ax);
    myra::out() << "  |Px-Ax| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  // Check P'*x = A'*x
    {
    Vector<double> Ptx = x; gemv_inplace(P,'T',Ptx);
    Vector<double> Atx = gemv(A,'T',x);
    double error = frobenius(Ptx-Atx);
    myra::out() << "  |P'x-A'x| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  }

// Tests inverse(Permutation) operation.
void test_Permutation3(int N)
  {
  // Check that a Permutation and its inverse() lead to inverse matrices.
  auto P = Permutation::random(N);
  auto Q = inverse(P);
  P.verify();
  Q.verify();
  auto A = P.make_SparseMatrix<double>().make_Matrix();
  auto B = Q.make_SparseMatrix<double>().make_Matrix();
  auto I = Matrix<double>::identity(N);
  double error = frobenius( gemm(A,B) - I );
  myra::out() << "  |I - P*inverse(P)| = " << error << std::endl;
  REQUIRE(error < 1.0e-9);
  }

// Tests Permutation*Permutation composition.
void test_Permutation4(int N)
  {
  // Make random Permutations P and Q.
  auto P = Permutation::random(N);
  auto Q = Permutation::random(N);
  P.verify();
  Q.verify();
  auto PQ = multiply(P,Q);

  // Define random Vector x, copy to y.
  auto x = Vector<double>::random(N);
  auto y = x;

  // Compare P*(Q*x) ..
  gemv_inplace(Q,x); 
  gemv_inplace(P,x);
  // .. versus PQ*y
  gemv_inplace(PQ,y);
  double error = euclidean(x-y);
  myra::out() << "  |P*(Q*x) - (P*Q)*x| = " << error << std::endl;
  REQUIRE(error < 1.0e-9);
  }

} // namespace

ADD_TEST("Permutation","[sparse]")
  {
  test_Permutation1(20);
  test_Permutation2(20);
  test_Permutation3(20);
  test_Permutation4(20);
  }

Results: [PASS]

  |PX-AX| = 0
  |P'X-A'X| = 0
  |XP-XA| = 0
  |XP'-XA'| = 0
  |Px-Ax| = 0
  |P'x-A'x| = 0
  |I - P*inverse(P)| = 0
  |P*(Q*x) - (P*Q)*x| = 0

Go back to Summary of /test programs.