doc/html/test_rcholesky2_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/rcholesky/SparseRCholeskySolver.h>

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

 using namespace myra;

 namespace {

 template<class Precision> void test(int I, int J, Precision tolerance)
   {
   myra::out() << typestring<Precision>() << std::endl;
   // Generate A.
   int N = I*J;
   auto A = laplacian2<Precision>(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(8).set_globsize(4).set_nthreads(1);
   // Reorder A using bisection and factor it.
   Permutation perm = bisect2(I,J);
   SparseRCholeskySolver<Precision> 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<Precision> Z = solver.inverse(Ri,Rj);
   int K = 15;
   // Test from the left.
     {
     auto B1 = Matrix<Precision>::random(ssize(Rj),K);
     auto B2 = Matrix<Precision>::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<Precision>::random(K,ssize(Ri));
     auto B2 = Matrix<Precision>::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("rcholesky2_partialsolve","[multifrontal][parallel]")
   {
   test<float >(20,20,1.0e-4f);
   test<double>(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< Precision >

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.

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

SparseRCholeskySolver.h
Sparse direct solver suitable for real symmetric positive definite systems.

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

myra::SparseRCholeskySolver
Sparse direct solver suitable for real symmetric positive definite systems.
Definition: SparseRCholeskySolver.h:61

random.h
Simplistic random number functions.

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.