The set of eigenvalues computed with `eigvals(A,B)` must be self-conjugated for real `A` and `B`

I came accross this rather minor issue when implementing software for spectrum assignment. Sorry if it has been already addressed in this forum.

It is known that for two real square matrices `A` and `B` of the same dimension, the set of generalized eigenvalues of the pair `(A,B)` is self conjugated (i.e., symmetric with respect to the real axis: thus, if `λ` is an eigenvalue then `conj(λ)` is also an eigenvalue). If the vector of eigenvalues represents an input parameter to some function (e.g., in my software for spectrum assignment), then the self-conjugation condition is usually tested to ensure real results for real input data. Unfortunately, the computed eigenvalues using the function `eigvals` do not usually satisfy the self-conjugation condition (due to roundoff errors) as can be simply seen from the following example, where the vector of eigenvalues `ev` has been generated with `ev = eigvals(rand(10,10),rand(10,10)) `:
```

julia> ev =
[-1.358678198251497 - 0.7515055758250914im
-1.3586781982514966 + 0.7515055758250913im
-0.5706054730423562 - 0.3751717714923926im
-0.5706054730423562 + 0.37517177149239256im
-0.1846452066231087 + 0.0im
0.4932925288129128 - 0.2915837932828266im
0.4932925288129128 + 0.2915837932828265im
0.9456289869357847 + 0.0888486972750476im
0.9456289869357848 - 0.08884869727504759im
5.861894498120641 + 0.0im]

julia> ev[imag.(ev) .> 0] == conj(ev[imag.(ev) .< 0]) # no exact fit
false

julia> ev[imag.(ev) .> 0] ≈ conj(ev[imag.(ev) .< 0]) # approximate fit accounts for the effect of roundoff
true
```

The symmetrization of `ev` is straightforward:
```
ev[imag.(ev) .> 0] = conj(ev[imag.(ev) .< 0]); ev
10-element Array{Complex{Float64},1}:
-1.358678198251497 - 0.7515055758250914im
-1.358678198251497 + 0.7515055758250914im
-0.5706054730423562 - 0.3751717714923926im
-0.5706054730423562 + 0.3751717714923926im
-0.1846452066231087 + 0.0im
0.4932925288129128 - 0.2915837932828266im
0.4932925288129128 + 0.2915837932828266im
0.9456289869357848 + 0.08884869727504759im
0.9456289869357848 - 0.08884869727504759im
5.861894498120641 + 0.0im

julia> ev[imag.(ev) .> 0] == conj(ev[imag.(ev) .< 0]) # exact fit achieved
true
```

Note that MATLAB does always ensure the self-conjugation condition in the case of real data. I think it would be nice to enforce this feature also in Julia.

The issue centers on the eigvals(A,B) entry point; begin by tracing how it handles real square matrix pairs and how eigenvalues are returned. Done means results for real A and B satisfy exact self-conjugation, with coverage for the reported generalized-eigenvalue example.

Written by the indexing model from the issue text.