Jordan canonical form

[John Swensen]

A recent email about logm for 4x4 affine matrices led me to discover that the 
Octave's logm routine is flawed.  The Octave logm function will only work for 
matrices that are diagonalizable.  In general, the matrix log is only defined 
for non-singular matrices, but the method used in logm.m is only valid for 
nonsingular, diagonizable matrices.

<snippet from the current Octave logm.m>
[V, D] = eig (A);
B = V * diag (log (diag (D))) * inv (V);

The correct way to do the matrix log is to decompose the matrix into its Jordan 
blocks and take the matrix log of each Jordan block.  As the the Jordan block 
for an eigenvalue \lambda_i is  nilpotent of order k, where k dimension of the 
Jordan block, then the matrix log of a Jordan block can be written as an 
infinite series where only the first k-1 terms are non-zero.  I have a working 
version of this method, with the assumption that  I have a Jordan decomposition 
of a matrix.  

However, it appears that Octave doesn't have the jordan function implemented.  
Is this correct, or am I missing something?  I couldn't find it in either 
OctaveForge or the main octave sources.  Does someone already have an efficient 
method for finding the Jordan canonical form of a matrix?  I can implement it, 
but I am sure mine is going to stink as the only way I can think of doing it is:
1) Find the eigenvalues
2) Look at the null space of (\lambda_i * I - A)^k for k = [1:n]
3) Use the difference in rank between successive powers to determine Jordan 
structure and the vectors in the null space as the eigenvectors and generalized 
eigenvectors.
4) Repeat 2-3 for each eigenvalue of the matrix

Does anyone have a better method?  Once I have a working jordan function, I 
will submit a patch with the updated logm function.  I put a check in to the 
logm function for singular matrices and a check to use the diagonalization 
method if the matrix is diagonalizable.

John Swensen