Linear Algebra

Kokkos provides the Kokkos::View class to store multi-dimensional arrays. Kokkos itself does not provide tools BLAS (e.g., Mat-Mat multiplication) or LAPACK (e.g., LU factorizations and solves) operations. MParT provides a small collection of c++ classes and functions for easily working with one- and two-dimensional Kokkos::View objects as if they were vectors and matrices, respectively.

Below is a comparison of the MParT linear algebra tools and Eigen. These examples focus on views defined in the Kokkos::HostSpace memory space, but these operators are defined for views in device memory (e.g., mpart::DeviceSpace) as well.

The eigen code snippets require #include <Eigen/Core> and #include <Eigen/Dense> while the MParT snippets require #include <Kokkos_Core.hpp>, #include "MParT/Utilities/LinearAlgebra.h", and the addition of using namespace mpart;.

Construction

For completeness, we first show how to construct two matrices \(A\) and \(B\) using either Eigen or Kokkos.

Matrix Construction (Eigen)

Eigen::MatrixXd A(2,2);
Eigen::MatrixXd B(2,2);

A << 1.0, 2.0,
     0.5, 1.0;

B << 1.0,   2.0,
     0.125, 0.5;

Matrix Construction (Kokkos)

Kokkos::View<double**, Kokkos::HostSpace> A("A", 2, 2);
Kokkos::View<double** ,Kokkos::HostSpace> B("B", 2, 2);

A(0,0) = 1.0; A(0,1) = 2.0;
A(1,0) = 0.5; A(1,1) = 1.0;

B(0,0) = 1.0;   B(0,1) = 2.0;
B(1,0) = 0.125; B(1,1) = 0.5;

Summation

We will now compute \(C=2A+B\) by first computing \(A+B\) and then adding an additional \(A\) to the result. This is simply to demonstrate the use of both the + and += operators.

Matrix Sum (Eigen)

Eigen::MatrixXd C;
C = A + B;
C += A;

Matrix Sum (Kokkos+MParT)

Kokkos::View<double** ,Kokkos::HostSpace> C;
C = A + B;
C += A;

Multiplication

The * operator is overloaded in both Eigen and MParT to implement matrix multiplication. Currently MParT does not provide any built-in tools for componentwise multiplication. The snippet below computes the product \(AB\).

Matrix Product (Eigen)

Eigen::MatrixXd C;
C = A*B;

Matrix Product (Kokkos+MParT)

StridedMatrix<double, Kokkos::HostSpace> C;
C = A*B;

Computing products with transposes is also straightforward. Below we compute both \(A^TB\) and \(A^TB^T\). Because MParT must use out-of-class implementations, the transpose syntax is slightly different between Eigen and MParT, but equally straightforward:

Transposed Matrix Product (Eigen)

Eigen::MatrixXd C;
C = A.transpose() * B; // A^T B
C = A.transpose() * B.transpose() // A^T B^T

Transposed Matrix Product (Kokkos+MParT)

StridedMatrix<double, Kokkos::HostSpace> C;
C = transpose(A) * B; // A^T B
C = tranpose(A) * transpose(B); // A^T B^T

LU Factorization and Solve

MParT provides a class similar to Eigen’s PartialPivLU class for LU factorizations. The Kokkos::HostSpace version of MParT’s implementation is actually just a thin wrapper around the Eigen implementation. The DeviceSpace version uses the cuSolver library.

The following snippets compute \(C=A^{-1} B\) using an LU factorization of the matrix \(A\). The determinant of \(A\) is also computed from the factorization.

LU Solve (Eigen)

Eigen::PartialPivLU<Eigen::MatrixXd> solver;
solver.compute(A); // Compute the LU factorization of A

Eigen::MatrixXd C;
C = solver.solve(B);

double det = solver.determinant();

LU Solve (Kokkos+MParT)

mpart::PartialPivLU solver;
solver.compute(A); // Compute the LU factorization of A

Kokkos::View<double**, Kokkos::HostSpace> C;
C = solver.solve(B);

double det = solver.determinant();

The MParT mpart::PartialPivLU class also has a solveInPlace function that computes \(A^{-1}B\) in place by overwriting the matrix \(B\). The use of solveInPlace can reduce memory allocations because additional space to store the matrix C is not needed. Here is an example:

In-place LU Solve (Eigen)

Eigen::PartialPivLU<Eigen::MatrixXd> solver;
solver.compute(A);

B = solver.solve(B);

In-place LU Solve (Kokkos+MParT)

mpart::PartialPivLU solver;
solver.compute(A);
solver.solveInPlace(B);

Classes and Functions

template<typename MemorySpace>
class PartialPivLU

Mimics the interface of the Eigen::PartialPivLU class, but using Kokkos::Views and CUBLAS/CUSOLVER linear algebra.

Note that the layout of the matrices used in this class is important. Cublas expects column major (layout left).

Public Functions

inline PartialPivLU()

inline PartialPivLU(Kokkos::View<const double**, Kokkos::LayoutLeft, MemorySpace> A)

void compute(Kokkos::View<const double**, Kokkos::LayoutLeft, MemorySpace> A)

Computes the LU factorization of a matrix A

void solveInPlace(Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace> x)

Computes A^{-1}x and stores the results in x.

Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace> solve(StridedMatrix<const double, MemorySpace> x)

