Linear Algebra

Files
file	linalg_custom.h
	Simple linear algebra routines.

Data Structures
struct	CholeskyInfo
	Stores info about a Cholesky decomposition. More...
struct	Matrix
	Represents a matrix of double-precision data. More...

Enumerations
enum	MatrixLinearAlgebraLibrary { libBLAS = 0, libMKL = 1, libEigen = 2, libInternal = 3 }
	List of supported linear algebra libraries. More...

Functions
CholeskyInfo	DefaultCholeskyInfo ()
	Construct a default CholeskyInfo object. More...
void	FreeFactorization (CholeskyInfo *cholinfo)
	Frees any data stored by the external library. More...
int	MatrixAddition (Matrix *A, Matrix *B, double alpha)
	Add two matrices of the same size, storing the result in B. More...
int	MatrixCholeskyFactorize (Matrix *mat)
	Compute the Cholesky decomposition on the matrix A. More...
int	MatrixCholeskyFactorizeWithInfo (Matrix *mat, CholeskyInfo *cholinfo)
	Compute the Cholesky decomposition on the matrix A. More...
int	MatrixCholeskySolve (Matrix *A, Matrix *b)
	Solve a linear system using a precomputed Cholesky factorization. More...
int	MatrixCholeskySolveWithInfo (Matrix *A, Matrix *b, CholeskyInfo *cholinfo)
	Solve a linear system using a precomputed Cholesky factorization. More...
void	MatrixMultiply (Matrix *A, Matrix *B, Matrix *C, bool tA, bool tB, double alpha, double beta)
	Matrix multiplication with scaling. More...
void	MatrixSymmetricMultiply (Matrix *Asym, Matrix *B, Matrix *C, double alpha, double beta)
	Matrix multiplication with a symmetric matrix A. More...
void	MatrixCopyDiagonal (Matrix *dest, Matrix *src)
	Copy just the diagonal element of src to the diagonal of dest. More...
enum MatrixLinearAlgebraLibrary	MatrixGetLinearAlgebraLibrary ()
	Get the linear algebra library currently being used. More...
void	MatrixPrintLinearAlgebraLibrary ()
	Prints which linear algebra library is being used to stdout. More...
void	MatrixLinAlgTimeStart ()
void	MatrixLinAlgTimeStop ()
void	MatrixLinAlgTimeReset ()
double	MatrixGetLinAlgTimeMilliseconds ()
Matrix	NewMatrix (int rows, int cols)
	Allocate a new matrix on the heap. More...
int	MatrixSetConst (Matrix *mat, double val)
	Sets all of the elements in a matrix to a single value. More...
int	FreeMatrix (Matrix *mat)
	Free the data for a matrix. More...
int	MatrixNumElements (const Matrix *mat)
	Get the number of elements in a matrix, i.e. m * n. More...
int	MatrixGetLinearIndex (const Matrix *mat, int row, int col)
	Get the linear index for a given row and column in the matrix. More...
double *	MatrixGetElementTranspose (const Matrix *mat, int row, int col, bool istranposed)
	Get the element of a matrix or its transpose. More...
int	MatrixSetElement (Matrix *mat, int row, int col, double val)
	The a matrix element to a given value. More...
double *	MatrixGetElement (const Matrix *mat, int row, int col)
	Get the element of a matrix given row, column indices. More...
int	MatrixCopyTranspose (Matrix *dest, Matrix *src)
	Copy a matrix to another matrix, transposed. More...
int	MatrixCopy (Matrix *dest, Matrix *src)
	Copy a matrix to another matrix. More...
int	MatrixScaleByConst (Matrix *mat, double alpha)
	Scale a matrix by a constant factor. More...
double	MatrixNormedDifference (Matrix *A, Matrix *B)
	Return the normed difference between 2 matrices of the same size. More...
int	MatrixFlatten (Matrix *mat)
	Flatten a 2D matrix to a column vector. More...
int	MatrixFlattenToRow (Matrix *mat)
	Flatten a 2D matrix to a row vector. More...
int	PrintMatrix (const Matrix *mat)
	Print the elements of a matrix to stdout. More...
int	PrintRowVector (const Matrix *mat)
	Print the entire matrix as a row vector. More...

Detailed Description

Enumeration Type Documentation

MatrixLinearAlgebraLibrary

enum MatrixLinearAlgebraLibrary

List of supported linear algebra libraries.

Function Documentation

DefaultCholeskyInfo()

CholeskyInfo DefaultCholeskyInfo

(

)

Construct a default CholeskyInfo object.

Returns

A CholeskyInfo with default values

