As mentioned in [another answer](https://mathoverflow.net/a/421259/124862) 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](https://en.wikipedia.org/wiki/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$.