Here is an algorithm, it may not be a good one. I will only explain how to find a basis $e_i$, $f_i$ such that $\langle e_i, e_j \rangle = \langle f_i, f_j \rangle =0$ and $\langle e_i, f_j \rangle = c_i \delta_{ij}$ for some constants $c_i$. I will punt on explaining how to make sure that $c_1$ divides $c_2$ which divides $c_3$ and so forth. This presentation is closely based on the algorithm in Wikipedia for computing Smith normal form.

Find $e$ and $f$ so that $\langle e,f \rangle \neq 0$ **(EDIT)** and such that the lattice spanned by $e$ and $f$ is the intersection of the whole lattice with the vector space spanned by $e$ and $f$.

Set $d := \langle e,f \rangle$. Complete $(e,f)$ to a basis $(e,f,g_1, g_2, \ldots, g_{2n})$ of our lattice.

**Case 1:** We have $\langle e,g_i \rangle = \langle f, g_i \rangle =0$ for all $i$. Take $e_1=e$, $f_1=f$ and apply our algorithm recursively to the sublattice spanned by the $g_i$.

**Case 2:** For all $i$, the integer $d$ divides $\langle e,g_i \rangle$ and $\langle f, g_i \rangle =0$. Replace $g_i$ by
$$g_i - \frac{1}{d}\langle g_i, f \rangle e- \frac{1}{d}\langle e, g_i \rangle f.$$
We have now reduced to Case 1.

**Case 3:** There is some $g_i$ so that $k:=\langle e, g_i \rangle$ is not divisible by $d$. Then, for some $q$, we have $0 < k-qd < d$. Set $f'=g_i-kf$ and $g'_i=f$. Then $(e,f',g_1,g_2,\ldots, g'_i, \ldots, g_n)$ is a basis, and $\langle e, f' \rangle$ is less than $d$. Return to the beginning with this new basis.

**Case 4:** There is some $g_i$ so that $k:=\langle f, g_i \rangle$ is not divisible by $d$. Just like Case 3, with the roles of $e$ and $f$ switched.

Since $d$ decreases at every step, eventually we will hit Case 1 and be able to reduce the dimension.