zldlh2_solver

Source: tests/multifrontal/zldlh/zldlh2_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>
#include <myramath/sparse/phasor.h>
#include <myramath/sparse/transpose.h>

// Solver type under test.
#include <myramath/multifrontal/zldlh/SparseZLDLHSolver.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;
  typedef std::complex<Precision> Number;
  int N = I*J;
  // Make laplacian A, shift it indefinite, phase shift to make it nontrivially hermitian.
  SparseMatrix<Number> A = laplacian2<Number>(I,J);  
  Number shift = 5.5;
  for (int n = 0; n < N; ++n)
    A(n,n) -= shift;
  phasor_hermitian(A);
  // Factor A = L*L'
  typedef SparseZLDLHSolver<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 SparseZLDLH.refine() automatically.
  // Test A*x = b, no refinement.
    {
    auto b = Vector<Number>::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<Number>::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 transpose(A)*x = b, using refine()
    {
    auto b = Vector<Number>::random(N);
    auto x = b;
    auto history = solver.refine(x.column(),'L','T');
    Precision residual = frobenius(transpose(A)*x-b) / frobenius(b);
    myra::out() << "  |A^T*x-b| (refine) = " << residual << std::endl;
    myra::out() << "  history = " << history << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test x*A = b, using refine()
    {
    auto b = Vector<Number>::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);
    }
  // Test x*tranpose(A) = b, using refine()
    {
    auto b = Vector<Number>::random(N);
    auto x = b;
    auto history = solver.refine(x.row(),'R','T');
    Precision residual = frobenius(x*transpose(A)-b) / frobenius(b);
    myra::out() << "  |x*A^T-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<Number>::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<Number>::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("zldlh2_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.000508248
  |A*x-b| (refine) = 8.13978e-06
  history = [ 0.000293949 6.34765e-06 ] (2)
  |A^T*x-b| (refine) = 7.67395e-06
  history = [ 0.000383531 5.95696e-06 ] (2)
  |x*A-b| (refine) = 9.24975e-06
  |x*A^T-b| (refine) = 6.34025e-06
  |inv(A)*b-x| (copy) = 8.43896e-06
  |A*x-b| (saved) = 7.02306e-06
double
inertia = {1531,2565}
  |A*x-b| (solve) = 1.03828e-12
  |A*x-b| (refine) = 1.82741e-14
  history = [ 1.12802e-12 1.40342e-14 ] (2)
  |A^T*x-b| (refine) = 2.68098e-14
  history = [ 1.92368e-12 2.13485e-14 ] (2)
  |x*A-b| (refine) = 2.66673e-14
  |x*A^T-b| (refine) = 2.94271e-14
  |inv(A)*b-x| (copy) = 2.10389e-14
  |A*x-b| (saved) = 2.44362e-14

Go back to Summary of /test programs.