ctrsmU

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

} // namespace

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

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

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

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

Results: [PASS]

std::complex<float>
  |U\B-X| = 1.44508e-05
  |transpose(U)\B-X| = 5.14509e-06
  |hermitian(U)\B-X| = 5.70594e-06
  |conjugate(U)\B-X| = 1.48461e-05
  |B/U-X| = 5.34548e-06
  |B/transpose(U)-X| = 1.49413e-05
  |B/hermitian(U)-X| = 1.40171e-05
  |B/conjugate(U)-X| = 5.26865e-06

Go back to Summary of /test programs.