zcholesky2_partialsolve

Source: tests/multifrontal/zcholesky/zcholesky2_partialsolve.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/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/sparse/Permutation.h>

// Algorithms.
#include <myramath/utility/random.h>
#include <myramath/utility/sortunique.h>
#include <myramath/dense/gemm.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 <tests/myratest.h>

using namespace myra;

namespace {

template<class Precision> void test(int I, int J, Precision tolerance)
  {
  myra::out() << typestring<Precision>() << std::endl;
  typedef std::complex<Precision> Number;
  // Generate A. 
  int N = I*J;
  SparseMatrix<Number> A = laplacian2<Number>(I,J);
  phasor_hermitian(A);
  // 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(8).set_globsize(4).set_nthreads(1);
  // Reorder A using bisection and factor it.
  Permutation perm = bisect2(I,J);
  SparseZCholeskySolver<Precision> solver(A,perm,options);
  // Make random restriction Ri.
  std::vector<int> Ri(20,0);
  for (int i = 0; i < Ri.size(); ++i)
    Ri[i] = random_int(N);
  sortunique(Ri);
  // Make random restriction Rj.
  std::vector<int> Rj(30,0);
  for (int j = 0; j < Rj.size(); ++j)
    Rj[j] = random_int(N);
  sortunique(Rj);
  // Form restricted inverse Z = Ri'inv(A)*Rj
  Matrix<Number> Z = solver.inverse(Ri,Rj);
  int K = 15;
  // Test from the left.
    {    
    auto B1 = Matrix<Number>::random(ssize(Rj),K);
    auto B2 = Matrix<Number>::random(ssize(Ri),K);
    Precision error_N = frobenius( gemm(Z,'N',B1) - solver.partialsolve(Ri,Rj,B1,'L','N') );
    Precision error_T = frobenius( gemm(Z,'T',B2) - solver.partialsolve(Ri,Rj,B2,'L','T') );
    Precision error_H = frobenius( gemm(Z,'H',B2) - solver.partialsolve(Ri,Rj,B2,'L','H') );
    Precision error_C = frobenius( gemm(Z,'C',B1) - solver.partialsolve(Ri,Rj,B1,'L','C') );
    myra::out() << "  error in Z^N * B = " << error_N << std::endl;
    myra::out() << "  error in Z^T * B = " << error_T << std::endl;
    myra::out() << "  error in Z^H * B = " << error_H << std::endl;
    myra::out() << "  error in Z^C * B = " << error_C << std::endl;
    REQUIRE(error_N < tolerance);
    REQUIRE(error_T < tolerance);
    REQUIRE(error_H < tolerance);
    REQUIRE(error_C < tolerance);
    }
  // Test from the right.
    {
    auto B1 = Matrix<Number>::random(K,ssize(Ri));
    auto B2 = Matrix<Number>::random(K,ssize(Rj));
    Precision error_N = frobenius( gemm(B1,Z,'N') - solver.partialsolve(Ri,Rj,B1,'R','N') );
    Precision error_T = frobenius( gemm(B2,Z,'T') - solver.partialsolve(Ri,Rj,B2,'R','T') );
    Precision error_H = frobenius( gemm(B2,Z,'H') - solver.partialsolve(Ri,Rj,B2,'R','H') );
    Precision error_C = frobenius( gemm(B1,Z,'C') - solver.partialsolve(Ri,Rj,B1,'R','C') );
    myra::out() << "  error in B * Z^N = " << error_N << std::endl;
    myra::out() << "  error in B * Z^T = " << error_T << std::endl;
    myra::out() << "  error in B * Z^H = " << error_H << std::endl;
    myra::out() << "  error in B * Z^C = " << error_C << std::endl;
    REQUIRE(error_N < tolerance);
    REQUIRE(error_T < tolerance);
    REQUIRE(error_H < tolerance);
    REQUIRE(error_C < tolerance);
    }
  }

} // namespace

ADD_TEST("zcholesky2_partialsolve","[multifrontal][parallel]")
  {
  test<float >(20,20,1.0e-4f);
  test<double>(20,20,1.0e-10);
  }

Results: [PASS]

float
  error in Z^N * B = 2.13941e-07
  error in Z^T * B = 2.28013e-07
  error in Z^H * B = 2.37269e-07
  error in Z^C * B = 2.18977e-07
  error in B * Z^N = 1.98323e-07
  error in B * Z^T = 2.2747e-07
  error in B * Z^H = 2.2248e-07
  error in B * Z^C = 2.30425e-07
double
  error in Z^N * B = 4.17196e-16
  error in Z^T * B = 3.94986e-16
  error in Z^H * B = 3.53956e-16
  error in Z^C * B = 3.96072e-16
  error in B * Z^N = 4.18703e-16
  error in B * Z^T = 4.50374e-16
  error in B * Z^H = 4.11423e-16
  error in B * Z^C = 4.22322e-16

Go back to Summary of /test programs.