lu2_inverse

Source: tests/multifrontal/lu/lu2_inverse.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/intRange.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.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/utility/ilinspace.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/inverse.h>
#include <myramath/dense/frobenius.h>
#include <myramath/sparse/laplacian2.h>

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

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

using namespace myra;

namespace {

template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Generate A. 
  int N = I*J;
  Pattern P = stencil2(I,J);
  auto A = SparseMatrix<Number>::random(P);
  // Options related to reordering. Test input is small, so need to use small 
  // blocks/no globbing to make sure the symbolic factorization is nontrivial
  typedef multifrontal::Options Options;
  Options options = Options::create().set_blocksize(4).set_globsize(4).set_nthreads(4);
  // Reorder A using bisection and factor it.
  Permutation perm = bisect2(I,J);
  SparseLUSolver<Number> solver(A,perm,options);
  // Compute inverse of A explicitly, for reference. 
  Matrix<Number> Z = inverse(A.make_Matrix());
  // Compare to solver.inverse(i,j) for each nonzero original pattern.
    {
    Precision error = 0;
    for (int j = 0; j < N; ++j)
      {
      std::vector<int> Pj = P.pattern(j);
      for (auto i = Pj.begin(); i != Pj.end(); ++i)
        {
        Number Zij = solver.inverse(*i,j);
        Precision error_ij = scalar_norm1(Z(*i,j)-Zij);
        if (error_ij > error)
          error = error_ij;
        }
      }
    myra::out() << "  |solver.inverse(i,j)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  // Check the inverse over some arbitrary 20x20 tensor block.
    {
    int i0 =  5; int i1 = 25;
    int j0 = 27; int j1 = 47;
    std::vector<int> i = ilinspace(i0,i1);
    std::vector<int> j = ilinspace(j0,j1);
    Matrix<Number> Z_block = solver.inverse(i,j,options);
    Precision error = frobenius(Z_block - Z.window(i0,i1,j0,j1)) / frobenius(Z_block);
    myra::out() << "  |solver.inverse(i_indices,j_indices)| = " << error << std::endl;
    REQUIRE(error < tolerance);
    }
  }
} // namespace

ADD_TEST("lu2_inverse","[multifrontal][parallel]")
  {
  int I = 10;
  int J = 10;
  test<NumberD>(I,J,1.0e-10);
  test<NumberZ>(I,J,1.0e-10);  
  }

Results: [PASS]

double
  |solver.inverse(i,j)| = 5.86198e-13
  |solver.inverse(i_indices,j_indices)| = 3.15979e-14
std::complex<double>
  |solver.inverse(i,j)| = 9.56597e-15
  |solver.inverse(i_indices,j_indices)| = 2.3699e-15

Go back to Summary of /test programs.