rcholesky2_inverse

Source: tests/multifrontal/rcholesky/rcholesky2_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/dense/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.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/rcholesky/SparseRCholeskySolver.h>

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

using namespace myra;

ADD_TEST("rcholesky2_inverse","[multifrontal][parallel]")
  {
  // Generate A. Need to break these tests apart, there's a wide dynamic 
  // range in times so using the same input size doesn't make the most sense.
  int I = 10;
  int J = 10;
  int N = I*J;
  SparseMatrix<double> A = laplacian2<double>(I,J);
  Pattern P = stencil2(I,J);
  // 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(2);
  // Reorder A using bisection and factor it.
  Permutation perm = bisect2(I,J);
  SparseRCholeskySolver<double> solver(A,perm,options);

  // Compute inverse of A explicitly, for reference. 
  Matrix<double> invA = inverse(A.make_Matrix());

  // Compare to solver.inverse(i,j) for each nonzero original pattern.
    {
    Timer timer;
    double error1 = 0;
    for (int j = 0; j < N; ++j)
      {
      std::vector<int> Pj = P.pattern(j);
      for (auto i = Pj.begin(); i != Pj.end(); ++i)
        {
        double z1 = invA(*i,j);
        double z2 = solver.inverse(*i,j);
        error1 = std::max(error1, scalar_norm1(invA(*i,j)-solver.inverse(*i,j)) );
        }
      }
    double t = timer.elapsed_time();
    myra::out() << "|solver.inverse(i,j)| = " << error1 << ", " << t << " sec" << std::endl;
    REQUIRE(error1 < 1.0e-12);
    }

  // Check the inverse over some arbitrary 20x20 tensor block.
    {
    Timer timer;
    Matrix<double> invA_block = solver.inverse(ilinspace(5,25),ilinspace(27,47));
    // Matrix<double> invA_block(20,20);
    double t = timer.elapsed_time();
    double error3 = frobenius(invA_block - invA.window(5,25,27,47));
    myra::out() << "|solver.inverse(i_indices,j_indices)| = " << error3 << ", " << t << " sec" << std::endl;
    REQUIRE(error3 < 1.0e-12);
    }

  // Check the inverse over a diagonal 20x20 tensor block.
    {
    Timer timer;
    LowerMatrix<double> invA_block = solver.inverse(ilinspace(5,25),options.set_nthreads(1));
    double t = timer.elapsed_time();
    double error4 = frobenius(invA_block.make_Matrix('T') - invA.window(5,25,5,25));
    myra::out() << "|solver.inverse(ij_indices)| = " << error4 << ", " << t << " sec" << std::endl;
    REQUIRE(error4 < 1.0e-12);
    }

  }

Results: [PASS]

|solver.inverse(i,j)| = 3.33067e-16, 0.154009 sec
|solver.inverse(i_indices,j_indices)| = 1.19586e-16, 0.020001 sec
|solver.inverse(ij_indices)| = 4.99847e-16, 0.002 sec

Go back to Summary of /test programs.