linalg.h

Note

These files are actually split among multiple header files, which are all included by linalg.h

Matrix Multiplication

enum slap_ErrorCode slap_MatMulAdd(Matrix C, Matrix A, Matrix B, sfloat alpha, sfloat beta)

Perform in-place Matrix multiplication with addition.

\[ C = \beta C + \alpha A B \]

The two input matrices can be aliased, but neither can be aliased with the output, doing so will result in undefined behavior.

Examples

Normal matrix multiplication:

slap_MatMulAdd(C, A, B, 1, 0);

With transpose:

slap_MatMulAdd(C, slap_Transpose(A), B, 1, 0);

Add result to the output, with aliasing

slap_MatMullAdd(C, A, slap_Transpose(A), 1, 1);

Transposed output and input scaling

slap_MatMulAdd(slap_Transpose(C), A, B, 0.5, 1);

See also: slap_MatMulAB(), slap_MatMulAtB()

Header File: slap/linalg.h

Parameters:

C – Destination matrix
A – [in] Left input matrix
B – [in] Right input matrix
alpha – [in] Scaling factor on output
beta – [in] Scaling factor on input

enum slap_ErrorCode slap_MatMulAB(Matrix C, Matrix A, Matrix B)

Simple in-place matrix multiplication.

Calculates

\[ C = A B \]

This method performs basic matrix multiplication.

This method ignores all metadata besides the size (including transpose state, strides, etc.) so should be used with care.

It’s provided mainly to give the lowest possible cost for matrix multiplication, by directly indexing into the underlying data without checking for transposes or striding.

See also: slap_MatMulAdd(), slap_MatMulAtB()

Header File: slap/linalg.h

Parameters:

C – [out] Destination matrix
A – [in] Left input matrix
B – [in] Right input matrix

enum slap_ErrorCode slap_MatMulAtB(Matrix C, Matrix A, Matrix B)

Simple matrix-transpose matrix multiplication.

Calculates

\[ C = A^T B \]

Similar to slap_MatMulAB(), this method ignores all metadata information. See that method for more details.

See also: slap_MatMulAdd(), slap_MatMulAB()

Header File: slap/linalg.h

Parameters:

C – [out] Destination matrix
A – [in] Left input matrix
B – [in] Right input matrix

enum slap_ErrorCode slap_MatMulAtB(Matrix C, Matrix A, Matrix B)

Simple matrix-transpose matrix multiplication.

Calculates

\[ C = A^T B \]

Similar to slap_MatMulAB(), this method ignores all metadata information. See that method for more details.

See also: slap_MatMulAdd(), slap_MatMulAB()

Header File: slap/linalg.h

Parameters:

C – [out] Destination matrix
A – [in] Left input matrix
B – [in] Right input matrix

enum slap_ErrorCode slap_LowerTriMulAdd(Matrix C, const Matrix L, const Matrix B, double alpha, double beta)

Multiply a matrix by a lower triangular matrix and add to another matrix.

Simply ignores entries in the upper triangle of L.

Parameters:

C – Output matrix (m x p)
L – Lower triangular matrix (m x n)
B – Input matrix (n x p)
alpha – Scaling factor for L*B
beta – Scaling factor for C

enum slap_ErrorCode slap_UpperTriMulAdd(Matrix C, const Matrix U, const Matrix B, double alpha, double beta)

Multiply a matrix by an upper triangular matrix and add to another matrix.

Simply ignores entries in the lower triangle of U.

Parameters:

C – Output matrix (m x p)
U – Upper triangular matrix (m x n)
B – Input matrix (n x p)
alpha – Scaling factor for U*B
beta – Scaling factor for C

Triangular Matrices

enum slap_ErrorCode slap_MakeLowerTri(Matrix A)

Make a matrix lower triangular.

Sets all of the entries above the diagonal to zero.

Parameters:: A – Matrix to make lower triangular. Can be square or rectangular.

enum slap_ErrorCode slap_MakeUpperTri(Matrix A)

Make a matrix upper triangular.

Sets all of the entries below the diagonal to zero.

Parameters:: A – Matrix to make upper triangular. Can be square or rectangular.

bool slap_CheckLowerTri(const Matrix A): Check if a matrix is lower triangular.

bool slap_CheckUpperTri(const Matrix A): Check if a matrix is upper triangular.

