lu2_solver

Source: tests/multifrontal/lu/lu2_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/Pattern.h>
#include <myramath/sparse/PatternRange.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/transpose.h>
#include <myramath/sparse/hermitian.h>
#include <myramath/sparse/conjugate.h>
#include <myramath/sparse/laplacian2.h>

// Solver under test.
#include <myramath/multifrontal/lu/SparseLUSolver.h>

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

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  myra::out() << typestring<Number>() << std::endl;
  typedef typename ReflectPrecision<Number>::type Precision;
  int N = I*J;
  // Make symmetric pattern, stuff with random values.
  Pattern P = stencil2(I,J);
  auto A = SparseMatrix<Number>::random(P);
  // Factor A = L*U
  typedef SparseLUSolver<Number> Solver;
  Solver solver(A);
  // Test A*x = b
    {
    auto b = Matrix<Number>::random(N,3);
    auto x = b;
    solver.refine(x,'L','N');
    Precision residual = frobenius(A*x-b) / frobenius(b);
    myra::out() << "  |A*b-x| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test transpose(A)*x = b
    {
    auto b = Matrix<Number>::random(N,3);
    auto x = b;
    solver.refine(x,'L','T');
    Precision residual = frobenius(transpose(A)*x-b) / frobenius(b);
    myra::out() << "  |A^T*x-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test hermitian(A)*x = b
    {
    auto b = Matrix<Number>::random(N,3);
    auto x = b;
    solver.refine(x,'L','H');
    Precision residual = frobenius(hermitian(A)*x-b) / frobenius(b);
    myra::out() << "  |A^H*x-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test conjugate(A)*x = b
    {
    auto b = Matrix<Number>::random(N,3);
    auto x = b;
    solver.refine(x,'L','C');
    Precision residual = frobenius(conjugate(A)*x-b) / frobenius(b);
    myra::out() << "  |A^C*b-x| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test x*A = b
    {
    auto b = Matrix<Number>::random(3,N);
    auto x = b;
    solver.refine(x,'R','N');
    Precision residual = frobenius(x*A-b) / frobenius(b);
    myra::out() << "  |x*A-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test x*transpose(A) = b
    {
    auto b = Matrix<Number>::random(3,N);
    auto x = b;
    solver.refine(x,'R','T');
    Precision residual = frobenius(x*transpose(A)-b) / frobenius(b);
    myra::out() << "  |x*A^T-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test x*hermitian(A) = b
    {
    auto b = Matrix<Number>::random(3,N);
    auto x = b;
    solver.refine(x,'R','H');
    Precision residual = frobenius(x*hermitian(A)-b) / frobenius(b);
    myra::out() << "  |x*A^H-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Test x*conjugate(A) = b
    {
    auto b = Matrix<Number>::random(3,N);
    auto x = b;
    solver.refine(x,'R','C');
    Precision residual = frobenius(x*conjugate(A)-b) / frobenius(b);
    myra::out() << "  |x*A^C-b| = " << residual << std::endl;
    REQUIRE(residual < tolerance);
    }
  // Make a deep copy of solver, test A*x = b again.
    {
    Solver copy = solver;
    auto b = Matrix<Number>::random(N,3);
    auto x = b;
    solver.refine(x,'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 = Matrix<Number>::random(N,3);
    auto x = b;
    saved.refine(x,'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("lu2_solver","[multifrontal][parallel]")
  {
  // Define problem size.
  int I = 64;
  int J = 64;
  int N = I*J;
  // Run tests.
  test<NumberD>(I,J,1.0e-8);
  test<NumberZ>(I,J,1.0e-8);
  }

Results: [PASS]

double
  |A*b-x| = 3.3134e-13
  |A^T*x-b| = 9.94901e-12
  |A^H*x-b| = 1.09132e-11
  |A^C*b-x| = 6.26222e-12
  |x*A-b| = 2.06765e-13
  |x*A^T-b| = 3.69549e-12
  |x*A^H-b| = 6.07387e-12
  |x*A^C-b| = 1.38604e-11
  |inv(A)*b-x| (copy) = 8.48731e-12
  |A*x-b| (saved) = 6.47891e-12
std::complex<double>
  |A*b-x| = 3.83262e-13
  |A^T*x-b| = 2.95353e-13
  |A^H*x-b| = 4.38668e-13
  |A^C*b-x| = 2.86627e-13
  |x*A-b| = 3.07921e-13
  |x*A^T-b| = 5.08965e-13
  |x*A^H-b| = 3.42837e-13
  |x*A^C-b| = 2.88567e-13
  |inv(A)*b-x| (copy) = 5.57157e-13
  |A*x-b| (saved) = 4.20228e-13

Go back to Summary of /test programs.