zLDLHSolver

Source: tests/dense/ZLDLHSolver.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/herk.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>
#include <myramath/dense/gesvd.h>

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

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

// Reporting.
#include <tests/myratest.h>
#include <myramath/io/FileStreams.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 + hermitian(G);

  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Number> L(A);
  typedef ZLDLHSolver<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.
    {
    // ZLDLHSolver is well suited for computing symmetric schur complements, 
    // S = Bh*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 (hermitian(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<Number> S2(J);
    solver.solveL(B,'L','N');
    herk_inplace(S2,B.top(n_plus),'H',one,one);
    herk_inplace(S2,B.bottom(n_minus),'H',-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.
    {
    // ZLTDHSolver can also do unsymmetric schur complements, S = Ch*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 = (hermitian(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),'H',B.top(n_plus),'N',one,one);
    gemm_inplace(S2,C.bottom(n_minus),'H',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<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) = conjugate(phasor);
  A(2,2) = one;
  // Restrict A to lower triangle, generate LDL' solver for it.
  LowerMatrix<Number> L(A);
  typedef ZLDLHSolver<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("cLDLHSolver","[dense][solver]")
  { test<float>(43,14,1.0e-3f); }

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

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

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

Results: [PASS]

std::complex<double>
  |inv(A)*B-X| = 4.05524e-13
  |inv(A^T)*B-X| = 4.14987e-13
  |inv(A^H)*B-X| = 4.31034e-13
  |inv(A^C)*B-X| = 4.68652e-13
  |B*inv(A)-X| = 3.55262e-13
  |B*inv(A^T)-X| = 4.06983e-13
  |B*inv(A^H)-X| = 4.0139e-13
  |B*inv(A^C)-X| = 4.33901e-13
  |B'*inv(A)*B - (L\B)'*I*(L\B)| = 5.47594e-13
  |C'*inv(A)*B - (L\C)'*I*(L\B)| = 4.93021e-13
  |inv(A)*B-X| (saved) = 4.07922e-13

Go back to Summary of /test programs.