strmmL

Source: tests/dense/trmm1.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/tril.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/trmm.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;
  // Take the lower triangle of a random square Matrix.
  auto A = Matrix<Number>::random(I,I);
  Matrix<Number> L = tril(A);
  char diag = 'N';
  Number alpha = random<Number>();
  // Case 1: multiplying L*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*L*X;
    trmm_inplace('L','L','N',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |L*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 2: multiplying transpose(L)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*transpose(L)*X;
    trmm_inplace('L','L','T',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |transpose(L)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 3: multiplying hermitian(L)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*hermitian(L)*X;
    trmm_inplace('L','L','H',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |hermitian(L)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 4: multiplying conjugate(L)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*conjugate(L)*X;
    trmm_inplace('L','L','C',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |conjugate(L)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 5: multiplying X*L = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*L;
    trmm_inplace('R','L','N',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*L-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 6: multiplying X*transpose(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*transpose(L);
    trmm_inplace('R','L','T',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*tranpose(L)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 7: multiplying X*hermitian(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*hermitian(L);
    trmm_inplace('R','L','H',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*hermitian(L)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 8: multiplying X*conjugate(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*conjugate(L);
    trmm_inplace('R','L','C',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*conjugate(L)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

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

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

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

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

Results: [PASS]

float
  |L*X-B| = 9.81828e-07
  |transpose(L)*X-B| = 8.12303e-07
  |hermitian(L)*X-B| = 8.21886e-07
  |conjugate(L)*X-B| = 8.90141e-07
  |X*L-B| = 0
  |X*tranpose(L)-B| = 9.61046e-07
  |X*hermitian(L)-B| = 8.80624e-07
  |X*conjugate(L)-B| = 0

Go back to Summary of /test programs.