phetrf2

Source: tests/pdense/phetrf2.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>
#include <myramath/dense/intRange.h>

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

// Parallel algorithms.
#include <myramath/pdense/phemm.h>
#include <myramath/pdense/phetrf.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 Precision> void test(int I, int J, Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;
  // Construct random hermitian matrix A
  auto A = Matrix<Number>::random(I,I);  
  A = A + hermitian(A);
  // A.transform([](const NumberZ& a){return NumberZ(a.real(),0);});
  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);
  // Test phetrf('L')
    {
    // Measure factor/solve times for serial routine.
    auto L_serial = L;
    auto X_serial = B;
    auto PQRI_serial = hetrf_inplace(L_serial);
    // This routine factors A = P'*L*Q'*I*Q*L'*P, several steps required for backsolution.
    swap_rows(PQRI_serial.P_swaps,X_serial);
    trsm_inplace('L', 'N', L_serial, X_serial);
    PQRI_serial.R.solve(X_serial,'L','N');
    swap_rows(PQRI_serial.Q_swaps,X_serial);
    X_serial.bottom(PQRI_serial.n_minus) *= -one;
    iswap_rows(PQRI_serial.Q_swaps,X_serial);
    PQRI_serial.R.solve(X_serial,'L','H');
    trsm_inplace('L', 'H', L_serial, X_serial);
    iswap_rows(PQRI_serial.P_swaps,X_serial);
    Precision e_serial = frobenius(A*X_serial-B) / frobenius(B);
    // Measure factor/solve times for serial routine.
    auto L_parallel = L;
    auto X_parallel = B;
    auto PQRI_parallel = phetrf_inplace(L_parallel,options);
    // This routine factors A = P'*L*Q'*I*Q*L'*P, several steps required for backsolution.
    swap_rows(PQRI_parallel.P_swaps,X_parallel);
    ptrsm_inplace('L', 'N', L_parallel, X_parallel, 'N', one, options);
    PQRI_parallel.R.solve(X_parallel,'L','N');
    swap_rows(PQRI_parallel.Q_swaps,X_parallel);
    X_parallel.bottom(PQRI_parallel.n_minus) *= -one;
    iswap_rows(PQRI_parallel.Q_swaps,X_parallel);
    PQRI_parallel.R.solve(X_parallel,'L','H');
    ptrsm_inplace('L', 'H', L_parallel, X_parallel, 'N', one, options);
    iswap_rows(PQRI_parallel.P_swaps,X_parallel);
    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("phetrf2","[pdense][parallel]")
  {
  int I = 512;
  int J = 256;
  test<float> (I,J,1.0e-3f);  
  test<double> (I,J,1.0e-8);
  }

Results: [PASS]

std::complex<float>
  |L*L'*X-B|, serial   = 1.8843e-05
  |L*L'*X-B|, parallel = 2.07508e-05
std::complex<double>
  |L*L'*X-B|, serial   = 3.37864e-14
  |L*L'*X-B|, parallel = 3.73965e-14

Go back to Summary of /test programs.