ppotrf1

Source: tests/pdense/ppotrf1.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>

// Serial algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/potrf.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/frobenius.h>

// Parallel algorithms.
#include <myramath/pdense/pgemm.h>
#include <myramath/pdense/ppotrf.h>
#include <myramath/pdense/ptrsm.h>
#include <myramath/pdense/Options.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;
typedef pdense::Options Options;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Construct random symmetric positive definite matrix A, random X, and image B=A*X
  auto G = Matrix<Number>::random(I,I);
  auto A = pgemm(G,G,'H');
  for (int i = 0; i < I; ++i)
    A(i,i) += 1;
  auto X = Matrix<Number>::random(I,J);  
  auto B = pgemm(A,X);
  Number one(1);
  // Initialize options.
  auto options = Options::create().set_nthreads(4).set_blocksize(128);
  // Testing/profiling ppotrf('L')
    {   
    // Measure factor/solve times for serial routines.
    auto L_serial = A;
    auto X_serial = B;
    potrf_inplace('L','R',L_serial);
    trsm_inplace('L','L','N',L_serial,X_serial);
    trsm_inplace('L','L','H',L_serial,X_serial); 
    Precision e_serial = frobenius(A*X_serial-B) / frobenius(B);
    // Measure factor/solve times for parallel routines.    
    auto L_parallel = A;
    auto X_parallel = B;
    ppotrf_inplace('L',L_parallel,options);
    ptrsm_inplace('L','L','N',L_parallel,X_parallel,'N',one,options);
    ptrsm_inplace('L','L','H',L_parallel,X_parallel,'N',one,options);
    Precision e_parallel = frobenius(A*X_parallel-B) / frobenius(B);
    // Report accuracy and time.
    myra::out() << "  |L*L'*X-B|, serial   = " << e_serial << std::endl;
    myra::out() << "  |L*L'*X-B|, parallel = " << e_parallel << std::endl;
    REQUIRE(e_serial < tolerance);
    REQUIRE(e_parallel < tolerance);
    }
  // Testing/profiling ppotrf('U')
    {   
    // Measure factor/solve times for serial routines.
    auto U_serial = A;
    auto X_serial = B;
    potrf_inplace('U','R',U_serial);
    trsm_inplace('L','U','N',U_serial,X_serial);
    trsm_inplace('L','U','H',U_serial,X_serial);
    Precision e_serial = frobenius(A*X_serial-B) / frobenius(B);
    // Measure factor/solve times for parallel routines.    
    auto U_parallel = A;
    auto X_parallel = B;
    ppotrf_inplace('U',U_parallel,options);
    ptrsm_inplace('L','U','N',U_parallel,X_parallel,'N',one,options);
    ptrsm_inplace('L','U','H',U_parallel,X_parallel,'N',one,options);
    Precision e_parallel = frobenius(A*X_parallel-B) / frobenius(B);
    // Report accuracy and time.
    myra::out() << "  |U*U'*X-B|, serial   = " << e_serial << std::endl;
    myra::out() << "  |U*U'*X-B|, parallel = " << e_parallel << std::endl;
    REQUIRE(e_serial < tolerance);
    REQUIRE(e_parallel < tolerance);
    }
  }

} // namespace

ADD_TEST("ppotrf1","[pdense][parallel]")
  {
  int I = 512;
  int J = 256;
  test<NumberS> (I,J,1.0e-4f);
  test<NumberD> (I,J,1.0e-8);
  test<NumberC> (I,J,1.0e-4f);
  test<NumberZ> (I,J,1.0e-8);
  }

Results: [PASS]

float
  |L*L'*X-B|, serial   = 3.8177e-07
  |L*L'*X-B|, parallel = 2.90744e-07
  |U*U'*X-B|, serial   = 3.62301e-07
  |U*U'*X-B|, parallel = 2.86831e-07
double
  |L*L'*X-B|, serial   = 6.62325e-16
  |L*L'*X-B|, parallel = 4.28738e-16
  |U*U'*X-B|, serial   = 6.24621e-16
  |U*U'*X-B|, parallel = 4.26244e-16
std::complex<float>
  |L*L'*X-B|, serial   = 3.82463e-07
  |L*L'*X-B|, parallel = 3.02693e-07
  |U*U'*X-B|, serial   = 3.61825e-07
  |U*U'*X-B|, parallel = 2.9542e-07
std::complex<double>
  |L*L'*X-B|, serial   = 6.65187e-16
  |L*L'*X-B|, parallel = 4.56632e-16
  |U*U'*X-B|, serial   = 6.26351e-16
  |U*U'*X-B|, parallel = 4.52292e-16

Go back to Summary of /test programs.