Tutorial for /dense routines

Overview

The /dense folder contains basic containers for tabulated numeric data like Matrix and Vector, as well as specialized containers for numeric data with exploitable nonzero structure, like DiagonalMatrix and LowerMatrix. Layered on top of these containers are various Range types, which are views into numeric data that is owned by some other underlying container. Algorithms complete the picture, by performing useful linear algebraic operations (arithmetic, decompositions, linear solution, eigensolution) upon Range's.

This layered approach (containers/ranges/algorithms) imitates the design of C++ Standard Template Library (containers/iterators/algorithms) and is intended to provide a measure of container-agnosticism to algorithms. Note that most of the heavy lifting within /dense is performed by the underlying BLAS and LAPACK libraries.

Below are links to the sections of this page:

• Matrix and MatrixRange
• BLAS1 and BLAS2 Routines
• BLAS3 Routines
• LAPACK Routines
• Solver Classes
• Appendix: Working with external data
• Appendix: A note about const-correctness.

Matrix and MatrixRange

Of all the containers in the library, Matrix is the one you will probably use the most. It tabulates numeric data into an IxJ rectangle. The internal storage is column-major to remain backwards compatible with Fortran, LAPACK and the BLAS. Matrix indices always start at zero. In fact, every container in the library uses a 0-based indexing convention. Matrix has many member functions for generating a subview (a MatrixRange) over it. Consider tutorial/dense/range.cpp, listed below, which shows some of them in action:

#include <myramath/dense/Matrix.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_range","[dense]") // int main()
  {
  // Make a Matrix A ..
  auto A = Matrix<double>::fill_rmajor(3,3,{1,2,3, 4,5,6, 7,8,9});
  // .. inspect various "views" of A ..
  std::cout << "A = "                 << A                 << std::endl; // all
  std::cout << "A.row(0) = "          << A.row(0)          << std::endl; // first row
  std::cout << "A.column(1) = "       << A.column(1)       << std::endl; // middle column
  std::cout << "A.left(2) = "         << A.left(2)         << std::endl; // left 2 columns
  std::cout << "A.window(1,3,1,3) = " << A.window(1,3,1,3) << std::endl; // trailing 2x2 submatrix
  std::cout << "A.top(2) = "          << A.top(2)          << std::endl; // top 2 rows
  std::cout << "A.top(2).left(2) = "  << A.top(2).left(2)  << std::endl; // leading 2x2 submatrix
  // .. modify A through a range ..
  A.top(2).left(2)(1,1) = 0;
  std::cout << "A.top(2).left(2) = " << A.top(2).left(2) << std::endl;
  std::cout << "A = "                << A                << std::endl;
  }

(Click here to view the output of this program)

A = size 3 by 3 Matrix of double:
[ 1 2 3 ]
[ 4 5 6 ]
[ 7 8 9 ]

A.row(0) = size 1 by 3 Matrix of double:
[ 1 2 3 ]

A.column(1) = size 3 by 1 Matrix of double:
[ 2 ]
[ 5 ]
[ 8 ]

A.left(2) = size 3 by 2 Matrix of double:
[ 1 2 ]
[ 4 5 ]
[ 7 8 ]

A.window(1,3,1,3) = size 2 by 2 Matrix of double:
[ 5 6 ]
[ 8 9 ]

A.top(2) = size 2 by 3 Matrix of double:
[ 1 2 3 ]
[ 4 5 6 ]

A.top(2).left(2) = size 2 by 2 Matrix of double:
[ 1 2 ]
[ 4 5 ]

A.top(2).left(2) = size 2 by 2 Matrix of double:
[ 1 2 ]
[ 4 0 ]

A = size 3 by 3 Matrix of double:
[ 1 2 3 ]
[ 4 0 6 ]
[ 7 8 9 ]

The last few lines of this example show that:

• MatrixRange's can be further partitioned into other MatrixRange's by method chaining.
• Accessors for peeking and poking scalar values use "function call syntax", A(i,j) = ..whatever..
• Modifying values in a MatrixRange causes change in the underlying Matrix. MatrixRange's are views, not deep copies.

BLAS1 and BLAS2 Routines

