dCholeskySolver

Source: tests/dense/RCholeskySolver.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>
#include <myramath/dense/MatrixRange.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/syrk.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/frobenius.h>

// Solver class under test.
#include <myramath/dense/RCholeskySolver.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)
  {  
  myra::out() << typestring<Precision>() << std::endl;
  // Construct random symmetric positive definite matrix, A = G'*G
  auto G = Matrix<Precision>::random(I,I);
  Matrix<Precision> A = gemm(G,'T',G) + 10*Matrix<Precision>::identity(I);
  // Restrict A to lower triangle, generate cholesky solver for it.
  LowerMatrix<Precision> L(A);
  typedef RCholeskySolver<Precision> Solver;
  Solver solver(L);
  // Solve A*X=B
    {
    auto B = Matrix<Precision>::random(I,J);
    auto X = B;
    solver.solve(X,'L','N');
    Precision error = frobenius(A*X-B)/frobenius(B);
    myra::out() << "  |A*X-B|/|B| = " << error << std::endl;  
    REQUIRE(error < tolerance);
    }
  // Solve X*A=B
    {
    auto B = Matrix<Precision>::random(J,I);
    auto X = B;
    solver.solve(X,'R','N');
    Precision error = frobenius(X*A-B)/frobenius(B);
    myra::out() << "  |X*A-B|/|B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Roundtrip solver to disk, then solve A*X=B again.
    {
    VectorStream inout;
    solver.write(inout);
    Solver saved(inout);
    auto B = Matrix<Precision>::random(I,J);
    auto X = B;
    saved.solve(X,'L','N');
    Precision error = frobenius(A*X-B)/frobenius(B);
    myra::out() << "  |inv(A)*B-X| (saved) = " << error << std::endl;  
    REQUIRE(error < tolerance);
        }
  };

} // namespace

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

ADD_TEST("dCholeskySolver","[dense][solver]")
  { test<double>(25,14,1.0e-10); }

Results: [PASS]

double
  |A*X-B|/|B| = 4.00191e-16
  |X*A-B|/|B| = 3.44249e-16
  |inv(A)*B-X| (saved) = 3.43516e-16

Go back to Summary of /test programs.