SymmAction

Source: tests/iterative/SymmAction.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/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>
#include <myramath/sparse/gemm.h>
#include <myramath/sparse/transpose.h>

// Iterative solve/action stuff.
#include <myramath/iterative/Action.h>
#include <myramath/iterative/SymmAction.h>

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

using namespace myra;

namespace {

// Test make_SymmAction() for a CMatrixRange
template<class Number> void test1(int IJ, int K, typename ReflectPrecision<Number>::type tolerance)
  {  
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Make a random symmetric Matrix A, wrap it up into a SymmAction
  auto A = Matrix<Number>::random(IJ,IJ);
  A = A + transpose(A);
  auto action = make_SymmAction(A);
  // Check action.make_Matrix()
    {
    Precision error = frobenius(action.make_Matrix()-A);
    myra::out() << "  |A(action) - A(dense)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Check A_action.multiply(X,B)
    {
    auto X = Matrix<Number>::random(IJ,K);
    auto B = Matrix<Number>::zeros(IJ,K);
    action.multiply(X,B);
    Precision error = frobenius(A*X-B);
    myra::out() << "  |A*X - B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Check A_action.multiply(X,B,alpha,beta)
    {
    auto X = Matrix<Number>::random(IJ,K);
    auto B = Matrix<Number>::random(IJ,K);
    auto alpha = random<Number>();
    auto beta  = random<Number>();
    auto C = alpha*A*X + beta*B;
    action.multiply(X,B,alpha,beta);
    Precision error = frobenius(B-C);
    myra::out() << "  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

// Test make_SymmAction() for a CSparseMatrixRange
template<class Number> void test2(int IJ, int N, int K, typename ReflectPrecision<Number>::type tolerance)
  {  
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  // Make a random symmetric SparseMatrix A, wrap it up into a SymmAction
  auto A = SparseMatrix<Number>::random(IJ,IJ,N);
  A = A + transpose(A);
  auto action = make_SymmAction(A);
  // Check action.make_Matrix()
    {
    Precision error = frobenius(action.make_Matrix()-A.make_Matrix());
    myra::out() << "  |A(action) - A(dense)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Check A_action.multiply(X,B)
    {
    auto X = Matrix<Number>::random(IJ,K);
    auto B = Matrix<Number>::zeros(IJ,K);
    action.multiply(X,B);
    Precision error = frobenius(A*X-B);
    myra::out() << "  |A*X - B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Check A_action.multiply(X,B,alpha,beta)
    {
    auto X = Matrix<Number>::random(IJ,K);
    auto B = Matrix<Number>::random(IJ,K);
    auto alpha = random<Number>();
    auto beta  = random<Number>();
    auto C = alpha*A*X + beta*B;
    action.multiply(X,B,alpha,beta);
    Precision error = frobenius(B-C);
    myra::out() << "  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

ADD_TEST("SymmAction","[iterative]")
  {
  // Test dense make_SymmAction()
  int IJ = 10;
  int K = 7;
  test1<NumberS>(IJ,K,1.0e-4f);
  test1<NumberD>(IJ,K,1.0e-12);
  test1<NumberC>(IJ,K,1.0e-4f);
  test1<NumberZ>(IJ,K,1.0e-12);
  // Test sparse make_SymmAction()
  int N = 20;
  test2<NumberS>(IJ,N,K,1.0e-4f);
  test2<NumberD>(IJ,N,K,1.0e-12);
  test2<NumberC>(IJ,N,K,1.0e-4f);
  test2<NumberZ>(IJ,N,K,1.0e-12);  
  }

Results: [PASS]

float
  |A(action) - A(dense)| = 0
  |A*X - B| = 6.97664e-07
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 8.91215e-08
double
  |A(action) - A(dense)| = 0
  |A*X - B| = 1.09628e-15
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 1.42849e-15
std::complex<float>
  |A(action) - A(dense)| = 0
  |A*X - B| = 2.53684e-06
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 2.4389e-06
std::complex<double>
  |A(action) - A(dense)| = 0
  |A*X - B| = 4.98497e-15
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[dense] | = 2.66469e-15
float
  |A(action) - A(dense)| = 0
  |A*X - B| = 0
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 1.23241e-07
double
  |A(action) - A(dense)| = 0
  |A*X - B| = 0
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 3.97415e-16
std::complex<float>
  |A(action) - A(dense)| = 0
  |A*X - B| = 0
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 5.69568e-07
std::complex<double>
  |A(action) - A(dense)| = 0
  |A*X - B| = 0
  | (alpha*A*X+beta*B)[action] - (alpha*A*X+beta*B)[sparse] | = 4.30113e-16

Go back to Summary of /test programs.