rldlt2_solver

Source: tests/multifrontal/rldlt/rldlt2_solver.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/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/sparse/Permutation.h>

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

// Solver under test.
#include <myramath/multifrontal/rldlt/SparseRLDLTSolver.h>

// For testing serialization roundtrip.
#include <myramath/io/VectorStream.h>

// Reporting.
#include <myramath/utility/stlprint.h>
#include <tests/myratest.h>

using namespace myra;
using namespace myra_stlprint;

namespace {

template<class Precision> void test(int I, int J, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  int N = I*J;
  // Make laplacian matrix A, shift it to make it indefinite.
  SparseMatrix<Precision> A = laplacian2<Precision>(I,J);  
  Precision shift = 5.5;
  for (int n = 0; n < N; ++n)
    A(n,n) -= shift;
  // Factor A = L*L'
  typedef SparseRLDLTSolver<Precision> Solver;
  Solver solver(A);
  myra::out() << "inertia = " << solver.inertia() << std::endl;
  // Note that this indefinite solver is not very accurate in
  // single precision, so we only check/throw residual errors 
  // when we're using backwards refinement with them. Note this 
  // refinement is baked into SparseRLDLT.refine() automatically.
  // Test A*x = b, no refinement.
    {
    auto b = Vector<Precision>::random(N);
    auto x = b;
    solver.solve(x.column(),'L','N');
    Precision residual = frobenius(A*x-b) / frobenius(b);
    myra::out() << "  |A*x-b| (solve) = " << residual << std::endl;
    // Note, we don't check the residual on this one.
    // REQUIRE(residual < tolerance);
    }
  // Test A*x = b, using refine()
    {
    auto b = Vector<Precision>::random(N);
    auto x = b;
    auto history = solver.refine(x.column(),'L','N');
    Precision residual = frobenius(A*x-b) / frobenius(b);
    myra::out() << "  |A*x-b| (refine) = " << residual << std::endl;
    myra::out() << "  history = " << history << std::endl;
    // Note, we do check the residual on this one (and all the rest)
    REQUIRE(residual < tolerance);
    }
  // Test x*A = b, using refine()
    {
    auto b = Vector<Precision>::random(N);
    auto x = b;
    auto history = solver.refine(x.row(),'R','N');
    Precision residual = frobenius(x*A-b) / frobenius(b);
    myra::out() << "  |x*A-b| (refine) = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Make a deep copy of solver, test A*x = b again.
    {
    Solver copy = solver;
    auto b = Vector<Precision>::random(N);
    auto x = b;
    auto history = solver.refine(x.column(),'L','N');   
    Precision residual = frobenius(A*x-b) / frobenius(b);
    myra::out() << "  |inv(A)*b-x| (copy) = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Roundtrip solver to file, test A*x = b again.
    {
    VectorStream inout;
    solver.write(inout);
    Solver saved(inout);
    auto b = Vector<Precision>::random(N);
    auto x = b;
    auto history = saved.refine(x.column(),'L','N');
    Precision residual = frobenius(A*x-b) / frobenius(b);
    myra::out() << "  |A*x-b| (saved) = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  }

} // namespace

ADD_TEST("rldlt2_solver","[multifrontal][parallel]")
  {
  // Define problem size.
  int I = 64;
  int J = 64;
  int N = I*J;
  // Run tests.
  test<float >(I,J,1.0e-4f);
  test<double>(I,J,1.0e-8);
  }

Results: [PASS]

float
inertia = {1531,2565}
  |A*x-b| (solve) = 0.00093264
  |A*x-b| (refine) = 1.24751e-05
  history = [ 0.00300359 1.25741e-05 ] (2)
  |x*A-b| (refine) = 5.98628e-06
  |inv(A)*b-x| (copy) = 3.85729e-06
  |A*x-b| (saved) = 8.37572e-06
double
inertia = {1531,2565}
  |A*x-b| (solve) = 4.18021e-13
  |A*x-b| (refine) = 4.91717e-15
  history = [ 2.30188e-13 4.88419e-15 ] (2)
  |x*A-b| (refine) = 1.23369e-14
  |inv(A)*b-x| (copy) = 1.14628e-14
  |A*x-b| (saved) = 1.33477e-14

Go back to Summary of /test programs.