ICholeskySolver

Source: tests/iterative/ICholeskySolver.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/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/iterative/ICholeskySolver.h>

// Algorithms.
#include <myramath/dense/frobenius.h>
#include <myramath/sparse/phasor.h>
#include <myramath/sparse/laplacian2.h>

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

using namespace myra;

namespace {

// Tests real-valued ICholeskySolver
template<class Precision> void test1(int I, int J, Precision tolerance)
  {
  typedef Precision Number;
  myra::out() << typestring<Number>() << std::endl;
  // Make a laplacian of an IxJ structured grid.
  auto A = laplacian2<Number>(I,J);
  ICholeskySolver<Number> solver(A);
  int N = I*J;
  // Solve an identity from the right and the left, make sure they are equal.
  auto eyeL = Matrix<Number>::identity(N);
  solver.solve(eyeL,'L');
  auto eyeR = Matrix<Number>::identity(N);
  solver.solve(eyeR,'R');
  Precision error = frobenius(eyeL-eyeR);
  myra::out() << "|L'\\L\\I - I/L/L'| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

// Tests complex-valued ICholeskySolver
template<class Precision> void test2(int I, int J, Precision tolerance)
  {
  typedef std::complex<Precision> Number;
  myra::out() << typestring<Number>() << std::endl;  
  // Make a laplacian of an IxJ structured grid.
  auto A = laplacian2<Number>(I,J);
  phasor_hermitian(A);
  ICholeskySolver<Number> solver(A);
  int N = I*J;
  // Solve an identity from the right and the left, make sure they are equal.
  auto eyeL = Matrix<Number>::identity(N);
  solver.solve(eyeL,'L');
  auto eyeR = Matrix<Number>::identity(N);
  solver.solve(eyeR,'R');
  Precision error = frobenius(eyeL-eyeR);
  myra::out() << "|L'\\L\\I - I/L/L'| = " << error << std::endl;
  REQUIRE(error < tolerance);
  }

} // namespace

ADD_TEST("ICholeskySolver","[iterative]")
  {
  // Real-valued.
  test1<float> (20,20,1.0e-4f);
  test1<double>(20,20,1.0e-8);  
  // Complex-valued.
  test2<float> (20,20,1.0e-4f);
  test2<double>(20,20,1.0e-8);
  }
  

Results: [PASS]

float
|L'\L\I - I/L/L'| = 1.88734e-07
double
|L'\L\I - I/L/L'| = 4.0476e-16
std::complex<float>
|L'\L\I - I/L/L'| = 1.99883e-07
std::complex<double>
|L'\L\I - I/L/L'| = 4.25093e-16

Go back to Summary of /test programs.