sLDLTSolver3

Source: tests/dense/RLDLTSolver.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/LowerMatrixRange.h>

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

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

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

// Reporting.
#include <myramath/utility/stlprint.h>
#include <tests/myratest.h>

using namespace myra;
using namespace myra_stlprint;

namespace {

template<class Precision> void test(int I, int J, Precision tolerance)
  {  
  myra::out() << typestring<Precision>() << std::endl;
  // Construct random symmetric matrix, A = G+G'
  auto G = Matrix<Precision>::random(I,I);
  auto A = G + transpose(G);
  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Precision> L(A);
  typedef RLDLTSolver<Precision> Solver;
  Solver solver(L);
  myra::out() << "inertia = " << solver.inertia() << std::endl;
  // Solve A*X=B
    {
    auto X = Matrix<Precision>::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 X*A=B
    {
    auto X = Matrix<Precision>::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);
    }
  // Generate symmetric schur complement.
    {
    // RLDLTSolver is well suited for computing symmetric schur complements, 
    // S = Bt*inv(A)*B. Try it out and compare to a reference calculation.
    // (In the case of sparse A/B, this approach is dramatically more efficient)
    // Reference answer.
    auto B = Matrix<Precision>::random(I,J);
    LowerMatrix<Precision> S1(transpose(B)*inverse(A)*B);
    // Structure exploiting answer, approximately half the work and inplace on B.
    Precision one(1);
    int n_plus  = solver.inertia().first;
    int n_minus = solver.inertia().second;
    LowerMatrix<Precision> S2(J);
    solver.solveL(B,'L','N');
    syrk_inplace(S2, B.top(n_plus), 'T', one, one);
    syrk_inplace(S2, B.bottom(n_minus), 'T',-one, one);
    // Compare/report.
    Precision error = frobenius(S1-S2) / frobenius(S1);
    myra::out() << "  |B'*inv(A)*B - (L\\B)'*I*(L\\B)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Generate unsymmetric schur complement.
    {
    // RLTDTSolver can also do unsymmetric schur complements, S = Ct*inv(A)*B.
    // Not especially efficient, but still tricky so it's demonstrated here.
    // (In the case of sparse A/B/C, this approach is dramatically more efficient)
    // Reference answer.
    auto B = Matrix<Precision>::random(I,J);
    auto C = Matrix<Precision>::random(I,J+3);
    Matrix<Precision> S1 = (transpose(C)*inverse(A)*B);
    // Structure exploiting answer.
    Precision one(1);
    int n_plus  = solver.inertia().first;
    int n_minus = solver.inertia().second;
    Matrix<Precision> S2(J+3,J);
    solver.solveL(B,'L','N');
    solver.solveL(C,'L','N');
    gemm_inplace(S2, C.top(n_plus) , 'T', B.top(n_plus) , 'N',  one, one);
    gemm_inplace(S2, C.bottom(n_minus), 'T', B.bottom(n_minus), 'N', -one, one);
    // Compare/report.
    Precision error = frobenius(S1-S2) / frobenius(S1);
    myra::out() << "  |C'*inv(A)*B - (L\\C)'*I*(L\\B)| = " << 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<Precision>::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| (saved) = " << error << std::endl;  
    REQUIRE(error < tolerance);
        }
  };

template<class Precision> void test3(Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Construct A = [0 1 0; 1 0 0; 0 0 1], a problematic example.
  auto A = Matrix<Precision>::zeros(3,3);
  Precision a10 = random<Precision>();
  A(1,0) = a10;
  A(0,1) = a10;
  A(2,2) = 1;
  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Precision> L(A);
  typedef RLDLTSolver<Precision> Solver;
  Solver solver(L);
  // Solve A*X=B
  auto X = Matrix<Precision>::random(3,2);
  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);
  }

} // namespace

ADD_TEST("sLDLTSolver","[dense][solver]")
  { test<float>(43,14,1.0e-3f); }

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

ADD_TEST("sLDLTSolver3","[dense][solver]")
  { test3<float>(1.0e-3f); }

ADD_TEST("dLDLTSolver3","[dense][solver]")
  { test3<double>(1.0e-10); }

Results: [PASS]

float
  |inv(A)*B-X| = 5.5718e-08

Go back to Summary of /test programs.