swaps

Source: tests/dense/swaps.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/Matrix.h>

// Algorithms.
#include <myramath/utility/ilinspace.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/swaps.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/intRange.h>

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

#include <algorithm>

using namespace myra;

namespace {

// Tests various operations w/ swap sequences.
void test_swaps(int I)
  {
  // Define a random permutation perm[] ..
  std::vector<int> perm = ilinspace(I);
  std::random_shuffle(perm.begin(),perm.end());
  // .. transform it into a swap sequence ..
  std::vector<int> swaps = perm2swaps(perm);
  // .. and also build it into a Matrix.
  Matrix<double> P = swaps2Matrix<double>(swaps);

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

  // Check P*X = swap_rows(X)
    {
    Matrix<double> PX = gemm(P,X);
    Matrix<double> X_copy = X;
    swap_rows(swaps,X_copy);
    double error = frobenius(PX-X_copy);
    myra::out() << "  |P*X-swap_rows(X)| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }
    
  // Check P'*X = iswap_rows(X)
    {
    Matrix<double> PtX = gemm(P,'T',X);
    Matrix<double> X_copy = X;
    iswap_rows(swaps,X_copy);
    double error = frobenius(PtX-X_copy);
    myra::out() << "  |P'*X-iswap_rows(X)| = " << error << std::endl;  
    REQUIRE(error < 1.0e-9);
    }
  
  // Check that swap_rows() and iswap_rows() are reverse operations.
    {
    Matrix<double> X_copy = X;
     swap_rows(swaps,X_copy);
    iswap_rows(swaps,X_copy);
    double error = frobenius(X_copy-X);
    myra::out() << "  |iswap_rows(swap_rows(X))-X| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

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

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

  // Check that swap_columns() and iswap_columns() are reverse operations.
    {
    Matrix<double> X_copy = X;
     swap_columns(swaps,X_copy);
    iswap_columns(swaps,X_copy);
    double error = frobenius(X_copy-X);
    myra::out() << "  |iswap_rows(swap_rows(X))-X| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  // Check that a perm and it's inverse_perm() lead to inverse matrices.
    {
    std::vector<int> iperm = inverse_perm(perm);
    std::vector<int> iswaps = perm2swaps(iperm);
    Matrix<double> iP = swaps2Matrix<double>(iswaps);
    double error = frobenius(gemm(P,iP)-Matrix<double>::identity(I));
    myra::out() << "  |I - P'inverse(P))| = " << error << std::endl;
    REQUIRE(error < 1.0e-9);
    }

  }

} // namespace

ADD_TEST("swaps","[dense]")
  {
  myra::out() << "Testing swaps.." << std::endl;
  test_swaps(20);
  myra::out() << std::endl;
  }

Results: [PASS]

Testing swaps..
  |P*X-swap_rows(X)| = 0
  |P'*X-iswap_rows(X)| = 0
  |iswap_rows(swap_rows(X))-X| = 0
  |X*P-swap_columns(X)| = 0
  |X*P'-iswap_columns(X)| = 0
  |iswap_rows(swap_rows(X))-X| = 0
  |I - P'inverse(P))| = 0

Go back to Summary of /test programs.