ppotrf2

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

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

// Parallel algorithms.
#include <myramath/pdense/pherk.h>
#include <myramath/pdense/phemm.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 hermitian positive definite matrix A
  auto G = Matrix<Number>::random(I,I);
  auto A = gemm(G,G,'H');
  LowerMatrix<Number> L(A);
  auto X = Matrix<Number>::random(I,J);
  auto B = phemm('L',L,X);
  Number one(1);
  // Initialize options.
  auto options = Options::create().set_nthreads(4).set_blocksize(128);
  // Measure factor/solve times for serial routines.
  auto L_serial = L;
  auto X_serial = B;
  potrf_inplace(L_serial);
  trsm_inplace('L','N',L_serial,X_serial);
  trsm_inplace('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 = L;
  auto X_parallel = B;
  ppotrf_inplace(L_parallel,options);
  ptrsm_inplace('L','N',L_parallel,X_parallel,'N',one,options);
  ptrsm_inplace('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);
  }

} // namespace

ADD_TEST("ppotrf2","[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.05581e-07
  |L*L'*X-B|, parallel = 2.92997e-07
double
  |L*L'*X-B|, serial   = 4.72538e-16
  |L*L'*X-B|, parallel = 4.25729e-16
std::complex<float>
  |L*L'*X-B|, serial   = 3.13721e-07
  |L*L'*X-B|, parallel = 3.03349e-07
std::complex<double>
  |L*L'*X-B|, serial   = 4.97572e-16
  |L*L'*X-B|, parallel = 4.59068e-16

Go back to Summary of /test programs.