Four MATLAB programs for the Jordan matrix decomposition
of complex floating point matrices:

1) jordan.m     top level program. uses Matlab Eigenvalue routine
2) jord.m       preferred starting point. Needs Eigenvalues.
3) nulld.m      robust null space sizer
4) jordlook.m   diagnostic program for tough null spaces



function [M,D] = jordan(a)

%Jordan  Compute the Jordan decomposition of square matrix A
%        such that A*M = M*D.  The elements of A may be real
%        or complex floating point numbers.
%
%        The eigenvalues of A are on the diagonal of D
%        and the generalized eigenvectors of A are the columns of M.
%
%        This routine uses the Jord function.
%
%usage: [M, D] = jordan(A)
%
%tested under 5.3.1
%
%see also: sym/jordan mhelp('jordan') eig hess schur

%Paul Godfrey
%pgodfrey@intersil.com
%2-05-01

small=2*sqrt(eps);

[r,c]=size(a);
if r~=c
   error('A matrix must be square!')
end
n=r;

if n==1
   D=a;
   M=1;
   return
end

if n<1
   D=[];
   M=[];
   return
end

[m,d]=eig(hess(a));
d=sort(diag(d));
tiny=norm(a)*eps;

%zero out zero evalues
p=find(abs(d)<=tiny);
if ~isempty(p)
   d(p)=0;
end

%A*M=M*D

[M,D]=jord(a,d,small);

return


function [M,D] = jord(a,d,small)
%Jord    Compute the Jordan decomposition of square matrix A
%        such that A*M = M*D.  The elements of A may be real
%        or complex floating point numbers.
%
%usage: [M, D] = jord(A, e, small)
%
%        The eigenvalues of A will appear on the diagonal of D
%        and the generalized eigenvectors of A are the columns of M.
%
%        The accuracy of this function is limited by the accuracy
%        in determining the exact eigenvalues of A.
%        For this reason this routine requires the user provide
%        the eigenvalues of A in the input vector, e.
%
%        This routine is used by the JORDAN function.
%        and uses the built in NULLD function.
%
%        Determining the Jordan form is a numerically ill conditioned
%        problem. If the condition number of A is large then you may
%        not get accurate results. Also,  if you have problems then
%        consider running the program again. A random combination of
%        Null Space vectors is used each time in an attempt to avoid
%        matrix singularity problems due to linearly dependent Null
%        Space subspace issues.
%        Also, consider changing the optional SMALL parameter.
%        Typically small ~= 2*sqrt(eps). If the program detects an error
%        then the program variables are saved in the JORD.MAT file for
%        inspection by the JordLook program attached below.
%
%
%tested under 5.3.1
%
%see also: sym/jordan eig hess schur

%Paul Godfrey
%pgodfrey@intersil.com
%Feb., 7, 2001

norma=sqrt(mean(mean(abs(a))));
tiny=norma*eps;

if nargin<3
   small=(1+norma)*sqrt(eps);
end

[r,c]=size(a);
if r~=c
   error('A matrix must be square!')
end
n=r;
I=eye(n);

if r==1
   D=a;
   M=1;
   return
end

if r<1
   D=[];
   M=[];
   return
end

condofa=cond(a);
if condofa>1e6
   Condition_number_of_A=condofa
   warning('Condition number of A is very large!');
end


d=d(:);
if size(d,1)~=n
   d=d
   error('Wrong number of eigenvalues provided!')
end

da=det(a);
dp=prod(d);
e=abs(abs(da)-abs(dp));
if e>sqrt(eps)
   disp(' ')
   warning('Eigenvalues are probably inaccurate!')
end

ds=flipud(sort(d));
sds=size(ds,1);

du=flipud(unique(ds));
sdu=size(du,1);

if sdu==sds
%  no repeated evalues at all
   [M,D]=eig(a);
   return
end

%  yes some repeated evalues
M=[];

for kk=1:sdu
       e=du(kk);
       ameig=sum(ismember(ds,e));% evalue alg mult
       a1=a-e*I;
    if ameig==1
       [u,s,v]=svd(a1);
       M=[M v(:,end)];
    else
       pp=0;
       ns=[];

           pp=pp+1;
           aa=I;

           for k=1:ameig
               aa=a1*aa;
               nn=size(nulld(aa,small),2);
               ns(k)=nn;
           end

       nsaa=[0; ns(:)]';
       dns=diff(nsaa);

       if max(dns)~=dns(1)
          Cond_of_A=cond(a)
          save jord
          M=I; D=I;
          error('Null space size error 1')
       end

