Linear Algebra

Defines
#define	mat_IDX(m, i, j)   array_INDEX2( m, double, i, j )
	index matrix at i,j.
#define	matrix_CHECK(flag, m)
	is m a matrix?
#define	matrix_CHECKSQR(flag, m)
	is m a square matrix?
#define	vec_IDX(v, i)   array_INDEX1( v, double, i )
	index vector at i.
#define	vector_CHECK(flag, v)
	is v a vector?
Functions
Array *	matrix_get_col (Array *m, int col)
	get a column of a matrix as vector.
Array *	matrix_get_row (Array *m, int row, bool alloc)
	get a row of a matrix as vector.
Array *	matrix_mean (Array *a, int dim)
	calculate the mean of matrix a along dimension dim.
Array *	matrix_mult (const Array *m1, const Array *m2)
	Matrix Multiplication.
Array *	matrix_pca (Array *X, Array **var, bool alloc)
	Principal Components analysis.
gsl_matrix *	matrix_to_gsl (Array *in, bool alloc)
	convert 2D- Array struct to gsl_matrix.
Array *	matrix_transpose (Array *m, bool alloc)
	Matrix transpose.

---

Detailed Description

STATUS: work in progress Linear Algebra.

matrices Vectors and Matrices

1D and 2D Array structs of DType DOUBLE are used to represent vectors and matrices, respectively. I.e. that each function starting with matrix_*( Array *m ) accepts only 2D Arrays of DType DOUBLE functions starting with vector_*( Array *m ) accept only 1D Arrays of DType DOUBLE. This is checked with the Macros matrix_CHECK() and vector_CHECK().

---

Define Documentation

#define mat_IDX	(		m,
			i,
			j
	)			array_INDEX2( m, double, i, j )

index matrix at i,j.

Parameters:

	m	matrix
	i,j	indices

Definition at line 89 of file linalg.h.

#define matrix_CHECK	(		flag,
			m
	)

Value:

if( (m)==NULL ){                                                                    \
     errprintf("got NULL as matrix\n");                                         \
     flag=FALSE;                                                                        \
  } else if(!( (m)->ndim==2 && (m)->dtype==DOUBLE )){                       \
     char *dts="";                                                                      \
     array_DTYPESTRING( dts, m->dtype );                                        \
     errprintf("not a matrix, ndim=%i, dtype=%s\n", m->ndim, dts ); \
     flag=FALSE;                                                                        \
  } else { flag=TRUE; }

is m a matrix?

Usage:

     bool ismatrix;
     matrix_CHECK( ismatrix, X );
     if( !ismatrix ) return NULL;

Parameters:

	flag	(output) (bool) set by macro
	m	(input) Array* to check

Definition at line 49 of file linalg.h.

#define matrix_CHECKSQR	(		flag,
			m
	)

Value:

if(!( (m)->ndim==2 && (m)->dtype==DOUBLE )){                                    \
     char *dts="";                                                                          \
     array_DTYPESTRING( dts, (m)->dtype );                                          \
     errprintf("not a matrix, ndim=%i, dtype=%s\n", (m)->ndim, dts );       \
     flag=FALSE;                                                                            \
  } else if( (m)->size[0]!=(m)->size[1] ) {                                     \
     errprintf("Matrix is not square, got (%i,%i)\n",(m)->size[0],(m)->size[1]); \
     flag=FALSE;                                                                            \
  } else { flag=TRUE; }

is m a square matrix?

Usage:

     bool ismatrix;
     matrix_CHECKSQR( ismatrix, X );
     if( !ismatrix ) return NULL;

Parameters:

	flag	(output) (bool) set by macro
	m	(input) Array* to check

Definition at line 72 of file linalg.h.

#define vec_IDX	(		v,
			i
	)			array_INDEX1( v, double, i )

index vector at i.

Parameters:

	v	vector
	i	index

Definition at line 98 of file linalg.h.

#define vector_CHECK	(		flag,
			v
	)

Value:

if( (v)==NULL ){                                                                        \
     errprintf("got NULL as vector\n");                                             \
     flag=FALSE;                                                                            \
  } else if(!( (v)->ndim==1 && (v)->dtype==DOUBLE ) ){                      \
     char *dts="";                                                                          \
     array_DTYPESTRING( dts, v->dtype );                                            \
     errprintf("not a vector, ndim=%i, dtype=%s\n", v->ndim, dts );     \
     flag=FALSE;                                                                            \
  } else { flag=TRUE; }

