rldlt2_schur

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

// Solver under test.
#include <myramath/multifrontal/rldlt/SparseRLDLTSolver.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.
  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.
  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);
  SparseRLDLTSolver<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);
    auto Cv = Bv;
    solver.solveD(Bv);
    Precision S2 = dot(Cv.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);
    solver.solveD(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("rldlt2_schur","[multifrontal][parallel]")
  {
  test<double>(1.0e-8);
  }

Results: [PASS]

|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = 2.47397e-14
|schur(B)-schur(BiBv)| = 2.71079e-14
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 1.11263e-14
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 2.28618e-14
|schur(B,C)-schur(BiBv,CiCv)| = 2.06378e-14
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 1.84857e-16

Go back to Summary of /test programs.