dtrmmLower

Source: tests/dense/trmm3.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/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)
  {
  // Useful typedef.
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Take the lower triangle L of a random square matrix A.
  int I = 35; 
  int J = 44;
  Matrix<Number> A = tril(Matrix<Number>::random(I,I));
  LowerMatrix<Number> L(A);
  Number alpha = random<Number>();
  char diag = 'N';
  // Case 1: multiply L*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*A*X;
    trmm_inplace('L','N',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |L*X (gemm) - L*X (trmm)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 2: multiply transpose(L)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*transpose(A)*X;
    trmm_inplace('L','T',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |transpose(L)*X (gemm) - transpose(L)*X (trmm)| = " << error << std::endl;  
    REQUIRE(error < tolerance);
    }
  // Case 3: multiply hermitian(L)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*hermitian(A)*X;
    trmm_inplace('L','H',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |hermitian(L)*X (gemm) - hermitian(L)*X (trmm)| = " << error << std::endl;  
    REQUIRE(error < tolerance);
    }
  // Case 4: multiply conjugate(L)*X=B
    {
    auto X = Matrix<Number>::random(I,J);
    Matrix<Number> B = alpha*conjugate(A)*X;
    trmm_inplace('L','C',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |conjugate(L)*X (gemm) - conjugate(L)*X (trmm)| = " << error << std::endl;  
    REQUIRE(error < tolerance);
    }
  // Case 5: multiply X*L = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*A;
    trmm_inplace('R','N',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*L (gemm) - X*L (trmm)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 6: multiply X*transpose(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*transpose(A);
    trmm_inplace('R','T',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*transpose(L) (gemm) - X*tranpose(L) (trmm)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 7: multiply X*hermitian(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*hermitian(A);
    trmm_inplace('R','H',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*hermitian(L) (gemm) - X*hermitian(L) (trmm)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Case 8: multiply X*conjugate(L) = B
    {
    auto X = Matrix<Number>::random(J,I);
    Matrix<Number> B = alpha*X*conjugate(A);
    trmm_inplace('R','C',L,X,diag,alpha);
    Precision error = frobenius(B-X);
    myra::out() << "  |X*conjugate(L) (gemm) - X*conjugate(L) (trmm)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }

} // namespace

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

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

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

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

Results: [PASS]

double
  |L*X (gemm) - L*X (trmm)| = 9.16888e-15
  |transpose(L)*X (gemm) - transpose(L)*X (trmm)| = 9.84126e-15
  |hermitian(L)*X (gemm) - hermitian(L)*X (trmm)| = 9.94787e-15
  |conjugate(L)*X (gemm) - conjugate(L)*X (trmm)| = 9.67413e-15
  |X*L (gemm) - X*L (trmm)| = 7.94848e-15
  |X*transpose(L) (gemm) - X*tranpose(L) (trmm)| = 8.23117e-15
  |X*hermitian(L) (gemm) - X*hermitian(L) (trmm)| = 7.67924e-15
  |X*conjugate(L) (gemm) - X*conjugate(L) (trmm)| = 9.38727e-15

Go back to Summary of /test programs.