psytrf2

Source: tests/pdense/psytrf2.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/sytrf.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/swaps.h>

// Parallel routines under test.
#include <myramath/pdense/pgemm.h>
#include <myramath/pdense/psytrf.h>
#include <myramath/pdense/ptrsm.h>
#include <myramath/pdense/psymm.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 matrix A
  auto A = Matrix<Number>::random(I,I);
  A = A+transpose(A);
  LowerMatrix<Number> L(A);
  auto X = Matrix<Number>::random(I,J);
  auto B = psymm('L',L,X);
  Number one(1);
  // Initialize options.
  auto options = Options::create().set_nthreads(4).set_blocksize(128);
  // Test psytrf('L')
    {
    // Measure factor/solve times for serial routine.
    auto L_serial = L;
    auto X_serial = B;
    auto PQRI_serial = sytrf_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','T');
    trsm_inplace('L', 'T', 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 = psytrf_inplace(L_parallel);
    // 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','T');
    ptrsm_inplace('L', 'T', 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("psytrf2","[pdense][parallel]")
  {
  int I = 512;
  int J = 256;
  test<NumberS> (I,J,1.0e-4f);
  test<NumberD> (I,J,1.0e-10);
  test<NumberC> (I,J,1.0e-4f);
  test<NumberZ> (I,J,1.0e-10);
  }

Results: [PASS]

float
  |L*L'*X-B|, serial   = 1.45642e-05
  |L*L'*X-B|, parallel = 1.40478e-05
double
  |L*L'*X-B|, serial   = 3.43104e-14
  |L*L'*X-B|, parallel = 2.70661e-14
std::complex<float>
  |L*L'*X-B|, serial   = 2.022e-05
  |L*L'*X-B|, parallel = 1.87708e-05
std::complex<double>
  |L*L'*X-B|, serial   = 3.37046e-14
  |L*L'*X-B|, parallel = 2.93221e-14

Go back to Summary of /test programs.