ctrsmLower

Source: tests/dense/trsm3.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>
#include <myramath/dense/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.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);
  A = tril(A);
  LowerMatrix<Number> L(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)*A*X;
    trsm_inplace('L','N',L,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(A)*X;
    trsm_inplace('L','T',L,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(A)*X;
    trsm_inplace('L','H',L,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(A)*X;
    trsm_inplace('L','C',L,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*A;
    trsm_inplace('R','N',L,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(A);
    trsm_inplace('R','T',L,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(A);
    trsm_inplace('R','H',L,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(A);
    trsm_inplace('R','C',L,B,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |B/conjugate(L)-X| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

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

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

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

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

Results: [PASS]

std::complex<float>
  |L\B-X| = 2.5539e-05
  |transpose(L)\B-X| = 5.80269e-06
  |hermitian(L)\B-X| = 5.66601e-06
  |conjugate(L)\B-X| = 2.53304e-05
  |B/L-X| = 6.02365e-06
  |B/transpose(L)-X| = 2.33278e-05
  |B/hermitian(L)-X| = 2.41674e-05
  |B/conjugate(L)-X| = 5.51361e-06

Go back to Summary of /test programs.