rcholesky2_schur

Source: tests/multifrontal/rcholesky/rcholesky2_schur.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/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/sparse/SparseMatrixBuilder.h>
#include <myramath/sparse/PatternRange.h>
#include <myramath/sparse/Permutation.h>

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

using namespace myra;

namespace {

// Generates random row indices.
std::vector<int> random_rows(int I, int N)
  {
  std::vector<int> answer;
  for (int n = 0; n < N; ++n)
    answer.push_back( random_int(0,I) );
  sortunique(answer);
  return answer;
  }

// Fills SparseMatrix B from (Bi,Bv) data.
template<class Number> SparseMatrix<Number> fill(int I, const std::vector<int>& Bi, const Matrix<Number>& Bv)
  {
  int J = Bv.size().second;
  SparseMatrixBuilder<Number> B_builder(I,J);
  for (int j = 0; j < J; ++j)
    for (int i = 0; i < Bi.size(); ++i)
      B_builder(Bi[i],j) = Bv(i,j);
  return B_builder.make_SparseMatrix();
  }

// Various schur() tests.
template<class Precision> void test(Precision tolerance)
  {
  int I = 20;
  int J = 25;
  int N = I*J;
  int K = 23;
  int K1 = K-1;
  int K2 = K+1;

  // Generate A.  
  auto A = laplacian2<Precision>(I,J);

  // Options related to reordering.
  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);
  SparseRCholeskySolver<Precision> solver(A,perm,options);

  // Testing symmetric sparse schur complement.
    {
    // Generate B, generate a schur complement B'*inv(A)*B
    auto B = SparseMatrix<Precision>::random(N,K,3*N);
    LowerMatrix<Precision> S1 = solver.schur(B,options);
    // Compare to reference/dense calculation.
    Matrix<Precision> S2 = gemm(B.make_Matrix(),'T',gemm(inverse(A.make_Matrix()),B.make_Matrix()));
    Precision S_error = frobenius(S1.make_Matrix('T')-S2)/frobenius(S2);
    myra::out() << "|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  // Testing symmetric restricted schur complement.
    {
    // Compute a symmetric schur complement Bv'*Bi'*inv(A)*Bi*Bv
    auto Bi = random_rows(N,10);
    auto Bv = Matrix<Precision>::random(ssize(Bi),K);
    auto B = fill(N,Bi,Bv);
    auto S1 = solver.schur(Bi,Bv,options);
    // Compare to reference/dense calculation.
    auto S2 = gemm(B.make_Matrix(),'T',gemm(inverse(A.make_Matrix()),B.make_Matrix()));       
    Precision S_error = frobenius(S1.make_Matrix('T')-S2)/frobenius(S2);
    myra::out() << "|schur(B)-schur(BiBv)| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  // Testing symmetric schur complement by a dense Vector, exercises write race resolution.
    {
    // Compute a symmetric schur complement Bv'*Bi'*inv(A)*Bi*Bv
    auto Bi = ilinspace(0,N);
    auto Bv = Matrix<Precision>::random(N,1);
    Precision S1 = solver.schur(Bi,Bv,options)(0,0);
    // Compare to reference/dense calculation.
    solver.solveL(Bv);
    Precision S2 = dot(Bv.vector(0),Bv.vector(0));
    Precision S_error = std::abs(S1-S2)/std::abs(S2);
    myra::out() << "|schur(Bi,Bv)-dot(L\\Bv,L\\Bv))| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  // Testing nonsymmetric schur complement.
    {
    // Generate C/D, compute schur complement C'*inv(A)*D  
    auto C = SparseMatrix<Precision>::random(N,K1,5*N);
    auto D = SparseMatrix<Precision>::random(N,K2,4*N);
    Matrix<Precision> S1 = solver.schur(C,D,options);
    // Compare to reference/dense calculation.
    Matrix<Precision> S2 = gemm(C.make_Matrix(),'T',gemm(inverse(A.make_Matrix()),D.make_Matrix()));
    Precision S_error = frobenius(S1-S2)/frobenius(S2);
    myra::out() << "|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  // Testing unsymmetric restricted schur complement.
    {
    // Compute a symmetric schur complement Bv'*Bi'*inv(A)*Ci*Cv
    auto Bi = random_rows(N,10);
    auto Bv = Matrix<Precision>::random(ssize(Bi),K1);
    auto Ci = random_rows(N,10);
    auto Cv = Matrix<Precision>::random(ssize(Ci),K2);  
    auto S1 = solver.schur(Bi,Bv,Ci,Cv,options);
    // Compare to reference/dense calculation.
    auto B = fill(N,Bi,Bv);
    auto C = fill(N,Ci,Cv);
    auto S2 = gemm(B.make_Matrix(),'T',gemm(inverse(A.make_Matrix()),C.make_Matrix()));   
    Precision S_error = frobenius(S1-S2)/frobenius(S2);
    myra::out() << "|schur(B,C)-schur(BiBv,CiCv)| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  // Testing unsymmetric schur complement by a dense Vector, exercises write race resolution.
    {
    auto Bi = ilinspace(0,N);
    auto Bv = Matrix<Precision>::random(N,1);
    auto Ci = ilinspace(0,N);
    auto Cv = Matrix<Precision>::random(N,1);
    Precision S1 = solver.schur(Bi,Bv,Ci,Cv,options)(0,0);
    // Compare to reference/dense calculation.
    solver.solveL(Bv);
    solver.solveL(Cv);
    Precision S2 = dot(Bv.vector(0),Cv.vector(0));
    Precision S_error = std::abs(S1-S2)/std::abs(S2);
    myra::out() << "|schur(Bi,Bv)-dot(L\\Bv,L\\Bv))| = " << S_error << std::endl;
    REQUIRE(S_error < tolerance);
    }

  }

} // namespace

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

Results: [PASS]

|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = 1.26579e-07
|schur(B)-schur(BiBv)| = 1.50978e-07
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 8.52305e-08
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 2.7761e-07
|schur(B,C)-schur(BiBv,CiCv)| = 1.49129e-07
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 5.43625e-07
|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = 2.23428e-16
|schur(B)-schur(BiBv)| = 3.53397e-16
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 6.16227e-16
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 6.0128e-16
|schur(B,C)-schur(BiBv,CiCv)| = 4.09679e-16
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 2.09545e-16

Go back to Summary of /test programs.