Cholesky

enum slap_ErrorCode slap_Cholesky(Matrix A)

Perform a Cholesky decomposition.

Performs a Cholesky decomposition on the square matrix A, storing the result in the lower triangular portion of A.

Header File: slap/linalg.h

Parameters:: A – a square symmetric matrix
Returns:: slap error code. There is a dedicated code SLAP_CHOLESKY_FAIL if the factorization fails due to a negative value on the diagonal (matrix isn’t positive definite).

enum slap_ErrorCode slap_TriSolve(Matrix L, Matrix b)

Solve a linear system of equation for a triangular matrix.

Uses back-substitution to solve a system of equations of the following form:

\[ L x = b \]

for a lower-triangular matrix $ L $, or

\[ L^T x = b \]

if slap_IsTransposed(L) is true, or

\[ R x = b \]

if slap_GetType(L) == slap_TRIANGULAR_UPPER.

Header File: slap/linalg.h

Parameters:

L – [in] A triangular matrix. Assumed to be lower-triangular unless its slap_MatrixType is slap_TRIANGULAR_UPPER
b – [inout] The right-hand-side vector. Stores the solution upon completion.

Returns:

slap error code

enum slap_ErrorCode slap_CholeskySolve(Matrix A, Matrix b)

Solve a linear system of equation with a precomputed Cholesky decomposition.

Here’s a simple derivation of how the decomposition works:

\[\begin{split}\begin{align} &Ax = b \\ &L L^T x = b \\ &\implies y = L^{-1} b \;\text{ solved using back-substitution} \\ &L^T x = y \\ &\implies x = L^{-T} y \;\text{ solved using back-substitution} \\ \end{align}\end{split}\]

Example

enum slap_ErrorCode err;
err = slap_Cholesky(A);
if (err == SLAP_CHOLESKY_FAIL) {
  printf("Matrix not positive definite!\n);
}
slap_CholeskySolve(A, x);

Header File: slap/linalg.h

Parameters:

A – [in] A square matrix whose Cholesky decomposition is stored in the lower triangular portion of the matrix
b – [inout] The right-hand-side vector. Stores the solution upon completion of the function.

Returns:

slap error code

QR

Functions

enum slap_ErrorCode slap_QR(Matrix A, Matrix betas, Matrix temp)

Performs QR decomposition.

Calculates the QR decomposition using Householder reflections. The R factor is stored in the upper triangular portion of A, while the reflection vectors are stored below the diagonal, which can be used to recover the Q matrix.

See also: slap_ComputeQ(), slap_Qtb(), slap_LeastSquares()

Header File: slap/qr.h

Parameters:

A – A square or skinny matrix
betas – [out] Stores scaling factors needed to recover Q. Must have the same number of rows as A.
temp – [in] A temporary vector with the same number of rows as A.

enum slap_ErrorCode slap_ComputeQ(Matrix Q, const Matrix R, const Matrix betas, Matrix Q_work)

Computes the Q matrix from a previously-computed QR decomposition.

See also: slap_QR(), slap_LeastSquares(), slap_Qtb()

Header File: slap/qr.h

Parameters:

A – [out] A square matrix with the same number of rows as R.
R – [in] A square or skinny matrix containing the results from a QR decomposition.
betas – [in] A vector of scaling factors needed to recover the Q matrix.
Q_work – A temporary matrix of the same size as Q for intermediate calculations

enum slap_ErrorCode slap_Qtb(const Matrix R, const Matrix betas, Matrix b)

Calculates $Q^T b$ where $Q$ is the orthogonal matrix from a QR decomposition.

This method is useful for solving linear systems with the QR decomposition, where the right side needs to be multiplied by $Q^T$.

The result is stored in b.

See also: slap_QR(), slap_LeastSquares()

Header File: slap/qr.h

Parameters:

R – [in]
betas – [in] Scaling factors for reflections
b – [inout] Vector to multiply in-place

enum slap_ErrorCode slap_LeastSquares(Matrix A, Matrix b, Matrix betas, Matrix temp)

Solve a least squares problem.

Solves $ \text{minimize} || A x - b || $, where $A$ has more rows ( $m$) than columns ( $n$).

Both A and b are over-written by this method. A is overwritten with the it’s QR decomposition, and the first $n$ rows of b contain the solution.