We can set up the problem as follows:
n = 3
jordan_type = (2,1)
assert sum(jordan_type) == n
import itertools
A_coeff_names = {(i,j) : 'a_{}_{}'.format(i,j) for (i,j) in itertools.combinations(range(n),2)}
R = PolynomialRing(QQ, names=A_coeff_names.values())
A_coeff = {idx : R(A_coeff_names[idx]) for idx in A_coeff_names}
A = Matrix(R, n, A_coeff)
J = block_diagonal_matrix([jordan_block(0, s) for s in jordan_type], subdivide=False)
A necessary condition for $gAg^{-1} = J$ is that the rank of $A^k$ equals the rank of $J^k$ for $1 \leq k < n$.
(And since $A$ has all zeros as eigenvalues, this will give the right Jordan type.)
In particular it is necessary that the $(\operatorname{rank}(J^k)+1)$-minors of $A^k$ vanish; these are polynomial equations upon the coefficients of $A$:
A_coeff_eqns = []
J_power = identity_matrix(QQ, n)
A_power = identity_matrix(R, n)
for k in range(n):
r = J_power.rank()
# need (r+1) x (r+1) minors to vanish
A_coeff_eqns.extend([eqn for eqn in A_power.minors(r+1) if eqn != 0])
J_power *= J
A_power *= A
We can solve these equations symbolically, substitute a solution into $A$, change to an exact ring (the fraction field of the polynomial ring in the new variables), calculate the Jordan form (to be sure), and change the names of the variables back to the originals where possible:
A_coeff_symb = map(SR, A_coeff_names.values())
A_coeff_eqns_symb = map(SR, A_coeff_eqns)
for sol in solve(A_coeff_eqns_symb, A_coeff_symb):
# substitute symbolic solution
A_new = A.change_ring(SR).apply_map(lambda x: x.subs(sol))
A_new_vars = list(set(sum([list(x.variables()) for x in A_new.list()], [])))
# change ring to fraction field of polynomial ring in the new variables
A_new_ring = PolynomialRing(QQ, names=A_new_vars).fraction_field()
A_new_rat = A_new.change_ring(A_new_ring)
if A_new_rat.jordan_form(subdivide=False) != J:
continue
# change names of variables back to the originals where possible
A_final_subs = {A_new_rat[i,j] : A_coeff[(i,j)] for (i,j) in A_coeff if A_new_rat[i,j] in A_new_vars}
A_final = A_new_rat.apply_map(lambda x: x.subs(A_final_subs))
print A_final
print
print A_final.jordan_form(subdivide=False, transformation=True)[1]
print
print '---'
print
Output:
[ 0 a_0_1 a_0_2]
[ 0 0 0]
[ 0 0 0]
[ a_0_1 0 0]
[ 0 1 1]
[ 0 0 a_0_1/(-a_0_2)]
---
[ 0 0 a_0_2]
[ 0 0 a_1_2]
[ 0 0 0]
[a_0_2 0 1]
[a_1_2 0 0]
[ 0 1 0]
---
Beware that this assumes Sage can solve the polynomial equations symbolically, which may not be true for large n. To give Sage an easier time (though the solutions may become less pretty), you can replace the list of equations by a Groebner basis for the same ideal:
A_coeff_eqns = R.ideal(A_coeff_eqns).groebner_basis()
Can you give a more explicit example?
Sure. Say I have an arbitrary uppertriangular matrix of Jordan type $(2,1)$.
$$ A = \begin{pmatrix} 0 & A_{0} & A_{1} \\ 0 & 0 & A_{2} \\ 0 & 0 & 0 \end{pmatrix} $$
By default
J_d, T_d = A.jordan_form(transformation=True)gives$$ J_d = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} \qquad T_d = \begin{pmatrix}A_{0} A_{2} & A_{1} & 0 \\ 0 & A_{2} & 0 \\ 0 & 0 & 1\end{pmatrix} $$
However, I expect $$J = \begin{pmatrix}0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{pmatrix}$$
In this case, I know that $A_2$ must be zero for $A$ to have Jordan type $(2,1)$ so I let $A' = A\big|_{A_2 = 0}$ and do
J, T = A'.jordan_form(transformation=True)getting the desired$$T = \begin{pmatrix}A_{0} & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & -\frac{A_{0}}{A_{1}}\end{pmatrix}$$
In general, the constraints may be more complicated than setting some entries to zero. Either way, I want Sage to figure out those constraints for me, and solve something like $g A g^{-1} = J$ for a $g$, given variable $A$ and fixed $J$.
If I try
A.jordan_form()in a quotient ring, I get aNotImplementedError:(An alternative approach is to write $A = g^{-1}Jg$ with $g$ symbolic, and solve for $g$ such that $A$ is strictly upper triangular ($n(n+1)/2$ polynomial equations for the $n^2$ entries of $g$). Would that suffice, or not? What quotient ring did you want to work in?
Yes, it might. How do I do it, exactly?
I would work in the ring $R/I$, where $R$ is the polynomial ring generated by matrix coefficients, and $I$ is the ideal of the Jordan form. For instance, if I have a 5 by 5 matrix $A$ of type $(3,2)$ then