% Computation of the Jordan:  C^{-1} A C = J decomposition

 n = 4
 A = diag(1:n); o = ones(n,1);
 W = full( spdiags([-o 2*o -o],-1:1, n,n) ); W(n,n)=1;  % Similarity

 A(r,r) = inv(W)*gallery('jordbloc',length(r),4)*W; A=round(A)
 [X0 D]=eig(A);
 EIgenvectors = X0, Single_eigenvalue = D(1,1)

disp('============================= eigs solve for Jordan ==============')
%%%%%%%%%%%%%%%%%% NOT accurate enough %%%%%%%%%%%%
opt.disp=0;  opt.tol=1.0e-12; opt.maxit=500;
[x1 d]=eigs( A,1,'sm',opt)
[x2 d]=eigs( (A - 4*eye(4))^2,1,'sm',opt)
  y = (A - 4*eye(4))*x2 ;
  [wp jp]=max(abs(x1)); [vn jn]=max(abs(y));
  x2 = x1(jp)/y(jn)*x2; 
[x3 d]=eigs( (A - 4*eye(4))^3,1,'sm',opt)
  y = (A - 4*eye(4))*x3 ;
  [wp jp]=max(abs(x2)); [vn jn]=max(abs(y));
  x3 = x2(jp)/y(jn)*x3; 
[x4 d]=eigs( (A - 4*eye(4))^4,1,'sm',opt)
  y = (A - 4*eye(4))*x4 ;
  [wp jp]=max(abs(x3)); [vn jn]=max(abs(y));
  x4 = x3(jp)/y(jn)*x4; 
A
X=[x1 x2 x3 x4]
Jordan_J = (X\A)*X

disp('============================= QR solve for Jordan ==============')

%_____________________ QR solve
x1=X0(:,3);
 [Q R]=qr( A-4*eye(n) );
 rhs = Q'*x1;  x2=R(1:n-1,:) \ rhs(1:n-1)
 rhs = Q'*x2;  x3=R(1:n-1,:) \ rhs(1:n-1)
 rhs = Q'*x3;  x4=R(1:n-1,:) \ rhs(1:n-1)
A
X=[x1 x2 x3 x4]
Jordan_J = (X\A)*X

