zldlt2_solver

Source: tests/multifrontal/zldlt/zldlt2_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/hermitian.h>

// Solver type under test.
#include <myramath/multifrontal/zldlt/SparseZLDLTSolver.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 complex-symmetric.
  SparseMatrix<Number> A = laplacian2<Number>(I,J);  
  Number shift = 5.5;
  for (int n = 0; n < N; ++n)
    A(n,n) -= shift;
  phasor_symmetric(A);
  // Factor A = L*L'
  typedef SparseZLDLTSolver<Precision> Solver;
  Solver solver(A);
  // 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 SparseZLDLT.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 hermitian(A)*x = b, using refine()
    {
    auto b = Vector<Number>::random(N);
    auto x = b;
    auto history = solver.refine(x.column(),'L','H');
    Precision residual = frobenius(hermitian(A)*x-b) / frobenius(b);
    myra::out() << "  |A^H*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*hermitian(A) = b, using refine()
    {
    auto b = Vector<Number>::random(N);
    auto x = b;
    auto history = solver.refine(x.row(),'R','H');
    Precision residual = frobenius(x*hermitian(A)-b) / frobenius(b);
    myra::out() << "  |x*A^H-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("zldlt2_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
  |A*x-b| (solve) = 0.000470822
  |A*x-b| (refine) = 1.12287e-05
  history = [ 0.000818973 1.12633e-05 ] (2)
  |A^H*x-b| (refine) = 8.28294e-06
  history = [ 0.000873155 8.3314e-06 ] (2)
  |x*A-b| (refine) = 8.83332e-06
  |x*A^H-b| (refine) = 4.43869e-06
  |inv(A)*b-x| (copy) = 1.25217e-05
  |A*x-b| (saved) = 9.45063e-06
double
  |A*x-b| (solve) = 1.59198e-12
  |A*x-b| (refine) = 1.40318e-14
  history = [ 3.83509e-12 1.40993e-14 ] (2)
  |A^H*x-b| (refine) = 1.74107e-14
  history = [ 2.37101e-12 1.75394e-14 ] (2)
  |x*A-b| (refine) = 8.6354e-15
  |x*A^H-b| (refine) = 1.92932e-14
  |inv(A)*b-x| (copy) = 1.91602e-14
  |A*x-b| (saved) = 1.33625e-14

Go back to Summary of /test programs.