We describe a matrix diagonalization algorithm for complex symmetric (not Hermitian) matrices, A = A^T , which is based on a two-step algorithm involving generalized Householder reflections based on the indefinite inner product _* = sum_i u_i v_i. This inner product is linear in both arguments and avoids complex conjugation. The complex symmetric input matrix is transformed to tridiagonal form using generalized Householder transformations (first step). An iterative, generalized QL decomposition of the tridiagonal matrix employing an implicit shift converges toward diagonal form (second step). The QL algorithm employs iterative deflation techniques when a machine-precision zero is encountered “prematurely” on the super-/sub-diagonal. The algorithm allows for a reliable and computationally efficient computation of resonance and antiresonance energies which emerge from complex-scaled Hamiltonians, and for the numerical determination of the real energy eigenvalues of pseudo–Hermitian and PT-symmetric Hamilton matrices. Numerical reference values are provided.