My research is in electron spin resonance, where we deal with the time-independent stochastic Liouville equation and end up with the task of finding the eigenvectors and eigenvalues of complex symmetric matrices (size ~ 10000x10000, ~30% elements non-zero).
We have been using Lanczos algorithm with conjugate gradients/QMR, but lose orthogonality of the Lanczos vectors (and hence the eigenvectors).
Note: Here orthogonality between complex vectors $a$, $b$ means $a^{T}b=0$, not $a^\dagger b=0$.
Therefore, the next step seems to be increasing the precision from double to quadruple, along with implementing Lanczos algorithm with partial reorthogonalization for sparse, complex symmetric matrices.
I want a link to a FORTRAN package that implements Lanczos algorithm with partial reorthogonalization for large, sparse, complex symmetric matrices. My Google search seemed futile.