FreeFactorization()

void FreeFactorization

(

CholeskyInfo *

cholinfo

)

Frees any data stored by the external library.

This should be called prior to creating a new factorization, which will allocate new memory for the factorization.

Parameters

cholinfo

FreeMatrix()

int FreeMatrix

(

Matrix *

mat

)

Free the data for a matrix.

Note this does NOT attempt to free the matrix object itself, only the data it wraps.

Parameters

mat

Postcondition

mat.data will be NULL.

Returns

0 if successful

MatrixAddition()

int MatrixAddition	(	Matrix *	A,
		Matrix *	B,
		double	alpha
	)

Add two matrices of the same size, storing the result in B.

Performs the following operation:

\[ B = B + \alpha A \]

Parameters

[in]	A	any matrix of size (m,n)
[in,out]	B	any matrix of size (m,n)
[in]	alpha	scalar factor on A

Returns

0 if successful

MatrixCholeskyFactorize()

int MatrixCholeskyFactorize

(

Matrix *

mat

)

Compute the Cholesky decomposition on the matrix A.

Not supported by all libraries. Prefer to use MatrixCholeskyFactorizeWithInfo().

Parameters

mat	A square, positive-definite matrix

Returns

0 if the call was successful

MatrixCholeskyFactorizeWithInfo()

int MatrixCholeskyFactorizeWithInfo	(	Matrix *	mat,
		CholeskyInfo *	cholinfo
	)

Compute the Cholesky decomposition on the matrix A.

Parameters

[in,out]	mat	A square, positive-definite matrix
	cholinfo	CholeskyInfo object for storing info about the factorization

Postcondition

cholinfo.success can be checked to see if the factorization was successful

Returns

0 if successful

MatrixCholeskySolve()

int MatrixCholeskySolve	(	Matrix *	A,
		Matrix *	b
	)

Solve a linear system using a precomputed Cholesky factorization.

Overwrite the input vector b. Prefer to use the more robust MatrixCholeskySolveWithInfo().

Parameters

[in]	A	A square matrix whose Cholesky decomposition has already been computed.
[in,out]	b	The right-hand-side vector. Stores the solution vector.

Returns

0 if successful

MatrixCholeskySolveWithInfo()

int MatrixCholeskySolveWithInfo	(	Matrix *	A,
		Matrix *	b,
		CholeskyInfo *	cholinfo
	)

Solve a linear system using a precomputed Cholesky factorization.

Parameters

[in]	A	A square matrix whose Cholesky decomposition has already been computed.
[in,out]	b	The right-hand-side vector. Stores the solution vector.
	cholinfo	Information about the precomputed Cholesky factorization in A.

Returns

MatrixCopy()

int MatrixCopy	(	Matrix *	dest,
		Matrix *	src
	)

Copy a matrix to another matrix.

Parameters

dest	a matrix of size (m,n)
src	a matrix of size (n,m)

Returns

0 if successful

MatrixCopyDiagonal()

void MatrixCopyDiagonal	(	Matrix *	dest,
		Matrix *	src
	)

Copy just the diagonal element of src to the diagonal of dest.

Parameters

dest	Destination matrix
src	Source matrix

MatrixCopyTranspose()

int MatrixCopyTranspose	(	Matrix *	dest,
		Matrix *	src
	)

Copy a matrix to another matrix, transposed.

Parameters

dest	a matrix of size (m,n)
src	a matrix of size (n,m)

Returns

0 if successful

MatrixFlatten()

int MatrixFlatten

(

Matrix *

mat

)

Flatten a 2D matrix to a column vector.

Changes the row and column data so that the matrix is now a column vector. The underlying data is unchanged.

Parameters

mat	Matrix to be flattened.

Returns

0 if successful

MatrixFlattenToRow()

int MatrixFlattenToRow

(

Matrix *

mat

)

Flatten a 2D matrix to a row vector.

Changes the row and column data so that the matrix is now a row vector. The underlying data is unchanged.

Parameters

mat	Matrix to be flattened

Returns

0 if successful

MatrixGetElement()

double* MatrixGetElement	(	const Matrix *	mat,
		int	row,
		int	col
	)

Get the element of a matrix given row, column indices.

Parameters

mat	Matrix of nonzero size
row	Row index
col	Column index

Returns

A pointer to the element of the matrix. NULL for invalid input.

MatrixGetElementTranspose()

double* MatrixGetElementTranspose	(	const Matrix *	mat,
		int	row,
		int	col,
		bool	istranposed
	)

Get the element of a matrix or its transpose.

