rldlt2_partialsolve

Source: tests/multifrontal/rldlt/rldlt2_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>

// Solver under test.
#include <myramath/multifrontal/rldlt/SparseRLDLTSolver.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;
  // Generate A. 
  int N = I*J;
  SparseMatrix<Precision> A = laplacian2<Precision>(I,J);
  Precision shift = 5.5;
  for (int n = 0; n < N; ++n)
    A(n,n) -= shift;
  // 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);
  SparseRLDLTSolver<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<Precision> Z = solver.inverse(Ri,Rj);
  int K = 15;
  // Test from the left.
    {    
    auto B1 = Matrix<Precision>::random(ssize(Rj),K);
    auto B2 = Matrix<Precision>::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<Precision>::random(K,ssize(Ri));
    auto B2 = Matrix<Precision>::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("rldlt2_partialsolve","[multifrontal][parallel]")
  { test<double>(20,20,1.0e-8); }

Results: [PASS]

double
  error in Z^N * B = 1.03766e-13
  error in Z^T * B = 1.22516e-13
  error in Z^H * B = 1.22516e-13
  error in Z^C * B = 1.03766e-13
  error in B * Z^N = 1.19607e-13
  error in B * Z^T = 1.27355e-13
  error in B * Z^H = 1.27355e-13
  error in B * Z^C = 1.19607e-13

Go back to Summary of /test programs.