The next container of interest is Vector, which has an associated VectorRange for representing subviews. The usual BLAS1 vector operations like inner products (dense/dot.h), norms (dense/euclidean.h) and addition (dense/axpy.h) are present and delegate directly to the BLAS. Note that all these algorithms accept VectorRange's, so it's easy to use "just a piece" of a Vector as an operand. Of course, any Vector is implicitly convertible to a VectorRange over all of its contents, so algorithms work with regular Vector's too. The listing below, tutorial/dense/axpy.cpp, shows these algorithms in action.

#include <myramath/dense/Vector.h>
#include <myramath/dense/dot.h>
#include <myramath/dense/axpy.h>
#include <myramath/dense/euclidean.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_axpy","[dense]") // int main()
  {
  // Make a few vectors ..
  auto x = Vector<double>::fill({1,2,3});
  auto y = Vector<double>::fill({4,5,6});
  // .. try out a few operations ..
  std::cout << "x = " << x << std::endl;
  std::cout << "x = " << y << std::endl;
  std::cout << "euclidean(x) = " << euclidean(x) << std::endl;
  std::cout << "dot(x,y) = " << dot(x,y) << std::endl;
  // .. compare axpy to overloaded operators ..
  auto z = Vector<double>::fill({0,0,0});
  axpy(2,x,z);
  axpy(3,y,z);
  std::cout << "z = " << z << std::endl;
  std::cout << "2*x+3*y = " << 2*x+3*y << std::endl;
  // .. try using dot() over ranges ..
  auto w = Vector<double>::fill({1,2,3,4,5,6});
  std::cout << "w = " << w << std::endl;
  std::cout << "w(0:3)'w(3:6) = " << dot(w.first(3),w.last(3)) << std::endl;
  // .. scale just a part of w ..
  w.first(4) *= 2;
  std::cout << "w = " << w << std::endl;
  }

(Click here to view the output of this program)

x = size 3 Vector of double:
[ 1 2 3 ]

x = size 3 Vector of double:
[ 4 5 6 ]

euclidean(x) = 3.74166
dot(x,y) = 32
z = size 3 Vector of double:
[ 14 19 24 ]

2*x+3*y = size 3 Vector of double:
[ 14 19 24 ]

w = size 6 Vector of double:
[ 1 2 3 4 5 6 ]

w(0:3)'w(3:6) = 32
w = size 6 Vector of double:
[ 2 4 6 8 5 6 ]

Vector's can be initialized from MatrixRange's that have size Nx1 or 1xN, such as the .row(i) and .column(j) subviews of a Matrix. When it comes to the question "is this a row Vector, or a column Vector?", the class is agnostic. It's just an N-vector. The "row-ness" or "column-ness" is only deduced from context. For instance, in the code excerpt below (tutorial/dense/gemv.cpp), the gemv(A,x) call will treat x as a column vector (size 3x1), while gemv(x,A) will treat x as a row vector (size 1x3).

#include <myramath/dense/Matrix.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/gemv.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_gemv","[dense]") // int main()
  {
  // Make a Matrix A and a Vector x ..
  auto A = Matrix<double>::fill_rmajor(3,3,{1,4,7, 2,5,8, 3,6,9});
  Vector<double> x(A.column(2));
  std::cout << "A = " << A << std::endl;
  std::cout << "x = " << x << std::endl;
  //
  //                   [1 4 7]   [7]   [102]
  // .. multiply A*x = [2 5 8] * [8] = [126]
  //                   [3 6 9]   [9]   [150]
  //
  std::cout << "A*x = " << gemv(A,x) << std::endl;
  //
  //                             [1 4 7]
  // .. multiply x*A = [7 8 9] * [2 5 8] = [50 122 194]
  //                             [3 6 9]
  //
  std::cout << "x*A = " << gemv(x,A) << std::endl;
  }

(Click here to view the output of this program)

A = size 3 by 3 Matrix of double:
[ 1 4 7 ]
[ 2 5 8 ]
[ 3 6 9 ]

x = size 3 Vector of double:
[ 7 8 9 ]

A*x = size 3 Vector of double:
[ 102 126 150 ]

x*A = size 3 Vector of double:
[ 50 122 194 ]

BLAS3 Routines

BLAS3 routines, which perform O(n^3) work upon O(n^2) data, typically provide the highest level of performance on modern CPU's due to improved data reuse from cache. The most well-known routine in this class is dense/gemm.h, which performs a fused multiply-add C=alpha*op(A)*op(B)+beta*C for MatrixRange's C, A and B. To avoid the need for data movement and temporaries, gemm() allows for various modifiers ("operations") to be applied to the A and B inputs. Each is denoted by a single character argument:

