zcholesky2_inverse

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

// Solver under test.
#include <myramath/multifrontal/zcholesky/SparseZCholeskySolver.h>

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

#include <algorithm>

using namespace myra;

ADD_TEST("zcholesky2_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<NumberZ> A = laplacian2<NumberZ>(I,J);
  phasor_hermitian(A);
  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(1);
  // Reorder A using bisection and factor it.
  Permutation perm = bisect2(I,J);
  SparseZCholeskySolver<double> solver(A,perm,options);
  // Compute inverse of A explicitly, for reference. 
  Matrix<NumberZ> 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)
        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<NumberZ> invA_block = solver.inverse(ilinspace(5,25), ilinspace(27,47));
    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<NumberZ> invA_block = solver.inverse(ilinspace(5,25));
    double t = timer.elapsed_time();
    double error4 = frobenius(invA_block.make_Matrix('H') - 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.88594e-16, 0.098005 sec
|solver.inverse(i_indices,j_indices)| = 1.58692e-16, 0.041003 sec
|solver.inverse(ij_indices)| = 5.03753e-16, 0.025001 sec

Go back to Summary of /test programs.