zcholesky2_schur

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

// Solver under test.
#include <myramath/multifrontal/zcholesky/SparseZCholeskySolver.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)
  {
  typedef typename ReflectComplex<Precision>::type Number;
  int I = 20;
  int J = 25;
  int N = I*J;
  int K = 23;
  int K1 = K-1;
  int K2 = K+1;

  // Generate A.  
  SparseMatrix<Number> A = laplacian2<Number>(I,J);
  phasor_hermitian(A);

  // 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);
  SparseZCholeskySolver<Precision> solver(A,perm,options);

  // Testing symmetric sparse schur complement.
    {
    // Generate B, generate a schur complement B'*inv(A)*B
    auto B = SparseMatrix<Number>::random(N,K,3*N);
    LowerMatrix<Number> S1 = solver.schur(B,options);
    // Compare to reference/dense calculation.
    Matrix<Number> S2 = gemm(B.make_Matrix(),'H',gemm(inverse(A.make_Matrix()),B.make_Matrix()));
    Precision S_error = frobenius(S1.make_Matrix('H')-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<Number>::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(),'H',gemm(inverse(A.make_Matrix()),B.make_Matrix()));       
    Precision S_error = frobenius(S1.make_Matrix('H')-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<Number>::random(N,1);
    Number S1 = solver.schur(Bi,Bv,options)(0,0);
    // Compare to reference/dense calculation.
    solver.solveL(Bv);
    Number 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<Number>::random(N,K1,5*N);
    auto D = SparseMatrix<Number>::random(N,K2,4*N);
    Matrix<Number> S1 = solver.schur(C,D,options);
    // Compare to reference/dense calculation.
    Matrix<Number> S2 = gemm(C.make_Matrix(),'H',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<Number>::random(ssize(Bi),K1);
    auto Ci = random_rows(N,10);
    auto Cv = Matrix<Number>::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(),'H',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<Number>::random(N,1);
    auto Ci = ilinspace(0,N);
    auto Cv = Matrix<Number>::random(N,1);
    Number S1 = solver.schur(Bi,Bv,Ci,Cv,options)(0,0);
    // Compare to reference/dense calculation.
    solver.solveL(Bv);
    solver.solveL(Cv);
    Number 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("zcholesky2_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.55229e-07
|schur(B)-schur(BiBv)| = 1.94353e-07
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 6.00614e-07
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 3.05684e-07
|schur(B,C)-schur(BiBv,CiCv)| = 1.92332e-07
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 6.13198e-07
|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = 2.46697e-16
|schur(B)-schur(BiBv)| = 4.39655e-16
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 5.00086e-16
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 6.05179e-16
|schur(B,C)-schur(BiBv,CiCv)| = 8.87199e-16
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 1.15114e-15

Go back to Summary of /test programs.