MATL, 29 bytes, non-competing

1$Yn1$ZQYotdUZS~0hXdF1hYSwXd+

Try it online!

This is my first MATL submission so there are bound to be improvements. I spent a while learning it and only at the end did I remember that this might not have worked using the MATL from May 7, 2016. Sure enough, I checked out the penultimate commit to that day and it did not run.

I would have liked to use diag but it seems MATL only supports the single argument version. The second argument would be needed to place values along the superdiagonal (or any other diagonal for different problems).

Octave, 60 bytes

@(a)diag(1-sign(diff(v=round(roots(poly(a)))).^2),1)+diag(v)

Try it online!

A port of my Mathematica solution.

J, 78 71 bytes

1(#\.|."_1#{."1],.2=/\,&_)@>@{[:p.@>[:-&.>/ .(+//.@(*/)&.>)],&.>[:-@=#\

Try it online!

The two challenging portions of this task, getting the eigenvalues and performing the diagonalization, ending up both taking about the same number of bytes. These were disallowed by the rules but in case any are curious, J has builtins for QR decomposition (128!:0) as well as LAPACK addons which could be used to find eigenvalues.

Explanation (Outdated)

There are two major portions to this verb: finding the eigenvalues and performing the diagonalization. First, in order to find the eigenvalues, the roots of the characteristic polynomial for the input matrix will have to be found. Using the same input matrix from the example,

   ] m =: _4 ]\ 5 4 2 1 0 1 _1 _1 _1 _1 3 0 1  1 _1 2
 5  4  2  1
 0  1 _1 _1
_1 _1  3  0
 1  1 _1  2

The characteristic polynomial for a matrix M can be found using |M - λI| = 0 where I is the identity matrix with the same dimensions as M. The expression M - λI can be modeled in J by boxing each element in M with a -1 if it is on the diagonal else otherwise a 0. Each box represents a polynomial in coefficient form.

   (],&.>[:-@=#\) m
┌────┬────┬────┬────┐
│5 _1│4 0 │2 0 │1 0 │
├────┼────┼────┼────┤
│0 0 │1 _1│_1 0│_1 0│
├────┼────┼────┼────┤
│_1 0│_1 0│3 _1│0 0 │
├────┼────┼────┼────┤
│1 0 │1 0 │_1 0│2 _1│
└────┴────┴────┴────┘

The determinant in J is -/ .* however, that operates on numbers, not boxed polynomials. Instead of multiplication, the polynomial product is needed which can be found using convolution ([:+//.*/). Folded subtraction is still used, and both these verbs need to operate within boxes so under (&.) unbox (>) is used.

   ([:-&.>/ .(+//.@(*/)&.>)],&.>[:-@=#\) m0
┌───────────────┐
│32 _64 42 _11 1│
└───────────────┘

These are the coefficients of the characteristic polynomial. The roots can be found using p. which converts the representation of a polynomial between coefficients and roots form.

   ([:p.@>[:-&.>/ .(+//.@(*/)&.>)],&.>[:-@=#\) m0
┌─┬───────┐
│1│4 4 2 1│
└─┴───────┘

The roots are [4, 4, 2, 1] and those are the eigenvalues of M.

Second, the diagonalization must be performed. Each continuous pair of values is tested for equality.

   (2=/\]) 4 4 2 1
1 0 0

A zero is appended and those values are columinized together with the eigenvalues.

   (],.0,~2=/\]) 4 4 2 1
4 1
4 0
2 0
1 0

Then each row is padded to the same length as the number of eigenvalues to form a square matrix.

   (#{."1],.0,~2=/\]) 4 4 2 1
4 1 0 0
4 0 0 0
2 0 0 0
1 0 0 0

Finally each row is shifted right with values falling off at the right and zeros being pushed in on the left. The first row is shifted zero times, the second once, the third twice, and so on.

   (-@i.@#|.!.0"_1#{."1],.0,~2=/\]) 4 4 2 1
4 1 0 0
0 4 0 0
0 0 2 0
0 0 0 1

The output is the Jordan decomposition of M.

Mathematica, 140 139 105 bytes

Total[DiagonalMatrix@@@{{#},{1-Sign[Differences@#^2],1}}]&@(x/.Solve[#~CharacteristicPolynomial~x==0,x])&

I just found the builtin DiagonalMatrix which allows for an easy way to place the 0s and 1s along the superdiagonal.

Usage

Sage, 79 bytes

lambda A:block_diagonal_matrix([jordan_block(*r)for r in A.charpoly().roots()])

Try it online

Since nobody else is posting solutions, I might as well go ahead and post one.

A.charpoly.roots() computes the roots (and algebraic multiplicities) of the characteristic polynomial of A (aka the eigenvalues and multiplicities). jordan_block constructs a Jordan block from the given root and multiplicity. block_diagonal_matrix forms a matrix with the Jordan blocks on the diagonal, which is exactly the definition of a Jordan normal form.