eigND - n-dimensional eigenvalues

Eigenvalues of a real non-symmetric N-dimensional matrix

eigND.c eigenvalues of a real non-symmetric N-dimensional matrix
E=eigND(A);
Calculates the eigenvalues of an N-dimensional 'square' matrix A.
if A is size [n x n x p x q x r x ...], then E are the eigenvalues of
each [n x n] square matrix from the first two dimensions of A.
E is size [n x 1 x p x q x r x ...]
Note that the second dimension of E always has size 1.

A must be real and is assumed to be non-symmetric. If A is complex then an error occurs, but if A is symmetric the computation is performed and the results may be less accurate than if symmetry was assumed. The lapack routine DGEEV is used.

Example:
A=[-100 -1e4;1 0];A(:,:,2)=A*2;
E=eigND(A)
E(:,:,1) =
-50.0000 +86.6025i
-50.0000 -86.6025i
E(:,:,2) =
1.0e+002 *
-1.0000 + 1.7321i
-1.0000 - 1.7321i
A is size [2 2 2] and the eigenvalues E are size [2 1 2].

Eigenvectors can be calculated by specifying a second output argement:
[V,E]=eigND(A);
This needs considerable memory space for V, and is barely faster than the built-in eig command. To save memory, E is NOT returned as a diagonal 'matrix' D, in contrast to the standard matlab eig command.

An optional 2nd input argument can be used. If this value is 1, then the input matrix A is OVERWRITTEN by lapack. This saves memory because the matrix A does not need to be duplicated prior to the lapack call.
This is only an issue when the input variable A is an actual variable rather than an evaluated expression. For example, eigND(A+0,1) does not modify A in the caller workspace, whereas eigND(A,1) does.

Tested on WIN32 and GLNXA64 platforms. Depending on the particular problem, the method seems to be 7-20% faster than using the standard matlab command EIG with nested loops for each of the additional dimensions of the data. The speed increase is presumably because there is no symmetry test and a number of memory allocations etc can be avoided. This assumes that the n-dimensional problem fits into memory.