cpotrf

Source: tests/dense/potrf.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/DiagonalMatrix.h>

// Data generators.
#include <myramath/dense/triu.h>
#include <myramath/dense/tril.h>
#include <myramath/dense/diag.h>

// Algorithms.
#include <myramath/dense/potrf.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/frobenius.h>

// Reporting.
#include <myramath/utility/Timer.h>
#include <tests/myratest.h>

using namespace myra;

namespace {

template<class Number> void test(int N, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  
  // Construct random symmetric positive definite A.
  auto G = Matrix<Number>::random(N,N);
  auto A = gemm(G,G,'H') + 10*Matrix<Number>::identity(N);
  
  // Check potrf('L','R'): it factors A = L*L', touching only tril(A)
    {
    auto L = A;
    Timer timer;
    potrf_inplace('L','R',L);
    double t = timer.elapsed_time();
    // Verify that triu(L) is untouched / still matches triu(A)
    Precision error1 = frobenius( (triu(L)-diag(L).make_Matrix()) - (triu(A)-diag(A).make_Matrix()) ) / frobenius(triu(A)-diag(A).make_Matrix());
    myra::out() << "  |triu(L)-triu(A)| = " << error1 << std::endl;
    REQUIRE(error1 < tolerance);
    // Verify that L*L' = A
    L = tril(L);
    Precision error2 = frobenius(gemm(L,L,'H')-A) / frobenius(A);
    myra::out() << "  |L*L'-A| = " << error2 << " (" << t << " sec)" << std::endl;
    REQUIRE(error2 < tolerance);
    }

  // Check potrf('L','L'): it factors A = L'*L, touching only tril(A)
    {
    auto L = A;
    Timer timer;
    potrf_inplace('L','L',L);
    double t = timer.elapsed_time();
    // Verify that triu(L) is untouched / still matches triu(A)
    Precision error1 = frobenius( (triu(L)-diag(L).make_Matrix()) - (triu(A)-diag(A).make_Matrix()) ) / frobenius(triu(A)-diag(A).make_Matrix());
    myra::out() << "  |triu(L)-triu(A)| = " << error1 << std::endl;
    REQUIRE(error1 < tolerance);
    // Verify that L*L' = A
    L = tril(L);
    Precision error2 = frobenius(gemm(L,'H',L)-A) / frobenius(A);
    myra::out() << "  |L'*L-A| = " << error2 << " (" << t << " sec)" << std::endl;
    REQUIRE(error2 < tolerance);
    }

  // Check potrf('U','R'): it factors A = U*U', touching only triu(A)
    {
    auto U = A;
    Timer timer;
    potrf_inplace('U','R',U);
    double t = timer.elapsed_time();
    // Verify that tril(U) is untouched / still matches tril(A)
    Precision error1 = frobenius( (tril(U)-diag(U).make_Matrix()) - (tril(A)-diag(A).make_Matrix()) ) / frobenius(tril(A)-diag(A).make_Matrix());
    myra::out() << "  |tril(U)-tril(A)| = " << error1 << std::endl;
    REQUIRE(error1 < tolerance);
    // Verify that U*U' = A
    U = triu(U);
    Precision error2 = frobenius(gemm(U,U,'H')-A) / frobenius(A);
    myra::out() << "  |U*U'-A| = " << error2 << " (" << t << " sec)" << std::endl;
    REQUIRE(error2 < tolerance);
    }

  // Check potrf('U','L'): it factors A = U'*U, touching only triu(A)
    {
    auto U = A;
    Timer timer;
    potrf_inplace('U','L',U);
    double t = timer.elapsed_time();
    // Verify that tril(U) is untouched / still matches tril(A)
    Precision error1 = frobenius( (tril(U)-diag(U).make_Matrix()) - (tril(A)-diag(A).make_Matrix()) ) / frobenius(tril(A)-diag(A).make_Matrix());
    myra::out() << "  |tril(U)-tril(A)| = " << error1 << std::endl;
    REQUIRE(error1 < tolerance);
    // Verify that U'*U = A
    U = triu(U);
    Precision error2 = frobenius(gemm(U,'H',U)-A) / frobenius(A);
    myra::out() << "  |U'*U-A| = " << error2 << " (" << t << " sec)" << std::endl;
    REQUIRE(error2 < tolerance);
    }
  }

} // namespace

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

ADD_TEST("dpotrf","[dense][lapack]")  
  { test<NumberD>(128,1.0e-9); }

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

ADD_TEST("zpotrf","[dense][lapack]")  
  { test<NumberZ>(128,1.0e-9); }

Results: [PASS]

std::complex<float>
  |triu(L)-triu(A)| = 0
  |L*L'-A| = 1.22153e-07 (0 sec)
  |triu(L)-triu(A)| = 0
  |L'*L-A| = 1.68518e-07 (0 sec)
  |tril(U)-tril(A)| = 0
  |U*U'-A| = 1.68872e-07 (0.001 sec)
  |tril(U)-tril(A)| = 0
  |U'*U-A| = 7.54459e-08 (0 sec)

Go back to Summary of /test programs.