is v a vector?

Usage:

     bool isvector;
     vector_CHECK( isvector, X );
     if( !isvector ) return NULL;

Parameters:

	flag	(output) (bool) set by macro
	v	(input) Array* to check

Definition at line 113 of file linalg.h.

---

Function Documentation

Array* matrix_get_col	(	Array *	m,
		int	col
	)

get a column of a matrix as vector.

Note:

memory copy is necessary, because the memory of a column is not aligned. It is therefore slower than calling matrix_get_row() with alloc=FALSE. If you need a lot of columns, you might want consider calling matrix_transpose() first and matrix_get_row() afterwards.

Parameters:

	m	the matrix
	col	index to the column

Returns:

a column vector

Definition at line 175 of file linalg.c.

Array* matrix_get_row	(	Array *	m,
		int	row,
		bool	alloc
	)

get a row of a matrix as vector.

Note:

if alloc is TRUE, the returned vector is independant. If alloc is FALSE, the memory is shared between m and the row vector but is owned by the matrix (an array_free(vec) call does not free the memory, but array_free(m) does).

Parameters:

	m	the matrix
	row	index to the row
	alloc	allocate memory for the row-data, or only pass the pointer to the new Array

Definition at line 138 of file linalg.c.

Array* matrix_mean	(	Array *	a,
		int	dim
	)

calculate the mean of matrix a along dimension dim.

Parameters:

	a	the matrix
	dim	along which dimension?

Returns:

vector of size a->size[1-dim]

Definition at line 99 of file linalg.c.

Array* matrix_mult	(	const Array *	m1,
		const Array *	m2
	)

Matrix Multiplication.

m1 and m2 must be matrices such that m1->size[1]==m2->size[0].

Parameters:

m1,m2

the matrices

Returns:

new matrix of size m1->size[1], m2->size[0]

Definition at line 200 of file linalg.c.

Array* matrix_pca	(	Array *	X,
		Array **	var,
		bool	alloc
	)

Principal Components analysis.

This implementation uses SVD to calculate the PCA. Example:

     Array *X = get_data_from_somewhere();
     Array *var, *X_pca;
     matrix_pca( X, &var, FALSE ); // overwrite X
     X_pca = matrix_pca( X, &var, TRUE ); // do not touch X

Parameters:

	X	a 2D array observations x variables containing the data
	var	output: vector of eigenvalues of XX^T in decreasing order; if you pass NULL, it is ignored;
	alloc	if true, new memory is allocated and returned. Else X is overwritten.

Returns:

pointer to the PCA'd matrix

Definition at line 282 of file linalg.c.

gsl_matrix* matrix_to_gsl	(	Array *	in,
		bool	alloc
	)

convert 2D- Array struct to gsl_matrix.

This function works also for non-double Array's, but the input array MUST be of DType DOUBLE if alloc is FALSE. If alloc is TRUE, the type is casted to double during the copy-process.

Note:

the overhead for alloc=FALSE is very small, for alloc=TRUE MxN copy processes are necessary.

Warning:

use alloc=FALSE only if you have a double-Array

Note:

if you free with gsl_matrix_free, the memory is NOT* free'd but still belongs to the Array!

Example:

     gsl_matrix *g;
     Array *a = array_new_dummy( DOUBLE, 2, 10, 11 );

     g = matrix_to_gsl( a, FALSES );
     gsl_matrix_free( g ); //does not free a->data (!)

     array_free( a ); // does free a->data

Parameters:

	in	the input 2D-array
	alloc	if true, copy the memory pointed to by in; if false, simply put the pointer into the gsl_matrix struct (bad stuff happens if it is not DType DOUBLE);

Returns:

the gsl_matrix

Definition at line 59 of file linalg.c.

Array* matrix_transpose	(	Array *	m,
		bool	alloc
	)

Matrix transpose.

Todo:

implement a faster way

Parameters:

	m	matrix
	alloc	if TRUE, return a new matrix that is the transpose of m; else do in-place transpose

Definition at line 238 of file linalg.c.