What method to solve a sparse complex symmetric (non-Hermitian) system?

I have a sparse system (about 78% of zero entries) that is complex and symmetric (but not Hermitian). The following figure shows the structure of the problem. The off-diagonal blocks are incidence matrices, the diagonal blocks are dense complex matrices with elements whose absolute value ranges from $$1 \times 10^{-14}$$ to $$1 \times 10^{3}$$. I've built a code using PETSc (made sure to configure --with-scalar-type=complex), it compiles and runs, but the solver either diverges or, if converged, gives nonsensical results.

The version of the code that uses dense routines (LAPACK, no PETSc) gives the expected results. The structure of the matrix generated by PETSc (see figure above) is as expected. The right hand side is zero except for a single entry (the first). The worst condition number I've seen for that system is $$1\times 10^{7}$$; but its, usually, $$1\times 10^{4}$$.

To solve that particular type of problem (sparse complex symmetric, not Hermitian), what is the recommended method (solver) to use? Can I use it with PETSc?

Such problems can be solved using an $$LDL^T$$ factorization (similar in memory and time cost to Cholesky). Your matrix is not very sparse so treating it as such may have limited benefit. I would recommend comparing sparse direct solvers (such as MUMPS and Umfpack) to the dense factorizations you've used. Complex-symmetric matrices are often difficult to precondition; I suggest checking the literature for your class of problems. Nonsymmetric Krylov methods (like GMRES) can be used, as can some symmetric methods (e.g., KSPCGSetType and MINRES variants) with appropriate preconditioners.