Basis for unstable manifold of a matrix

Question

Given a square matrix A, how can I generate a basis for the generalized eigenspace corresponding to all eigenvectors $\lambda_i$ such that $\vert \lambda_i \vert > 1$? I.e., if $A$ is $n \times n$ and acts on $\mathbb{R}^n$, then $\mathbb{R}^n = E^u \oplus E^1 \oplus E^s$, corresponding to the action of eigenvalues of modulus greater than / less than / equal to 1.

If A is diagonalizable, this can easily be accomplished with the Eigensystem command:

{evals, evecs}=Eigensystem[A];
Take[evecs, Length@Select[evals, Abs[#] > 1 &]]

This approach relies on the Eigenvalue command returning the eigenvalues (and hence eigenvectors) in decreasing order of their moduli. Also note that MMA returns the eigenvalues with repetition according to their algebraic multiplicity, and that the eigenvectors are returned as "row vectors" (technically, a list of vectors).

My problem arises when A is not diagonalizable. The natural function to look at would be JordanDecomposition, which retuns the generalized eigenvectors (N.B. as columns in a matrix!) and the Jordan Canonical Form. However, the ordering is a bit...abstract. As best as I've been able to tell from taking the JordanDecomposition of random matrices, the ordering of the Jordan blocks seems to be some combination of the order in which roots are returned by the Root command for irreducible factors of the characteristic polynomial, the ordinary Greater sorting, the size of the Jordan block, and whether it corresponds to a complex eigenvalue.

The kludge that I've been able to come up with is the following:

{s, j} = JordanDecomposition[A];
bigevals=Union@Select[Eigenvalues[A], Abs[#] > 1 &];
cols= Flatten[ Position[Diagonal[j], #]& /@ bigevals];
s[[All, cols]]

This seems to return the appropriate generalized eigenvectors.

The Questions

1. Surely there's a better way? The above kludge seems...kludgy. Perhaps something like NullSpace[MatrixPower[A-lambda IdentityMatrix[First@Dimensions@A], n]] upon detecting that an eigenvalue lambda has algebraic multiplicity n, i.e., going directly to the definition?
2. For my application, the matrix A will only ever be real valued. Thus all complex eigenvalues occur as complex conjugates (perhaps with multiplicity!). I would prefer to have a real basis. For example, if A={{0,2},{-2,0} then the eigenvalues are $\lambda_{\pm} = \pm 2 i$. Then $\left\{ \begin{bmatrix} 1 \\ 0 \end{bmatrix} , \begin{bmatrix} 0 \\ 1 \end{bmatrix} \right\}$ is a basis for the unstable manifold of $A$. Nullspace[MatrixPower[A,2]+4 IdentityMatrix[2]] gives the correct result here. Note that this question is a tertiary concern at best; there are other ways of dealing with this in "post-production."

Answer 1

Remember that the diagonal elements of the Jordan normal form are the eigenvalues. There is no need to compute the eigenvalues and the decomposition. Instead, the following code creates a picker list from the diagonal elements.

A = 1/8 {{27, 48, 81}, {-6, 0, 0}, {1, 0, 3}};
{s, j} = JordanDecomposition[A];
Pick[Transpose[s], Map[Abs[N[#]] > 1 &, Transpose[j, {1, 1}]]]

From there, you might want to use Orthogonalize …