Returns a view containing A^{-1}x.

double determinant() const

Returns the determinant of the matrix A based on its LU factorization.

template<typename MemorySpace>
struct TransposeObject

Wrapper around a Kokkos view also storing a bool indicating whether the view should be transposed in matrix multiplications.

Template Parameters:

MemorySpace –

Public Functions

inline TransposeObject(StridedMatrix<const double, MemorySpace> viewIn, bool isTransposedIn = false)

inline int rows() const

inline int cols() const

Public Members

bool isTransposed

StridedMatrix<const double, MemorySpace> view

template<typename ...Traits>
TransposeObject<typename Kokkos::View<Traits...>::memory_space> mpart::transpose(Kokkos::View<Traits...> view)

Returns a transpose object around a view.

template<typename MemorySpace>
void mpart::dgemm(double alpha, TransposeObject<MemorySpace> A, TransposeObject<MemorySpace> B, double beta, StridedMatrix<double, MemorySpace> C)

BLAS-like function computing \( C = \alpha AB + \beta C\), \(C=\alpha A^TB + \beta C\), \(C=\alpha AB^T+\beta C\) or \(C=\alpha A^TB^T + \beta C\).

The matrix C must be preallocated.

Warning

doxygenfunction: Unable to resolve function “mpart::AddInPlace” with arguments (Kokkos::View<double**, Traits1…>, Kokkos::View<const double**, Traits2…>) in doxygen xml output for project “mpart” from directory: /home/runner/work/MParT/MParT/build/docs/doxygen/xml. Potential matches:

- template<typename ViewType1, typename ViewType2> void AddInPlace(ViewType1 x, ViewType2 y)

Warning

doxygenfunction: Unable to resolve function “mpart::AddInPlace” with arguments (Kokkos::View<double*, Traits1…>, Kokkos::View<const double*, Traits2…>) in doxygen xml output for project “mpart” from directory: /home/runner/work/MParT/MParT/build/docs/doxygen/xml. Potential matches:

- template<typename ViewType1, typename ViewType2> void AddInPlace(ViewType1 x, ViewType2 y)

template<typename ScalarType, typename ...Traits1, typename ...Traits2>
Kokkos::View<ScalarType**, typename Kokkos::View<ScalarType**, Traits1...>::memory_space> mpart::operator+(Kokkos::View<const ScalarType**, Traits1...> x, Kokkos::View<const ScalarType**, Traits2...> y)

template<typename ScalarType, typename ...Traits1, typename ...Traits2>
Kokkos::View<ScalarType*, typename Kokkos::View<ScalarType*, Traits1...>::memory_space> mpart::operator+(Kokkos::View<const ScalarType*, Traits1...> x, Kokkos::View<const ScalarType*, Traits2...> y)

z = x + y

Warning

doxygenfunction: Unable to resolve function “mpart::operator+=” with arguments (Kokkos::View<ScalarType**, Traits1…>, Kokkos::View<const ScalarType**, Traits2…>) in doxygen xml output for project “mpart” from directory: /home/runner/work/MParT/MParT/build/docs/doxygen/xml. Potential matches:

- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType**, Traits1...> &operator+=(Kokkos::View<ScalarType**, Traits1...> &x, Kokkos::View<ScalarType**, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType**, typename Kokkos::View<ScalarType**, Traits1...>::memory_space> &operator+=(Kokkos::View<ScalarType**, Traits1...> &x, Kokkos::View<const ScalarType**, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType*, Traits1...> &operator+=(Kokkos::View<ScalarType*, Traits1...> &x, Kokkos::View<ScalarType*, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType*, typename Kokkos::View<ScalarType*, Traits1...>::memory_space> &operator+=(Kokkos::View<ScalarType*, Traits1...> &x, Kokkos::View<const ScalarType*, Traits2...> const &y)

Warning

doxygenfunction: Unable to resolve function “mpart::operator+=” with arguments (Kokkos::View<ScalarType*, Traits1…>, Kokkos::View<const ScalarType*, Traits2…>) in doxygen xml output for project “mpart” from directory: /home/runner/work/MParT/MParT/build/docs/doxygen/xml. Potential matches:

- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType**, Traits1...> &operator+=(Kokkos::View<ScalarType**, Traits1...> &x, Kokkos::View<ScalarType**, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType**, typename Kokkos::View<ScalarType**, Traits1...>::memory_space> &operator+=(Kokkos::View<ScalarType**, Traits1...> &x, Kokkos::View<const ScalarType**, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType*, Traits1...> &operator+=(Kokkos::View<ScalarType*, Traits1...> &x, Kokkos::View<ScalarType*, Traits2...> const &y)
- template<typename ScalarType, typename ...Traits1, typename ...Traits2> Kokkos::View<ScalarType*, typename Kokkos::View<ScalarType*, Traits1...>::memory_space> &operator+=(Kokkos::View<ScalarType*, Traits1...> &x, Kokkos::View<const ScalarType*, Traits2...> const &y)

template<typename MemorySpace>
StridedMatrix<double, MemorySpace> mpart::operator*(TransposeObject<MemorySpace> A, TransposeObject<MemorySpace> B)

Returns C = A*B, or A^T*B or A^T*B^T or A*B^T