• 'T'ranspose an input
• 'C'onjugate an input
• 'H'ermitian transpose an input (that is, transpose and conjugate it)
• Do 'N'othing to an input, just use it as-is

For the record, it's worth pointing out that this library uses different conventions from the BLAS for "op" modifiers: the BLAS has no option for just conjugating, and instead uses 'C' to denote conjugate transpose (where as this library calls that op='H', for hermitian). Use care. The example below (tutorial/dense/gemm.cpp) exercises all of these options for complex operands and serves as your reference for expected behavior.

#include <myramath/dense/Matrix.h>
#include <myramath/dense/gemm.h>

#include <tests/myratest.h>

#include <complex>
#include <iostream>

using namespace myra;

ADD_TEST("dense_gemm","[dense]") // int main()
  {
  typedef std::complex<double> Z;
  auto A = Matrix<Z>::fill_cmajor(2,2,{Z(1,2),Z(3,4),Z(5,6),Z(7,8)});
  auto B = Matrix<Z>::fill_cmajor(2,2,{Z(9,10),Z(11,12),Z(13,14),Z(15,16)});
  std::cout << "A = " << A << std::endl;
  std::cout << "B = " << B << std::endl;
  std::cout << "A*B = " << gemm(A,'N',B,'N') << std::endl;
  std::cout << "transpose(A)*B = " << gemm(A,'T',B,'N') << std::endl;
  std::cout << "hermitian(A)*B = " << gemm(A,'H',B,'N') << std::endl;
  std::cout << "conjugate(A)*B = " << gemm(A,'C',B,'N') << std::endl;
  }

(Click here to view the output of this program)

A = size 2 by 2 Matrix of std::complex<double>:
[ (1,2) (5,6) ]
[ (3,4) (7,8) ]

B = size 2 by 2 Matrix of std::complex<double>:
[ (9,10) (13,14) ]
[ (11,12) (15,16) ]

A*B = size 2 by 2 Matrix of std::complex<double>:
[ (-28,154) (-36,210) ]
[ (-32,238) (-40,326) ]

transpose(A)*B = size 2 by 2 Matrix of std::complex<double>:
[ (-26,108) (-34,148) ]
[ (-34,276) (-42,380) ]

hermitian(A)*B = size 2 by 2 Matrix of std::complex<double>:
[ (110,-16) (150,-24) ]
[ (278,-8) (382,-16) ]

conjugate(A)*B = size 2 by 2 Matrix of std::complex<double>:
[ (156,-14) (212,-22) ]
[ (240,-10) (328,-18) ]

Another noteworthy BLAS3 routine is dense/syrk.h, which is worth exploring because it leads to a new container type. Short for "symmetric rank-K" update, B=syrk(A) is semantically equivalent to B=gemm(A,A,'T') but performs roughly half the work by exploiting the fact that the output is known to be symmetric. The BLAS-provided syrk() routines write into a square output matrix, but will touch only one of its triangles ('U'pper or 'L'ower) as specified by the caller.

Typically the untouched triangle of syrk()'s output just goes along for the ride as wasted space. However LAPACK working note #199 describes an alternative that is waste-free and BLAS3-friendly called "rectangular packed format", which is encapsulated here as a specialized container named LowerMatrix. The code excerpt below (tutorial/dense/syrk.cpp) explores the relationships between calling (i) gemm(A,A'T') on a Matrix, (ii) syrk() on the 'L'ower triangle of a Matrix, and (iii) syrk() directly on a LowerMatrix:

#include <myramath/dense/Matrix.h>
#include <myramath/dense/LowerMatrix.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/syrk.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_syrk","[dense]") // int main()
  {
  // Make a Matrix A and empty outputs B ..
  auto A = Matrix<double>::fill_rmajor(3,3,{1,4,7, 2,5,8, 3,6,9});
  auto B1 = Matrix<double>::zeros(3,3);
  auto B2 = Matrix<double>::zeros(3,3);  
  auto B3 = LowerMatrix<double>::zeros(3);
  // .. compare gemm_inplace() ..
  gemm_inplace(B1,A,'N',A,'T');
  std::cout << "B1 = gemm_inplace(A,'N',A,'T') = " << B1 << std::endl;
  // .. to syrk_inplace() onto the lower triangle of a Matrix..
  syrk_inplace(B2,'L',A);
  std::cout << "B2 = syrk_inplace('L',A) = " << B2 << std::endl;
  // .. and to syrk_inplace() onto a LowerMatrix ..  
  syrk_inplace(B3,A);
  std::cout << "B3 = syrk_inplace(A) = " << B3 << std::endl;
  // Note the return type of (non-inplace) syrk() is a LowerMatrix ..
  std::cout << "syrk(A) = " << syrk(A) << std::endl;
  // .. but can expand it back out into a Matrix, both triangles filled.
  std::cout << "syrk(A).make_Matrix('T') = " << syrk(A).make_Matrix('T') << std::endl;
  }

