rldlt2_saddle1

Source: tests/multifrontal/rldlt/rldlt2_saddle1.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/sparse/Pattern.h>
#include <myramath/sparse/PatternRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/sparse/Permutation.h>

// Algorithms.
#include <myramath/utility/sortunique.h>
#include <myramath/dense/euclidean.h>
#include <myramath/sparse/laplacian2.h>
#include <myramath/sparse/transpose.h>
#include <myramath/sparse/bothcat.h>
#include <myramath/sparse/gemv.h>
#include <myramath/multifrontal/detail/restriction.h>
#include <myramath/multifrontal/symbolic/reorder.h>
#include <myramath/multifrontal/symbolic/mmd.h>

// Solver under test.
#include <myramath/multifrontal/rldlt/SparseRLDLTSolver.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;

namespace detail {

template<class Precision>
void test_saddle1(const SparseMatrix<Precision>& A, const Permutation& P, multifrontal::Options options, Precision tolerance)
  {
  SparseRLDLTSolver<Precision> solver(A,P,options);
  // Check A*x = b solution.
  int N = solver.size();
  auto b = Vector<Precision>::random(N);
  auto x = b;
  solver.solve(x.column());
  Precision error = euclidean( gemv(A,x) - b );
  myra::out() << "|A*x-b| = " << error << std::endl;
  myra::out() << "A-flops = " << solver.tree().n_work_llt() << std::endl;
  REQUIRE( error < tolerance );
  }

template<class Precision>
void test_saddle1(int I, int J, int K, Precision tolerance)
  {
  // Generate system with saddle point structure, A = [L R; R' 0], where L is 
  // a graph laplacian and R is some restriction of dirichlet constraints.
  // .. laplacian part L, sized IxJ
  int N = I*J;
  SparseMatrix<Precision> L = laplacian2<Precision>(I,J);
  // .. restriction part R, enforces K random constraints
  std::vector<int> Ri;
  for (int k = 0; k < K; ++k)
    Ri.push_back( random_int(N) );
  sortunique(Ri);
  K = ssize(Ri);
  using namespace multifrontal::detail;
  SparseMatrix<Precision> R = restriction<Precision>(N,Ri);
  // .. transpose part R'
  SparseMatrix<Precision> Rt = transpose(R);
  // .. zero part 0
  SparseMatrix<Precision> Z = SparseMatrix<Precision>::zeros(K,K);  
  // Concatenate system.
  SparseMatrix<Precision> A = bothcat(L,R,Rt,Z);

  // Check the flop count to solve just L.
  typedef multifrontal::Options Options;
  Options options = Options::create().set_blocksize(1).set_globsize(1).set_nthreads(1); 
  Permutation P = reorder(L.pattern());
  myra::out() << "L-flops = " << AssemblyTree(L.pattern(),P,options).n_work_llt() << std::endl;

  // Run the test with various permutations ..

  // .. P0, natural order
  std::vector<int> natural(N+K,-1);
  for (int n = 0; n < N+K; ++n)
    natural[n] = n;
  Permutation P0 = Permutation::from_iperm(natural);
  test_saddle1(A,P0,options,tolerance);

  // .. P1, reversed order
  std::vector<int> reversed(N+K,-1);
  for (int n = 0; n < N+K; ++n)
    reversed[n] = N+K-1-n;
  Permutation P1 = Permutation::from_iperm(reversed);
  test_saddle1(A,P1,options,tolerance);

  // .. P2, from reorder()
  Permutation P2 = reorder(A.pattern());
  test_saddle1(A,P2,options,tolerance);

  // .. P3, from mmd()
  Permutation P3 = mmd(A.pattern());
  test_saddle1(A,P3,options,tolerance);

  }

} // namespace detail

ADD_TEST("rldlt2_saddle1","[multifrontal]")
  {
  int IJK = 20;
  ::detail::test_saddle1<float >(IJK,IJK,IJK, 1.0e-4f);  
  ::detail::test_saddle1<double>(IJK,IJK,IJK, 1.0e-10);
  }

Results: [PASS]

L-flops = 33577
|A*x-b| = 2.60667e-06
A-flops = 98773
|A*x-b| = 2.62035e-06
A-flops = 82951
|A*x-b| = 2.15171e-06
A-flops = 36449
|A*x-b| = 2.086e-06
A-flops = 47682
L-flops = 34879
|A*x-b| = 5.19611e-15
A-flops = 100000
|A*x-b| = 5.75119e-15
A-flops = 83614
|A*x-b| = 4.60797e-15
A-flops = 42835
|A*x-b| = 4.1312e-15
A-flops = 54197

Go back to Summary of /test programs.