zLUSolver

Source: tests/dense/LUSolver.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/intRange.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/conjugate.h>

// Solver class under test.
#include <myramath/dense/LUSolver.h>

// Serialization testing.
#include <myramath/io/VectorStream.h>

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Construct random matrix and a solver for it.
  const auto A = Matrix<Number>::random(I,I);
  typedef LUSolver<Number> Solver;
  Solver solver(A);
  // Total of 8 solve() variations to validiate, side['L',R'] x op['N','T','H','C']
  // Solve A*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    auto B = A*X;
    solver.solve(B,'L','N');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |inv(A)*B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve transpose(A)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    auto B = transpose(A)*X;
    solver.solve(B,'L','T');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |inv(A^T)*B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve hermitian(A)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    auto B = hermitian(A)*X;
    solver.solve(B,'L','H');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |inv(A^H)*B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve conjugate(A)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    auto B = conjugate(A)*X;
    solver.solve(B,'L','C');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |inv(A^C)*B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve X*A=B
    {
    auto X = Matrix<Number>::random(J,I);
    auto B = X*A;
    solver.solve(B,'R','N');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |B*inv(A)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve X*transpose(A)=B
    {
    auto X = Matrix<Number>::random(J,I);
    auto B = X*transpose(A);
    solver.solve(B,'R','T');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |B*inv(A^T)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve X*hermitian(A)=B
    {
    auto X = Matrix<Number>::random(J,I);
    auto B = X*hermitian(A);
    solver.solve(B,'R','H');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |B*inv(A^H)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Solve X*conjugate(A)=B
    {
    auto X = Matrix<Number>::random(J,I);
    auto B = X*conjugate(A);
    solver.solve(B,'R','C');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |B*inv(A^C)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Save solver to disk, then solve A*X=B again.
    {
    VectorStream inout;
    solver.write(inout);
    Solver saved(inout);
    auto X = Matrix<Number>::random(I,J);
    auto B = A*X;
    saved.solve(B,'L','N');
    Precision error = frobenius(B-X)/frobenius(X);
    myra::out() << "  |inv(A)*B-X| = " << error << std::endl;  
    REQUIRE(error < tolerance);
        }
  }

} // namespace

ADD_TEST("sLUSolver","[dense][solver]")
  { test<NumberS>(200,120,1.0e-2f); }

ADD_TEST("dLUSolver","[dense][solver]")
  { test<NumberD>(200,120,1.0e-10); }

ADD_TEST("cLUSolver","[dense][solver]")
  { test<NumberC>(200,120,1.0e-2f); }

ADD_TEST("zLUSolver","[dense][solver]")
  { test<NumberZ>(200,120,1.0e-10); }

Results: [PASS]

std::complex<double>
  |inv(A)*B-X| = 6.50334e-14
  |inv(A^T)*B-X| = 6.92375e-14
  |inv(A^H)*B-X| = 6.28447e-14
  |inv(A^C)*B-X| = 6.13314e-14
  |B*inv(A)-X| = 6.36775e-14
  |B*inv(A^T)-X| = 6.53646e-14
  |B*inv(A^H)-X| = 6.61393e-14
  |B*inv(A^C)-X| = 6.7714e-14
  |inv(A)*B-X| = 6.65871e-14

Go back to Summary of /test programs.