(Click here to view the output of this program)

B1 = gemm_inplace(A,'N',A,'T') = size 3 by 3 Matrix of double:
[ 66 78 90 ]
[ 78 93 108 ]
[ 90 108 126 ]

B2 = syrk_inplace('L',A) = size 3 by 3 Matrix of double:
[ 66 0 0 ]
[ 78 93 0 ]
[ 90 108 126 ]

B3 = syrk_inplace(A) = size 3 by 3 LowerMatrix of double:
[ 66 * * ]
[ 78 93 * ]
[ 90 108 126 ]

syrk(A) = size 3 by 3 LowerMatrix of double:
[ 66 * * ]
[ 78 93 * ]
[ 90 108 126 ]

syrk(A).make_Matrix('T') = size 3 by 3 Matrix of double:
[ 66 78 90 ]
[ 78 93 108 ]
[ 90 108 126 ]

The BLAS3 routine dense/trsm.h is very relevant for linear system solution using LU-like decompositions. It performs triangular system backsolution across multiple right hand sides. Just like syrk(), trsm() can operate upon either the 'U'pper/'L'ower triangle of a Matrix, or a LowerMatrix directly. Note that trsm() is highly reconfigurable, being able to solve L*X=B "from the left", X*L=B "from the right", and can also apply an op modifier ('N','T','H' or 'C') to L during the process.

The remaining BLAS3 routines are:

• Hermitian rank-K updates, see dense/herk.h, very similar to syrk()
• Multiplication by a triangular matrix, see dense/trmm.h, essentially the reverse operation to trsm()
• A specialized multiplication routine for a symmetric matrix, see dense/symm.h
• A specialized multiplication routine for a hermitian matrix, see dense/hemm.h

LAPACK Routines

Building on the high performance basic blocks of the BLAS, LAPACK provides a collection of matrix decompositions that can be used for (i) linear system solution, (ii) least squares solution and (iii) spectral representations. LAPACK provides a huge amount of functionality, this library really only exposes the routines that are relevant for sparse linear system solution. Many of the routines operate "in-place" and overwrite their input data. The example program below (tutorial/dense/potrf.cpp) uses potrf_inplace() ("positive triangular factorization", the Cholesky algorithm) to solve an A*X=B system. LowerMatrix plays a large role in this example because the system A is symmetric.

#include <myramath/dense/Matrix.h>
#include <myramath/dense/LowerMatrix.h>
#include <myramath/dense/syrk.h>
#include <myramath/dense/symm.h>
#include <myramath/dense/potrf.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/frobenius.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_potrf","[dense]") // int main()
  {
  // Make arbitrary G ..
  auto G = Matrix<double>::fill_rmajor(3,3,{1,1,1, 1,2,3, 1,3,6});
  std::cout << "G = " << G << std::endl;
  // .. form a symmetric positive A using A=G*G' (uses syrk)..
  LowerMatrix<double> A = syrk(G);
  std::cout << "A = syrk(G) = " << A << std::endl;
  // .. form arbitrary solution X ..
  auto X = Matrix<double>::fill_rmajor(3,2,{4,6, 3,2, 5,1});
  std::cout << "X = " << X << std::endl;
  // .. form the right hand sides B using B=A*X (uses symm)..
  Matrix<double> B = symm('L',A,X);
  std::cout << "B = symm(A,X) = " << B << std::endl;
  // .. factor A = L*L', overwriting A with L (uses potrf) ..
  potrf_inplace(A);
  std::cout << "L = potrf(A) = " << A << std::endl;
  // .. solve L*L'*X = B, overwriting B with X (uses trsm, twice) ..
  trsm_inplace('L','N',A,B);
  trsm_inplace('L','T',A,B);
  // .. check result
  std::cout << "A\\B = L\\L'\\B = " << B << std::endl;
  std::cout << "frobenius(A\\B-X) = " << frobenius(B-X) << std::endl;
  }

