High speed computations

In order to write fast matrix computations (to compute eigenvalues for example), you need to use LAPACK. As it is not that easy to run a simple example on mac os X, here’s a little howto with some explanations on what the parameters to the functions are.

First of all, the LAPACK (high level matrix computations, system of equation solving) and BLAS (basic matrix/vector operations) are part of the Accelerate framework. You do not need to download or install anything.

BLAS

Let’s start with a simple example. Find:

Let’s call the matrix on the right “T”. So we want to compute T’ (transpose of T) multiply T. The convention for the function to call is

• s  working with real numbers
• ge general matrices (not symmetric or Hermitian matrix)
• mm matrix multiplication

Our function is sgemm (see cblass quick ref). What are the parameters ? The prototype of this function is :

sgemm (Order, TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)

Where :

Order	Can be CblasRowMajor or CblasColMajor (if we iterate through the elements in memory, do we go through rows one by one or columns)
TRANSA	Can be CblasNoTrans (do not transpose before operation), CblasTrans (transpose before operation) or CblasConjTrans (conjugate transpose before operation)
TRANSB	Same as TRANSA but for the second matrix.
M	Number of rows of A (size of A is MxK).
N	Number of columns of B.(size of B is KxN).
K	Number of columns of A / rows of B.
ALPHA	Operation is ALPHA * AB + BETAC. This is 1.0 to have plain AB
A	Pointer to the first value of the matrix A. (const float *)
LDA	Leading dimension of A. Since we are in RowMajor, this is the number of columns of A = K.
B	Pointer to the first value of matrix B.
LDB	Leading dimension of B. In our case it’s the columns of B = N.
BETA	See ALPHA. Set to 0.0 to have AB.
C	Pointer to the first value of resulting matrix C.
LDC	Leading dimension of C. (N)

Note that M is the columns of A if TRANSA is CblasTrans. You mus view M,N and K as the sizes just before multiplication (once transpositions have been applied).

Our function call will be :

sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,2,2,3,1,T,2,T,2,0, A,2).

597 octets

matrix Simple matrix multiplication using BLAS.

This program is compiled and run with :

# gcc matrix.c -o matrix -framework Accelerate
# ./matrix

LAPACK

Let’s find the eigenvectors of this simple matrix:

2	-3
-3	6

The function to find this result is *ssyevr (SEP). This name is made of:

s	Single precision real value (not complex).
sy	Main matrix is symmetric.
ev	Eigenvector calculation.
r	Relatively Robust Representation version of this calculation.

The prototype for this function is (API):

ssyevr_( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO )

Yes, this is not a joke ! It is this long. The underscore at the end of the function was added by f2c. Let’s go through all these parameters one by one but before diving into these parameters, one must know that all matrices are in column major and there is no parameter to change this.

JOBZ	‘N’: compute eigenvalues only, ‘V’: eigenvalues and eigenvectors.
RANGE	‘A’: find all eigenvalues, ‘V’: find all eigenvalues in the interval ]VL,VU], ‘I’: find the eigenvalues with index in [IL,IU].
UPLO	‘U’: symmetric matrix is only stored in the upper triangle. ‘L’: for lower.
N	Order of the matrix A (size of A is NxN).
A	Source array.
LDA	Leading dimension of A (=N).
VL	See RANGE. (L stands for Low)
VU	See RANGE. (U stands for Up)
IL	See RANGE.
IU	See RANGE.
ABSTOL	Absolute error tolerance for the eigenvalues. If set to zero, machine precision will be used during convergence test. Higher values = faster but less precise.
M	Output of the number of eigenvalues found.
W	Output eigenvalues in ascending order.
Z	Output array. First M columns contain the orthonormal eigenvectors. i-th column contains the eigenvector associated with the i-th eigenvalue.
LDZ	Leading dimension of Z (=N).
ISUPPZ	Output array of integers telling which eigenvectors are nonzero. ??
WORK	Workspace (real array)
LWORK	Size of the workspace. Should be >= (NB+6)N where NB is the maximal blocksize for SSYTRD and SORMTR returned by ILAENV
IWORK	Workspace (integer array)
LIWORK	Size of IWORK array. Should be >= 10N
INFO	Result information. 0 = success, -i: i-th argument had an illegal value, > 0 internal error

To build our call, we first used ILAENV to find the maximal blocksize for “SSYTRD” and “SORMTR” (we have been told to do so by the documentation on SSYEVR. The result depends on platform implementation. The default would be 32.

function NAME, OPTS, N1, N2, N3, N4 ">ILAENV

We then pass all arguments to SSYEVR. Please note that passing the arguments directely will not work (you have to set them into variables and pass the values as references since this is what a Fortran routine wants).

JOBZ	‘V’	RANGE	‘A’	UPLO	‘U’
N	2	LDA	2	VL	0 (not used)
VU	0 (not used)	IL	0 (not used)	IU	0 (not used)
ABSTOL	0.01	M	&value_count	W	&eigenvalues
Z	&eigenvectors	LDZ	2	ISUPPZ	&isuppz
WORK	&work	LWORK	(32+6) *2 = 76	IWORK	&iwork
LIWORK	20	INFO	&info

Here is the sample code :

4 Kb	eigen Compute eigenvectors using LAPACK.

To compile and execute it, you need to do:

# gcc -o eigen -framework Accelerate eigen.c
# ./eigen