lu2_schur

Source: tests/multifrontal/lu/lu2_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/sparse/Pattern.h>
#include <myramath/sparse/PatternRange.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/lu/SparseLUSolver.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();
  }

template<class Number> void test(typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  // Generate A, structured stencil but random values.
  int I = 20;
  int J = 25;
  int N = I*J;
  int K = 23;
  int K1 = K-1;
  int K2 = K+1;

  // Generate A.  
  Pattern P = stencil2(I,J);
  auto A = SparseMatrix<Number>::random(P);

  // 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);
  SparseLUSolver<Number> solver(A,perm,options);
  
    {
    // Generate C/D, compute a schur complement C'*inv(A)*D
    auto C = SparseMatrix<Number>::random(N,K1,5);
    auto D = SparseMatrix<Number>::random(N,K2,4);
    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 < 1.0e-10);
    }

    {
    // Compute a restricted schur complement Bv'*Bi'*inv(A)*Ci*Cv
    auto Bi = random_rows(N,10);
    auto Bv = Matrix<Number>::random(N,K1);
    auto Ci = random_rows(N,10);
    auto Cv = Matrix<Number>::random(N,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);
    }

  }

} // namespace

ADD_TEST("lu2_schur","[multifrontal][parallel]")
  {
  test<NumberD>(1.0e-8);
  test<NumberZ>(1.0e-8);
  }

Results: [PASS]

|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 5.03812e-14
|schur(B,C)-schur(BiBv,CiCv)| = 6.22951e-14
|C'*inv(A)*D [dense] - C'*inv(A)*D [sparse]| = 5.61056e-15
|schur(B,C)-schur(BiBv,CiCv)| = 4.53647e-15

Go back to Summary of /test programs.