(Click here to view the output of this program)

G = size 3 by 3 Matrix of double:
[ 1 1 1 ]
[ 1 2 3 ]
[ 1 3 6 ]

A = syrk(G) = size 3 by 3 LowerMatrix of double:
[ 3 * * ]
[ 6 14 * ]
[ 10 25 46 ]

X = size 3 by 2 Matrix of double:
[ 4 6 ]
[ 3 2 ]
[ 5 1 ]

B = symm(A,X) = size 3 by 2 Matrix of double:
[ 80 40 ]
[ 191 89 ]
[ 345 156 ]

L = potrf(A) = size 3 by 3 LowerMatrix of double:
[ 1.73205 * * ]
[ 3.4641 1.41421 * ]
[ 5.7735 3.53553 0.408248 ]

A\B = L\L'\B = size 3 by 2 Matrix of double:
[ 4 6 ]
[ 3 2 ]
[ 5 1 ]

frobenius(A\B-X) = 1.42081e-12

For sake of completeness, a nonsymmetric example (tutorial/dense/getrf.cpp) based on getrf_inplace() (a pivoted LU decomposition) is shown below. The notable changes are (i) the pivoting step results in row swaps that have to accounted for during backsolution (ii) the lack of symmetry of the problem means that the both the 'L'ower and 'U'pper triangles of the decomposition must be stored, leading to different looking calls to trsm_inplace(). Since pivoting requires additional integer-valued storage to record the row swaps, MyraMath introduces another view-like type (an intRange) to represent this parameter.

#include <myramath/dense/Matrix.h>
#include <myramath/dense/intRange.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/getrf.h>
#include <myramath/dense/swaps.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/frobenius.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_getrf","[dense]") // int main()
  {
  // Make arbitrary A ..
  auto A = Matrix<double>::fill_rmajor(3,3,{1,4,9, 3,6,7, 5,8,2});
  // .. form arbitrary solution X ..
  auto X = Matrix<double>::fill_rmajor(3,2,{10,13, 11,14, 12,15});
  std::cout << "A = " << A << std::endl;
  std::cout << "X = " << X << std::endl;
  // .. form the right hand sides B using B=A*X (uses gemm)..
  Matrix<double> B = gemm(A,X);
  std::cout << "B = A*X = " << B << std::endl;
  // .. factor A = P'*L*U, overwriting A with L/U and returning P (uses getrf_inplace)
  std::vector<int> P(3);
  getrf_inplace(A,P);
  // solve P'*L*U*X=B, overwriting B with X
  swap_rows(P,B);
  trsm_inplace('L','L','N',A,B,'U');
  trsm_inplace('L','U','N',A,B);
  // .. check result
  std::cout << "A\\B = P\\L\\U\\B = " << B << std::endl;
  std::cout << "frobenius(A\\B-X) = " << frobenius(B-X) << std::endl;
  }

(Click here to view the output of this program)

A = size 3 by 3 Matrix of double:
[ 1 4 9 ]
[ 3 6 7 ]
[ 5 8 2 ]

X = size 3 by 2 Matrix of double:
[ 10 13 ]
[ 11 14 ]
[ 12 15 ]

B = A*X = size 3 by 2 Matrix of double:
[ 162 204 ]
[ 180 228 ]
[ 162 207 ]

A\B = P\L\U\B = size 3 by 2 Matrix of double:
[ 10 13 ]
[ 11 14 ]
[ 12 15 ]

frobenius(A\B-X) = 2.7231e-14

The other noteworthy decompositions that are wrapped from LAPACK are:

