ctrmmU

Source: tests/dense/trmm2.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/triu.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> U = triu(A);
  char diag = 'N';
  Number alpha = random<Number>();
  // Case 1: multiplying U*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*U*X;
    trmm_inplace('L','U','N',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |U*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 2: multiplying transpose(U)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*transpose(U)*X;
    trmm_inplace('L','U','T',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |transpose(U)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 3: multiplying hermitian(U)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*hermitian(U)*X;
    trmm_inplace('L','U','H',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |hermitian(U)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 4: multiplying conjugate(U)*X = B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*conjugate(U)*X;
    trmm_inplace('L','U','C',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |conjugate(U)*X-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 5: multiplying X*U = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*U;
    trmm_inplace('R','U','N',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*U-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 6: multiplying X*transpose(U) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*transpose(U);
    trmm_inplace('R','U','T',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*tranpose(U)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 7: multiplying X*hermitian(U) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*hermitian(U);
    trmm_inplace('R','U','H',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*hermitian(U)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 8: multiplying X*conjugate(U) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*conjugate(U);
    trmm_inplace('R','U','C',A,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*conjugate(U)-B| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

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

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

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

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

Results: [PASS]

std::complex<float>
  |U*X-B| = 1.0123e-05
  |transpose(U)*X-B| = 8.99374e-06
  |hermitian(U)*X-B| = 8.93056e-06
  |conjugate(U)*X-B| = 1.05462e-05
  |X*U-B| = 9.87578e-06
  |X*tranpose(U)-B| = 9.4355e-06
  |X*hermitian(U)-B| = 9.79563e-06
  |X*conjugate(U)-B| = 9.64413e-06

Go back to Summary of /test programs.