cLDLTSolver

Source: tests/dense/ZLDLTSolver.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/frobenius.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/conjugate.h>
#include <myramath/dense/inverse.h>

// Solver class under test.
#include <myramath/dense/ZLDLTSolver.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 hermitian matrix, A = G+G'
  auto G = Matrix<Number>::random(I,I);
  auto A = G + transpose(G);
  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Number> L(A);
  typedef ZLDLTSolver<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(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);
    }
  // Generate symmetric schur complement.
    {
    // ZLDLTSolver 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<Number>::random(I,J);
    LowerMatrix<Number> S1 (transpose(B)*inverse(A)*B);
    // Structure exploiting answer, approximately half the work and inplace on B.
    Number one(1);
    int n_plus  = solver.inertia().first;
    int n_minus = solver.inertia().second;
    LowerMatrix<Number> 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)'*(L\\B)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    } 
  // Generate unsymmetric schur complement.
    {
    // ZLTDTSolver 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<Number>::random(I,J);
    auto C = Matrix<Number>::random(I,J+3);
    Matrix<Number> S1 = (transpose(C)*inverse(A)*B);
    // Structure exploiting answer.
    Number one(1);
    int n_plus  = solver.inertia().first;
    int n_minus = solver.inertia().second;
    Matrix<Number> 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)'*(L\\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 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| (saved)= " << error << std::endl;  
    REQUIRE(error < tolerance);
        }
  }

template<class Precision> void test3(Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;
  // Construct A = [0 1 0; 1 0 0; 0 0 1], a problematic example.
  auto A = Matrix<Number>::zeros(3,3);
  Precision one(1);
  Precision two(2);
  Precision pi = std::acos(-one);
  Precision phase = two*pi*random<Precision>();
  Number phasor(std::cos(phase),std::sin(phase));
  A(1,0) = phasor;
  A(0,1) = phasor;
  A(2,2) = one;
  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Number> L(A);
  typedef ZLDLTSolver<Precision> Solver;
  Solver solver(L);
  // Solve A*X=B
  auto X = Matrix<Number>::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("cLDLTSolver","[dense][solver]")
  { test<float>(43,14,1.0e-3f); }

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

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

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

Results: [PASS]

std::complex<float>
  |inv(A)*B-X| = 1.07733e-05
  |inv(A^T)*B-X| = 1.40538e-05
  |inv(A^H)*B-X| = 1.10574e-05
  |inv(A^C)*B-X| = 1.15433e-05
  |B*inv(A)-X| = 1.14482e-05
  |B*inv(A^T)-X| = 1.20242e-05
  |B*inv(A^H)-X| = 1.19082e-05
  |B*inv(A^C)-X| = 1.24889e-05
  |B'*inv(A)*B - (L\B)'*(L\B)| = 1.22627e-05
  |C'*inv(A)*B - (L\C)'*(L\B)| = 1.20489e-05
  |inv(A)*B-X| (saved)= 1.15767e-05

Go back to Summary of /test programs.