zldlt2_schur

Source: tests/multifrontal/zldlt/zldlt2_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>
#include <myramath/sparse/phasor.h>

// Solver under test.
#include <myramath/multifrontal/zldlt/SparseZLDLTSolver.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);
  double shift = 5.5;
  for (int n = 0; n < N; ++n)
    A(n,n) -= shift;
  phasor_symmetric(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);
  SparseZLDLTSolver<double> 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(),'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<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(),'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<Number>::random(N,1);
    Number S1 = solver.schur(Bi,Bv,options)(0,0);
    // Compare to reference/dense calculation.
    solver.solveL(Bv);
    auto Cv = Bv;
    solver.solveD(Cv);
    Number S2 = dotu(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);
    }

  // 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(),'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<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(),'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<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);
    solver.solveD(Cv);
    Number S2 = dotu(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("zldlt2_schur","[multifrontal][parallel]")
  {
  test<double>(1.0e-8);
  }

Results: [PASS]

|B'*inv(A)*B [dense] - B'*inv(A)*B [sparse]| = 3.71176e-14
|schur(B)-schur(BiBv)| = 3.44232e-14
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 1.52267e-14
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 3.44319e-14
|schur(B,C)-schur(BiBv,CiCv)| = 3.54289e-14
|schur(Bi,Bv)-dot(L\Bv,L\Bv))| = 1.60637e-15

Go back to Summary of /test programs.