doc/html/test_lu2_partialsolve.html

 // ========================================================================= //
 // 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/SparseMatrix.h>
 #include <myramath/sparse/SparseMatrixRange.h>
 #include <myramath/sparse/Permutation.h>

 // Algorithms.
 #include <myramath/utility/random.h>
 #include <myramath/utility/sortunique.h>
 #include <myramath/dense/gemm.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 {

 template<class Number> void test(int I, int J, typename ReflectPrecision<Number>::type tolerance)
   {
   myra::out() << typestring<Number>() << std::endl;
   typedef typename ReflectPrecision<Number>::type Precision;
   // Generate A.
   int N = I*J;
   SparseMatrix<Number> A = laplacian2<Number>(I,J);
   // Options related to reordering. Test input is small, so need to use small
   // blocks/no globbing to make sure the symbolic factorization is nontrivial
   typedef multifrontal::Options Options;
   Options options = Options::create().set_blocksize(4).set_globsize(4).set_nthreads(1);
   // Reorder A using bisection and factor it.
   Permutation perm = bisect2(I,J);
   SparseLUSolver<Number> solver(A,perm,options);
   // Make random restriction Ri.
   std::vector<int> Ri(20,0);
   for (int i = 0; i < Ri.size(); ++i)
     Ri[i] = random_int(N);
   sortunique(Ri);
   // Make random restriction Rj.
   std::vector<int> Rj(30,0);
   for (int j = 0; j < Rj.size(); ++j)
     Rj[j] = random_int(N);
   sortunique(Rj);
   // Form restricted inverse Z = Ri'inv(A)*Rj
   Matrix<Number> Z = solver.inverse(Ri,Rj);
   int K = 15;
   // Test from the left.
     {
     auto B1 = Matrix<Number>::random(ssize(Rj),K);
     auto B2 = Matrix<Number>::random(ssize(Ri),K);
     Precision error_N = frobenius( gemm(Z,'N',B1) - solver.partialsolve(Ri,Rj,B1,'L','N') );
     Precision error_T = frobenius( gemm(Z,'T',B2) - solver.partialsolve(Ri,Rj,B2,'L','T') );
     Precision error_H = frobenius( gemm(Z,'H',B2) - solver.partialsolve(Ri,Rj,B2,'L','H') );
     Precision error_C = frobenius( gemm(Z,'C',B1) - solver.partialsolve(Ri,Rj,B1,'L','C') );
     myra::out() << "  error in Z^N * B = " << error_N << std::endl;
     myra::out() << "  error in Z^T * B = " << error_T << std::endl;
     myra::out() << "  error in Z^H * B = " << error_H << std::endl;
     myra::out() << "  error in Z^C * B = " << error_C << std::endl;
     REQUIRE(error_N < tolerance);
     REQUIRE(error_T < tolerance);
     REQUIRE(error_H < tolerance);
     REQUIRE(error_C < tolerance);
     }
   // Test from the right.
     {
     auto B1 = Matrix<Number>::random(K,ssize(Ri));
     auto B2 = Matrix<Number>::random(K,ssize(Rj));
     Precision error_N = frobenius( gemm(B1,Z,'N') - solver.partialsolve(Ri,Rj,B1,'R','N') );
     Precision error_T = frobenius( gemm(B2,Z,'T') - solver.partialsolve(Ri,Rj,B2,'R','T') );
     Precision error_H = frobenius( gemm(B2,Z,'H') - solver.partialsolve(Ri,Rj,B2,'R','H') );
     Precision error_C = frobenius( gemm(B1,Z,'C') - solver.partialsolve(Ri,Rj,B1,'R','C') );
     myra::out() << "  error in B * Z^N = " << error_N << std::endl;
     myra::out() << "  error in B * Z^T = " << error_T << std::endl;
     myra::out() << "  error in B * Z^H = " << error_H << std::endl;
     myra::out() << "  error in B * Z^C = " << error_C << std::endl;
     REQUIRE(error_N < tolerance);
     REQUIRE(error_T < tolerance);
     REQUIRE(error_H < tolerance);
     REQUIRE(error_C < tolerance);
     }
   }

 } // namespace

 ADD_TEST("lu2_partialsolve","[multifrontal][parallel]")
   {
   test<NumberD>(20,20,1.0e-10);
   test<NumberZ>(20,20,1.0e-10);
   }
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.

myra::multifrontal::Options
Options pack for routines in /multifrontal.
Definition: Options.h:24

myra::Permutation
Represents a Permutation matrix, used to reorder rows/columns/etc of various numeric containers...
Definition: Permutation.h:34

myra::Matrix
Tabulates an IxJ matrix. Allows random access, has column major layout to be compatible with BLAS/LAP...
Definition: bdsqr.h:20

frobenius.h
Routines for computing Frobenius norms of various algebraic containers.

myra::Matrix::random
static Matrix< Number > random(int I, int J)
Generates a random Matrix of specified size.
Definition: Matrix.cpp:353

sortunique.h
Reduces a std::vector to its unique entries, and sorts it.

SparseMatrix.h
General purpose compressed-sparse-column (CSC) container.

myra::sortunique
void sortunique(std::vector< T > &v)
Reduces a std::vector to its unique entries, and sorts it.
Definition: sortunique.h:20

myra
Definition: syntax.dox:1

Number.h
Various utility functions/classes related to scalar Number types.

SparseLUSolver.h
Sparse direct solver suitable for symmetric-pattern nonsymmetric-value A.

Matrix.h
General purpose dense matrix container, O(i*j) storage.

myra::SparseLUSolver
Sparse direct solver suitable for symmetric-pattern nonsymmetric-valued A.
Definition: SparseLUSolver.h:57

myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33

Permutation.h
Aggregates a (perm, iperm, swaps) triple into a vocabulary type.

random.h
Simplistic random number functions.

myra::SparseMatrix
Stores an IxJ matrix A in compressed sparse column format.
Definition: bothcat.h:23

laplacian2.h
Helper routines for reordering/filling 2D structured grids. Used by many unit tests.

gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.

SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.

intRange.h
Interface class for representing subranges of contiguous int&#39;s.