lu2_constructors

Source: tests/multifrontal/lu/lu2_constructors.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/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/dense/euclidean.h>
#include <myramath/sparse/gemv.h>
#include <myramath/sparse/laplacian2.h>

// Solver under test.
#include <myramath/multifrontal/lu/SparseLUSolver.h>

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

using namespace myra;

ADD_TEST("lu2_constructors","[multifrontal][parallel]")
  {
  typedef double Number;
  typedef double Precision;
  // Make laplacian matrix A.  
  int I = 64;
  int J = 64;
  int N = I*J;
  Pattern P = stencil2(I,J);
  auto A = SparseMatrix<Number>::random(P);
  // Define solver options.
  typedef SparseLUSolver<Number> Solver;
  typedef ::multifrontal::Options Options;
  Options options = Options::create().set_blocksize(16).set_globsize(4).set_nthreads(1); 
  // Three construction options ..
  Permutation perm = bisect2(I,J);
  Solver solver1(A,options);                // .. factor A, using internal reorder()'ing
  Solver solver2(A,perm,options);           // .. factor A, using caller-provided ordering
  Solver solver3(A,solver1.tree(),options); // .. factor A, reusing symbolic factorization
  // Make random b and three solution vectors.
  auto b = Vector<Precision>::random(N);
  auto x1 = b;
  auto x2 = b;
  auto x3 = b;
  // Solve system using each solver.
  solver1.refine(x1.column());
  solver2.refine(x2.column());
  solver3.refine(x3.column());
  // Check residuals.
  Precision error1 = euclidean(A*x1-b);
  Precision error2 = euclidean(A*x2-b);
  Precision error3 = euclidean(A*x3-b);
  myra::out() << "|A*x1-b| = " << error1 << std::endl;
  myra::out() << "|A*x2-b| = " << error2 << std::endl;
  myra::out() << "|A*x3-b| = " << error3 << std::endl;
  }

Results: [PASS]

|A*x1-b| = 3.00701e-12
|A*x2-b| = 3.22009e-12
|A*x3-b| = 3.00701e-12

Go back to Summary of /test programs.