# 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