Vectorized complex division in Eigen uses the naive algorithm, which is susceptible to over- and underflow.

The complex division packet op pdiv uses the naive algorithm to compute x + iy = (a + ib) / (c + id) in the following manner:

denom = c**2 + d**2
x = (a * c + b * d )/ denom
y = (b * c - a * d) / denom

for each complex number in the output packet. This algorithm is prone to over- and underflow when c and d are very large or very small, respectively.

This causes certain tests to fail in the LAPACK linear solver test suite when run through the wrappers in the eigen/lapack directory to test the solvers implemented in Eigen.

It is worth noting that the failure disappears if vectorization is disabled, because the C standard implements a more robust algorithm for scalar complex division, as described on page 21 in: http://www.open-std.org/jtc1/sc22/wg11/docs/n435.pdf

I would like to propose that we adopt an algorithm at least as robust as the one in the C standard. For a survey of a variety of algorithms for complex division, see: https://arxiv.org/pdf/1210.4539.pdf

Update: Apparently GCC does not implement the more robust algorithm for scalar division either, so we should probably provide our own implementation in this case.