• Bunch-Kauffman (LDL') factorization of a symmetric indefinite matrix, see dense/sytrf.h
• Bunch-Kauffman (LDL') factorization of a hermitian indefinite matrix, see dense/hetrf.h
• QR decomposition of a general matrix dense/geqrf.h
• LQ decomposition of a general matrix dense/gelqf.h
• Eigen-decomposition of a real symmetric matrix, see dense/syev.h
• Eigen-decomposition of a real symmetric tridiagonal matrix, see dense/stev.h
• Eigen-decomposition of a complex hermitian matrix, see dense/heev.h
• Eigen-decomposition of a general matrix, see dense/geevd.h
• Singular value decomposition of a general matrix, see dense/gesvd.h

This list may grow as needs require.

Solver Classes

LAPACK and the BLAS expose a fixed set of functions that can be pipelined together for solving linear algebra problems. That design is (for the most part) imitated by the routines in /dense. However, one area where slight embellishments have been made is with the introduction of "Solver" classes which encapsulate all the steps (factorization, pivot reordering, backsolution) associated with a given decomposition. Consider tutorial/dense/solver.cpp, below. It solves the same linear system as tutorial/dense/getrf.cpp but at a slightly higher level (using LUSolver() instead of getrf_inplace(), trsm_inplace(), etc).

#include <myramath/dense/Matrix.h>
#include <myramath/dense/LUSolver.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/gemm.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_solver","[dense]") // int main()
  {
  // Make arbitrary A ..
  auto A = Matrix<double>::fill_rmajor(3,3,{1,4,9, 3,6,7, 5,8,2});
  // .. form arbitrary solution X ..
  auto X = Matrix<double>::fill_rmajor(3,2,{10,13, 11,14, 12,15});
  std::cout << "A = " << A << std::endl;
  std::cout << "X = " << X << std::endl;
  // .. form the right hand sides B using B=A*X (uses gemm)..
  Matrix<double> B = gemm(A,X);
  std::cout << "B = A*X = " << B << std::endl;
  // .. form solver for A ..
  LUSolver<double> solver(A);
  // .. solve A*X = B ..
  solver.solve(B);
  // .. check result
  std::cout << "A\\B = " << B << std::endl;  
  std::cout << "frobenius(A\\B-X) = " << frobenius(B-X) << std::endl;
  }

(Click here to view the output of this program)

A = size 3 by 3 Matrix of double:
[ 1 4 9 ]
[ 3 6 7 ]
[ 5 8 2 ]

X = size 3 by 2 Matrix of double:
[ 10 13 ]
[ 11 14 ]
[ 12 15 ]

B = A*X = size 3 by 2 Matrix of double:
[ 162 204 ]
[ 180 228 ]
[ 162 207 ]

A\B = size 3 by 2 Matrix of double:
[ 10 13 ]
[ 11 14 ]
[ 12 15 ]

frobenius(A\B-X) = 2.7231e-14

On toy examples the benefits of LUSolver vs getrf_inplace() are admittedly debatable, but it could in principle save you few moments thinking about (i) "what was the sequence of trsm() calls, again?" and (ii) "how do I undo the row interchanges, again?"

The following Solver classes are available:

• RCholeskySolver, for real symmetric positive definite A, within dense/RCholeskySolver.h
• ZCholeskySolver, for complex hermitian positive definite A, within dense/ZCholeskySolver.h
• RLDLTSolver, for real symmetric indefinite A, within dense/RLDLTSolver.h
• ZLDLHSolver, for complex hermitian indefinite A, within dense/ZLDLHSolver.h
• ZLDLTSolver, for complex symmetric A, within dense/ZLDLTSolver.h
• LUSolver, for general square A, within dense/LUSolver.h
• QRSolver, for full-rank least-squares solution ("tall and skinny" A), within dense/QRSolver.h

The real motivation for introducing the Solver concept now is that it sets the stage for all the sparse /multifrontal solvers to come (where there is no universally accepted "LAPACK/BLAS"-like fixed function layer). Each /multifrontal solver encapsulates related functionality (reordering, factoring, solving, refinement, and more) into a single Solver entity.

Appendix: Working with external data

To expand on an idea mentioned earlier, the main purpose of the "iterator" concept in the C++ STL is to decouple std::algorithms from the std::containers upon which they operate. All of the "Range" classes in this library are intended to serve a similar role. They allow algorithms (like gemm() or gemv()) to operate upon data other than our own Matrix and Vector containers, provided that the external data is "layout compatible" with the assumptions of the Range (column major, fixed leading dimension between columns). In practice, this is almost always the case, because it's really BLAS/LAPACK/Fortran interoperability that is sets these layout conventions for everybody. MyraMath just complies with the conventions, and provides MatrixRange to help you get your own compliant data into MyraMath algorithms.

This appendix example (tutorial/dense/external.cpp) shows how to use MatrixRange to feed "external" data into MyraMath's gemm_inplace() wrapper. Here, a C-style array of double[]'s serves as mock "external" data. To instantiate a MatrixRange, you must provide essentially the same parameters that you'd pass into the BLAS anyway: (i) a pointer to a buffer of numbers, (ii) size metadata, (iii) "leading dimension" metadata.

#include <myramath/dense/Matrix.h>
#include <myramath/dense/gemm.h>

#include <tests/myratest.h>

#include <iostream>

using namespace myra;

ADD_TEST("dense_external","[dense]") // int main()
  {
  // Make a random Matrix A
  auto A = Matrix<double>::fill_rmajor(3,3,{1,3,5,4,6,8,9,7,2}); 
  // Someone just gave us some input values B, intended
  // to represent a 3x3 matrix stored in column-major order.
  double B_values[9] = {0.1, 1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9};  
  // We also have some output values C, it's all zero.
  double C_values[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};  
  // Wrap up B_values and C_values into MatrixRange's
  MatrixRange<double> B(B_values,3,3,3); // (pointer,I,J,LD)
  MatrixRange<double> C(C_values,3,3,3); // (pointer,I,J,LD)
  // Examine inputs.
  std::cout << "A = " << A << std::endl;
  std::cout << "B = " << B << std::endl;
  std::cout << "C = " << C << std::endl;
  // Multiply C = A*B
  gemm_inplace(C,A,'N',B,'N');
  // Examine output.
  std::cout << "C = gemm(A,B) = " << C << std::endl;
  // Examine C_values, verify that the MatrixRange "writes-through" to it.
  for (int c = 0; c < 9; ++c)
    std::cout << "C_values[" << c << "] = " << C_values[c] << std::endl;
  }

(Click here to view the output of this program)

A = size 3 by 3 Matrix of double:
[ 1 3 5 ]
[ 4 6 8 ]
[ 9 7 2 ]

B = size 3 by 3 Matrix of double:
[ 0.1 3.4 6.7 ]
[ 1.2 4.5 7.8 ]
[ 2.3 5.6 8.9 ]

C = size 3 by 3 Matrix of double:
[ 0 0 0 ]
[ 0 0 0 ]
[ 0 0 0 ]

C = gemm(A,B) = size 3 by 3 Matrix of double:
[ 15.2 44.9 74.6 ]
[ 26 85.4 144.8 ]
[ 13.9 73.3 132.7 ]

C_values[0] = 15.2
C_values[1] = 26
C_values[2] = 13.9
C_values[3] = 44.9
C_values[4] = 85.4
C_values[5] = 73.3
C_values[6] = 74.6
C_values[7] = 144.8
C_values[8] = 132.7

If there was no intermediary like MatrixRange, and gemm_inplace() only accepted concrete Matrix objects, this wouldn't have been possible. At least, not without copying to/from temporaries for B and C. Similar idioms are used for getting "external" data into de-facto-standard formats for sparse matrices (like compressed-sparse-column) into MyraMath's /sparse routines and /multifrontal solvers. Furthermore, the /dense MatrixRange's and VectorRange's still have roles to play even in /sparse-land, because they serve as the interface types for right-hand-sides and solutions.

Appendix: A note about const-correctness.

As mentioned earlier, a MatrixRange is a view into data "owned" by some other container (a slice/window of a Matrix for example, or possibly external data). MatrixRange is a "mutable" view which permits modification of the underlying data. For "immutable" (const) views, there is a complementary class called CMatrixRange. This design is similar to the containers from the C++ STL, which present both (mutable) iterator's and (immutable) const_iterator's. Typically both range types are passed around as const&, when reading signatures just remember that the mutability of the underlying data is baked into the type itself (MatrixRange vs. CMatrixRange), not the const& qualifier. As an example, consider the signature of gemm_inplace(), which updates C+=A*B:

void gemm_inplace(
  const  MatrixRange<float>& C,  // C's data is mutable
  const CMatrixRange<float>& A,  // A's data is const
  const CMatrixRange<float>& B); // B's data is const

Because C is modified, it must be a (mutable) MatrixRange. But since A and B are passed as (immutable) CMatrixRange's, the routine is not permitted to modify them. Note that a MatrixRange can always be implicitly converted into a CMatrixRange, because it's always safe to "add constness" to a variable. The reverse conversion, "casting away const" from CMatrixRange to MatrixRange, is stylistically discouraged but can be accomplished by explicitly calling the .remove_const() member function of CMatrixRange.

Continue to Tutorial for /sparse routines, or go back to API Tutorials