As mentioned in another answer the fact that $A$ and $A^T$ are similar corresponds to a commutative diagram of modules over the polynomial ring $F[\lambda]$: $$ \begin{array}{c} 0 & \to & F[\lambda]\otimes V & & \stackrel{\lambda I - A}{\longrightarrow} & F[\lambda]\otimes V & \to & V_{A} & \to & 0\\ && \downarrow~R &&& Q~\downarrow & & \downarrow P \\ 0 & \to & F[\lambda]\otimes V & & \stackrel{\lambda I - A^T}{\longrightarrow} & F[\lambda]\otimes V & \to & V_{A^T} & \to & 0 \end{array} $$ where:

• $V$ denotes that $n$-dimensional column vector space over $F$
• For a square matrix $B$ we use $V_B$ to denote $V$ as a module over $F[\lambda]$ on which $\lambda$ acts by $B$.
• The maps $P$, $Q$, $R$ are isomorphisms of $F[\lambda]$ modules.

The matrix $P$ is what we need to determine.

The idea is that once we determine $Q$ and $R$, then $P$ can be easily calculated.

The Wikipedia page for the Smith Normal Form sketches an algorithm from which one can compute invertible matrices $F$ and $G$ over $F[\lambda]$ such that $F(\lambda u - A)G$ is a diagonal matrix over $F[\lambda]$. It follows that $G^T(\lambda u -A^T)F^T$ is the same diagonal matrix.

We can now take $R^{-1}=G(F^T)^{-1}$ and $Q=(G^{T})^{-1}F$ to get the identity $Q(\lambda I - A)=(\lambda I - A^{T})R$ as required. As mentioned above, now that we have calculated $Q$ and $R$, we can calculate $P$.

As mentioned in another answer the fact that $A$ and $A^T$ are similar corresponds to a commutative diagram of modules over the polynomial ring $F[\lambda]$: $$ \begin{array}{c} 0 & \to & F[\lambda]\otimes V & & \stackrel{\lambda I - A}{\longrightarrow} & F[\lambda]\otimes V & \to & V_{A} & \to & 0\\ && \downarrow~R &&& Q~\downarrow & & \downarrow P \\ 0 & \to & F[\lambda]\otimes V & & \stackrel{\lambda I - A^T}{\longrightarrow} & F[\lambda]\otimes V & \to & V_{A^T} & \to & 0 \end{array} $$ where:

• $V$ denotes that $n$-dimensional column vector space over $F$
• For a square matrix $B$ we use $V_B$ to denote $V$ as a module over $F[\lambda]$ on which $\lambda$ acts by $B$.
• The maps $P$, $Q$, $R$ are isomorphisms of $F[\lambda]$ modules.

The matrix $P$ is what we need to determine.

The idea is that once we determine $Q$ and $R$, then $P$ can be easily calculated.

The Wikipedia page for the Smith Normal Form sketches an algorithm from which one can compute invertible matrices $F$ and $G$ over $F[\lambda]$ such that $F(\lambda u - A)G$ is a diagonal matrix over $F[\lambda]$. It follows that $G^T(\lambda u -A^T)F^T$ is the same diagonal matrix.

We can now take $R^{-1}=G(F^T)^{-1}$ and $Q=(G^{T})^{-1}F$ to get the identity $Q(\lambda I - A)=(\lambda I - A^{T})R$ as required. As mentioned above, now that we have calculated $Q$ and $R$, we can calculate $P$.

Update: Since it may not be immediately obvious how to calculate the matrix of $P$ with entries in $F$, here is an explicit calculation which uses the matrix $Q$ over $F[\lambda]$.

Note that if $e_i$ denotes the $i$-th standard basis vector of $V$, then under the horizontal map from $F[\lambda]\otimes V$ to $V_A$ or $V_B$, $1\otimes e_i$ goes to $e_i$.

Under $Q$, the image of $1\otimes e_i$ in $F[\lambda]\otimes V$ is $E_i=\sum_{k=1}^n Q_{k,i} \otimes e_k$, where $Q_{k,i}$ is the $(k,i)$-th entry of $Q$ and is a polynomial in $\lambda$. The image of $E_i$ in $V_B$ is $f_i=\sum_{k=1}^n Q_{k,i}(A)e_k$. We then write $f_i=\sum_{l=1} P_{l,i} e_l$ to get the matrix $P$. (It is a bit curious that it only depends on $Q$!) Note that to actually calculate $P$, we need to calculate the matrix entries of $Q_{k,i}(A)$ for all $i$ and $k$.