ctrsmL

Source: tests/dense/trsm1.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/MatrixRange.h>

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/dense/potrf.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/hermitian.h>
#include <myramath/dense/conjugate.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

template<class Number> void test(typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;  
  int I = 27; 
  int J = 34;
  // Compute the cholesky triangle of a random SPD matrix to get a suitably conditioned L.
  auto G = Matrix<Number>::random(I,I);
  Matrix<Number> A = G * hermitian(G);
  potrf_inplace('L','R',A);
  Matrix<Number> L = tril(A);
  char diag = 'N';
  Number alpha = random<Number>();
  Number one(1);
  // Case 1: solving L*X=B
    {
    auto X = Matrix<Number>::random(I,J);  
    Matrix<Number> B = (one/alpha)*L*X;
    trsm_inplace('L','L','N',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |L\\B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 2: solving transpose(L)*X = B
    {
    auto X = Matrix<Number>::random(I,J);    
    Matrix<Number> B = (one/alpha)*transpose(L)*X;
    trsm_inplace('L','L','T',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |transpose(L)\\B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 3: solving hermitian(L)*X = B
    {
    auto X = Matrix<Number>::random(I,J);    
    Matrix<Number> B = (one/alpha)*hermitian(L)*X;
    trsm_inplace('L','L','H',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |hermitian(L)\\B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 4: solving conjugate(L)*X = B
    { 
    auto X = Matrix<Number>::random(I,J);    
    Matrix<Number> B = (one/alpha)*conjugate(L)*X;
    trsm_inplace('L','L','C',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |conjugate(L)\\B-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 5: solving X*L = B
    {    
    auto X = Matrix<Number>::random(J,I);  
    Matrix<Number> B = (one/alpha)*X*L;
    trsm_inplace('R','L','N',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |B/L-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 6: solving X*transpose(L) = B
    {
    auto X = Matrix<Number>::random(J,I);  
    Matrix<Number> B = (one/alpha)*X*transpose(L);
    trsm_inplace('R','L','T',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |B/transpose(L)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 7: solving X*hermitian(L) = B
    {
    auto X = Matrix<Number>::random(J,I);  
    Matrix<Number> B = (one/alpha)*X*hermitian(L);
    trsm_inplace('R','L','H',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |B/hermitian(L)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 8: solving X*conjugate(L) = B
    {
    auto X = Matrix<Number>::random(J,I);  
    Matrix<Number> B = (one/alpha)*X*conjugate(L);
    trsm_inplace('R','L','C',A,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |B/conjugate(L)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

ADD_TEST("strsmL","[dense][lapack]")
  { test<NumberS>(1.0e-2f); }

ADD_TEST("dtrsmL","[dense][lapack]")
  { test<NumberD>(1.0e-10); }

ADD_TEST("ctrsmL","[dense][lapack]")
  { test<NumberC>(1.0e-2f); }

ADD_TEST("ztrsmL","[dense][lapack]")
  { test<NumberZ>(1.0e-10); }

Results: [PASS]

std::complex<float>
  |L\B-X| = 1.23347e-05
  |transpose(L)\B-X| = 5.50482e-06
  |hermitian(L)\B-X| = 6.47073e-06
  |conjugate(L)\B-X| = 1.32608e-05
  |B/L-X| = 6.05496e-06
  |B/transpose(L)-X| = 1.1193e-05
  |B/hermitian(L)-X| = 1.10228e-05
  |B/conjugate(L)-X| = 5.45893e-06

Go back to Summary of /test programs.