Jordan matrices

mt_assign1.m

%% Matrix Theory - Assignment 1
%  Johan Astborg, HT14
%  joastbg@gmail.com

%% Illustrates the usage
%  Feel free to change the test matrix
function draft()
   
    % Test matrices
    
    %{
    M =[-9   1 -21    63 -252;
        70 -69 141  -421 1684;
        -575 575 -1149 3451 -13801;
        3891 -3891 7782 -23345 93365;
        1024 -1024 2048 -6144 24572];
    %}
    M = compan(poly(1:10));
    %M = eye(4);
    %M = zeros(3);
    
    J = jordanmatris(M);
    disp(J);
end

%% Computes the Jordan normal form of an integer matrix
%  Builds the jordan matrix using direct sum in iterative steps, where
%  the jordan block matrices are constructed using the eigenvalues and
%  their multiplicity.
%
%  Usage:
%  [ev, mult] = jordanmatris(A)
%  [ev, mult] = jordanmatris(A, tol)
%
%  See also: heltalsev, diag, jordan
function J = jordanmatris(A, tol)

    % Sets the tolerance to some proper default
	DEFAULT_TOL = 1e-6;

    % Uses the default when no tolerance value is provided as argument
    if nargin < 2
        tol = DEFAULT_TOL;
    end

    % Call the function defined below to calculate eigenvalues and 
    % the corresponding multiplicity. 
	[ev, mult] = heltalsev(A, tol);
    
    % Build the jordan matrix using direct sum in iterative steps, where
    % the jordan block matrices are constructed using the eigenvalues and
    % their multiplicity. A multiplicity > 1, means a block matrix of dim
    % greater than 1. Here we use a shifted identity matrix and the
    % particular eigenvalue on the diagonal (jordan matrix).
    J = [];
    for j=1:length(ev)
        % Creates the block jordan matrix
        b = [[zeros(mult(j)-1,1) eye(mult(j)-1)]' zeros(mult(j),1)]' + diag(ones(1,mult(j)) * ev(j));
        % Builds the jordan matrix by adding the block matrices
        J = blkdiag(J, b);
    end   
end

%% Computes the integer eigenvalues of an square integer matrix
%  heltalsev - Determines if a given integer matrix A has
%  integer eigenvalues, and then returns the eigenvalues (ev)
%  and their algebraic multiplicity (mult)%   Usage:
%  [ev, mult] = heltalsev(A)
%  [ev, mult] = heltalsev(A, tol)
%
%  See also: poly, roots, jordanmatris, eig
%
function [ev, mult] = heltalsev(A, tol)

    % Sets the tolerance to some proper default
	DEFAULT_TOL = 1e-6;

    % Uses the default when no tolerance value is provided as argument
    if nargin < 2
        tol = DEFAULT_TOL;
    end
    
    % Verify A is square matrix
    if diff(size(A)) ~= 0 || ~ismatrix(A)
        throw(MException('heltalsev:badarg', ' A must be square matrix.'));
    end
    
    % No need to recompute this all the time
    detA = det(A);
    
    % Verify, using the determinant of A, that A is an integer matrix.
    % This works because the determinant is integer for integer matrices.
    if abs(round(detA) - detA) > tol
        throw(MException('heltalsev:badarg', 'Only integer matrices supported.'));
    end
    
    % Adjust for the case when all eigenvalues are the same or when the 
    % determinant and trace are both zero. The determinant is the product 
    % of all eigenvalues and the trace is the sum. If either equals zero, 
    % the eigenvalues are all zero (with multiplicity according to dim(A)).
    len = length(A);
    if (detA^(1/len) == trace(A)/len) || (detA == 0 || trace(A) == 0)  
        ev = trace(A)/len;
        mult = len;
        return;
    end
    
    % Check the case when A is diagonal. Then there is no need to 
    % calculate the eigenvalues (below), just return the diagonal elements.
    if isdiag(A)
        sumA = sum(A);
        ev = sumA;
        mult = arrayfun(@(x)sum(sumA==x), unique(sumA));
        return;
    end
    
    % Characteristic polynomial of A
    p = poly(A);
    
    % Finds the roots of p (i.e. the eigenvalues)
    r = roots(p);

    % Check so eigenvalues meet the tolerance critiria (are integers)
    if sum(abs(round(r) - r)) > tol * length(r)
        throw(MException('heltalsev:error', 'Eigenvalues are not integers (you may try to change tol).'));
    end
    
    % The eigenvalues are rounded (this is fine, hence passed test above)
    ev = sort(round(r));
    
    % Calculates the multiplicities
    mult = arrayfun(@(x)sum(round(r) == x), unique(round(r)))';
end