djacobi

Source: tests/dense/jacobi.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>
#include <myramath/dense/DiagonalMatrixRange.h>

// Algorithms.
#include <myramath/dense/dimm.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/jacobi.h>
#include <myramath/dense/transpose.h>
#include <myramath/dense/frobenius.h>

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

using namespace myra;

namespace {

template<class Precision> void test(int N, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  // Fill symmetric A.
  auto A = Matrix<Precision>::random(N,N);
  A = A+transpose(A);
  // Compute (X,D) using jacobi()
  const auto XD = jacobi(A);
  const auto& X = XD.first;
  const auto& D = XD.second;
  // Check errors (orthogonality, residual, convergence to diagonal form).
  Precision I_error = frobenius(gemm(X,'T',X)-Matrix<Precision>::identity(N));
  myra::out() << "  |X'X-I| = " << I_error << std::endl;  
  Precision R_error = frobenius(A*X-X*D);
  myra::out() << "  |A*X-X*D| = " << R_error << std::endl;
  REQUIRE(I_error < tolerance);
  REQUIRE(R_error < tolerance);
  }

} // namespace

ADD_TEST("sjacobi","[dense]")
  { test<float> (10, 1.0e-3f); }

ADD_TEST("djacobi","[dense]")
  { test<double> (20, 1.0e-8); }

Results: [PASS]

double
  |X'X-I| = 1.0386e-14
  |A*X-X*D| = 4.10424e-14

Go back to Summary of /test programs.