%      determine the number and lengths of the eigen-chains
       clear ec;
       ec(1:dns(1))=1; % number of chains
       for k=2:length(dns)
           ec(1:dns(k))=ec(1:dns(k))+1;
       end

       if sum(ec)~=ameig
          Cond_of_A=cond(a)
          save jord
          M=I;D=I;
          error('Null space size error 2')
       end

%     ec(1) contains the length of the longest chain(s)
      k=1;
      
      clear jv;
      while k<= dns(1)
          p=find(ec==ec(k));
%         could be more than 1 chain of a given length
          if isempty(p)
             Cond_of_A=cond(a)
             save jord
             M=I; D=I;
             error('Null space size error 3');
          end
          aa=I;
          for m=1:ec(k)
             aa=aa*a1;
          end

%         [u,s,v]=svd(aa);
%         last cols of v are the nullspace

          pp=max(size(p));
          v=nulld(aa, small);         
          jv(:,p)=v*(rand(size(v,2),pp)-0.5)*16;
%         take a random(+-8) LC of the NSVs to make sure
%         that we use no repeats from v.--> Large cond(M)
%         This can happen since N1cN2cN3...  
          k=k+pp;
      end
      clear v;

      for k=1:dns(1) % do each chain
          v(:,1)=jv(:,k);% termination vector for kth chain
          for p=2:ec(k)
              v(:,p)=a1*v(:,p-1);
          end
          vv=fliplr(v(:,1:ec(k)));
          M=[M vv];
      end
   end
end

k=abs(det(M))^(-1/n);
M=k*M; % make det(m)=+-1
% if we did our job right then
% m is square and always non-singular
Mi=inv(M);
D=Mi*a*M;
d0=diag(D);
d1=diag(D,1);
D=diag(d0)+diag(d1,1);

return


function Z = nulld(A, small)
%NULLD   Null space.
%   Z = NULL(A) is an orthonormal basis for the null space of A obtained
%   from the singular value decomposition.  That is,  A*Z has negligible
%   elements, size(Z,2) is the nullity of A, and Z'*Z = I.
%
%   Uses a modified singular value tolerance formula and method.
%
%
%   See also SVD, ORTH, RANK, RREF.

 
%Paul Godfrey
%pgodfrey@intersil.com
%Feb. 7, 2001

% small is ~ 2*sqrt(eps);
% used to change zero values into small*nonzero_min()
% and used to detect all zero singular values <small


norma=sqrt(mean(mean(abs(A))));
tiny=norma*eps;

if nargin<2
   small=(1+norma)*sqrt(eps);
end

[m,n] = size(A);
if m~=n
   error('Matrix must be square!')
end

p=find(abs(A)<tiny);
if ~isempty(p)
   A(p)=0;
end

[U,S,V] = svd(A,0);
S=diag(S);
s=S;
norma=max(s);

% normalize s
smax=max(s);
if smax==0;smax=1;end
s=s/smax;
snorm=s;

% if s is all zeros
t=find(s>0);
if isempty(t);Z=V;return;end

p=find(s<tiny);
if ~isempty(p)
   s(p)=tiny;
end

%duplicate smallest non-zero value*small
%so log(s) won't goof up
%not needed because of above code
p=find(s==0);
if ~isempty(p)
   s(p)=small*min(s(t));
end

% find the largest change between consecutive s values
% this should be the spot where the zero singular values begin.
logs=log10(s);
sdifflog=-diff(logs);
%good place to break for debug
smax=max(sdifflog);
r=find(sdifflog==smax);

% use the last index
% or 0 for Inf or NaN
if min(size(r))>0
   r=r(end);
else
   r=0;
end

% the above method can still fail for matrices with
% uniformly decreasing s values but it's more consistent
% than the Matlab null command.

% the null space vectors are found
% in the last few cols of V
Z = V(:,r+1:n);

% Fix degenerate cases to match Matlab's Null output

% Scaled Identity matrix has no null space at all
if snorm==ones(n,1);  Z=zeros(n,0); end

% All zero matrix has V as it's null space
%Is this really an all zero vector? 
if max(S)<=small; Z=V; end

% save data for debugging in the nulld.mat file
%save nulld

return


%JordLook  JORD.MAT data analyses program
%
%see also: Jordan, Jord

%Paul Godfrey
%pgodfrey@intersil.com
%Feb. 7, 2001

clc
clear all
load jord
 
aa=I;
for k=1:n
    aa=aa*a1;
    nsv=nulld(aa);
    sk(1,k)=size(nsv,2);
    sk(2:(n+1),k)=svd(aa);
end

disp('Null Space size and singular values') 
Data=sk
 
return