zCholeskySolver

Source: tests/dense/ZCholeskySolver.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/LowerMatrix.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/herk.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/ZCholeskySolver.h>

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

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

using namespace myra;

namespace {

template<class Precision> void test(int I, int J, Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;
  // Construct random symmetric positive definite matrix, A = G'*G
  Matrix<Number> G = Matrix<Number>::random(I,I);
  Matrix<Number> A = gemm(G,'H',G);
  // Restrict A to lower triangle, generate cholesky solver for it.
  LowerMatrix<Number> L = herk(G,'H');
  typedef ZCholeskySolver<Precision> Solver;
  Solver solver(L);
  // 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(X-B)/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(X-B)/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(X-B)/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(X-B)/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(X-B)/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(X-B)/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(X-B)/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(X-B)/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(X-B)/frobenius(X);
    myra::out() << "  |inv(A)*B-X| (saved) = " << error << std::endl;  
    REQUIRE(error < tolerance);
        }
  }

} // namespace

ADD_TEST("cCholeskySolver","[dense][solver]")
  { test<float>(25,14,1.0e-2f); }

ADD_TEST("zCholeskySolver","[dense][solver]")
  { test<double>(82,45,1.0e-10); }

Results: [PASS]

std::complex<double>
  |inv(A)*B-X| = 5.11768e-13
  |inv(A^T)*B-X| = 5.78896e-13
  |inv(A^H)*B-X| = 5.26593e-13
  |inv(A^C)*B-X| = 5.13542e-13
  |B*inv(A)-X| = 5.32481e-13
  |B*inv(A^T)-X| = 5.09338e-13
  |B*inv(A^H)-X| = 5.25113e-13
  |B*inv(A^C)-X| = 5.55489e-13
  |inv(A)*B-X| (saved) = 5.37226e-13

Go back to Summary of /test programs.