If istransposed is false, then this method acts just like MatrixGetElement(). Otherwise, it is equalivalent to flipping the row and col arguments to MatrixGetElement().

Parameters

mat	Matrix with nonzero size and initialized data
row	Row index
col	Column index
istranposed	Are the indicies for the transpose of A?

Returns

Pointer to the data at the given element.

MatrixGetLinearAlgebraLibrary()

enum MatrixLinearAlgebraLibrary MatrixGetLinearAlgebraLibrary

(

)

Get the linear algebra library currently being used.

Its value is determined by the build system and cannot be changed at runtime.

Returns

The linear algebra library being used by the system.

MatrixGetLinearIndex()

int MatrixGetLinearIndex	(	const Matrix *	mat,
		int	row,
		int	col
	)

Get the linear index for a given row and column in the matrix.

Converts a cartesian index of row and column into a linear index for accessing an element of the underlying data.

Parameters

mat	Matrix with nonzero size and initialized data
row	Row index
col	Column index

Returns

Linear index corresponding to row and col. Returns -1 for a bad input.

MatrixMultiply()

void MatrixMultiply	(	Matrix *	A,
		Matrix *	B,
		Matrix *	C,
		bool	tA,
		bool	tB,
		double	alpha,
		double	beta
	)

Matrix multiplication with scaling.

Perform the computation

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

Parameters

[in]	A	Matrix of size (m,n)
[in]	B	Matrix of size (n,p)
[in,out]	C	Output matrix of size (m,p)
[in]	tA	Should A be transposed
[in]	tB	Should B be transposed
[in]	alpha	scalar on the \( A B \) term
[in]	beta	scalar on the \( C \) term. Set to zero for pure matrix multiplication.

MatrixNormedDifference()

double MatrixNormedDifference	(	Matrix *	A,
		Matrix *	B
	)

Return the normed difference between 2 matrices of the same size.

Returns \( \sqrt{\sum_{i=0}^{m-1} \sum_{j=0}^{n-1} (A_{ij} - B_{ij})^2 } \)

Parameters

A	A matrix of dimension (m,n)
B	A matrix of dimension (m,n)

Returns

MatrixNumElements()

int MatrixNumElements

(

const Matrix *

mat

)

Get the number of elements in a matrix, i.e. m * n.

Parameters

mat	Any matrix

Returns

Number of elements in the matrix

MatrixPrintLinearAlgebraLibrary()

void MatrixPrintLinearAlgebraLibrary

(

)

Prints which linear algebra library is being used to stdout.

MatrixScaleByConst()

int MatrixScaleByConst	(	Matrix *	mat,
		double	alpha
	)

Scale a matrix by a constant factor.

Parameters

mat	Fully initialized matrix of non-zero size. Values will be modified.
alpha	scalar by which to multiply the matrix

Returns

0 if successsful

MatrixSetConst()

int MatrixSetConst	(	Matrix *	mat,
		double	val
	)

Sets all of the elements in a matrix to a single value.

Parameters

mat	Matrix to be modified
val	Value to which each element will be set

Returns

0 if successful

MatrixSetElement()

int MatrixSetElement	(	Matrix *	mat,
		int	row,
		int	col,
		double	val
	)

The a matrix element to a given value.

Parameters

mat	Matrix with nonzero size and initialized data
row	Row index
col	Column index
val	Value to which the element should be set

Returns

0 if successful

MatrixSymmetricMultiply()

void MatrixSymmetricMultiply	(	Matrix *	Asym,
		Matrix *	B,
		Matrix *	C,
		double	alpha,
		double	beta
	)

Matrix multiplication with a symmetric matrix A.

Perform the following computation

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

For a symmetric matrix \( A \).

Parameters

[in]	Asym
[in]	B
[in,out]	C
[in]	alpha
[in]	beta

NewMatrix()

Matrix NewMatrix	(	int	rows,
		int	cols
	)

Allocate a new matrix on the heap.

Data will not be initialized. Wrapper around a call to malloc. Must be followed by a call to FreeMatrix.

Parameters

rows	number of rows in the matrix
cols	number of columns in the matrix

Returns

A new matrix

PrintMatrix()

int PrintMatrix

(

const Matrix *

mat

)

Print the elements of a matrix to stdout.

Precision of the printing can be controlled by the global variable PRECISION.

Parameters

mat	Matrix to be printed

Returns

0 if successful

PrintRowVector()

int PrintRowVector

(

const Matrix *

mat

)

Print the entire matrix as a row vector.

Same result as calling PrintMatrix() after a call to MatrixFlattenToRow().

Parameters

mat	Matrix to be printed

Returns

0 if successful