Fix of bug in gelsd! for complex arrays and some more.

[andreasnoack]

Added some eps methods for complex floats. Fixed a bug in dgesld for complex types and added input argument for singularity value. Updated \ accordingly. Added some methods to handle right and left hand side of different type. Restricted methods in operators apply for Numbers only.

[ViralBShah]

The xgelsd stuff looks fine. What is the rationale for changing the behaviour of eps() for Complex? About the changes in operators.jl, I think it is actually a good idea to restrict some of those to Number, but it would also be annoying to implement each operator for every new type. Also, I wonder if there might be unintended side-effects.

@JeffBezanson @StefanKarpinski What are your thoughts on the changes to operators.jl?

It would have been easier to merge this without the changes to operators.jl. Let's see what others have to say at this point, before reverting or creating a new pull request.

[nolta]

The operator stuff breaks the test suite:

     * linalg
no method .>(Array{Float64,1},Float64)

[andreasnoack]

Sorry for not testing before sending the request. It passes all tests now.

The error was caused by .> not being defined for arrays after my changes. I have moved the definition .> (x,y) = y.<x next to > (x,y) = y < x which I think makes sense. If the restrictions in operators are problematic I can remove them, but I extrapolated from the discussion about \ and restricted the operators that did not make sense for matrices.

The changes to eps should not break anything. I changed abs to real for efficiency and added methods similar to those of eps(Float64) and eps(Float32).

[ViralBShah]

eps(z::Complex) should be reverted to the original definition using abs. It is not clear what eps should do for complex numbers, and how meaningful that is. While we do have it, I think it should be matlab compatible, which is what the original behaviour was.

[andreasnoack]

But I do not think the behaviour has changed. The output is the same for eps(real(x::Complex)) and eps(abs(x::Complex)). The only difference is that the former avoids two squarings and a square root.

[StefanKarpinski]

They're not the same:

julia> eps(abs(0.+1im))
2.220446049250313e-16

julia> eps(real(0.+1im))
5.0e-324

Matlab gives the former answer, not the latter.

[andreasnoack]

I buy that argument! I'll revert it.

[nolta]

Defining eps(::Complex) seems like a matlab feature we can do without.

@andreasnoackjensen, would you mind splitting your patch into 3 parts?:

1. the xgelsd stuff,
2. a narrowly tailored fix for Error: stack overflow when solving linear system for different types #1418, which also adds randn(3,3)\complex(randn(3)) to the test suite,
3. the rest of the operator changes.

[andreasnoack]

Fair enough. I'll close this one and make some new ones.

[StefanKarpinski]

Yes, eps(::Complex) is a bit iffy since it's not clear what it means, whereas eps for floating-point is very well-defined except for the weird edge cases. We disagree with Matlab at realmax() – Matlab has eps(2.2250738585072014e-308) == 4.9407e-324, whereas we give Inf. I think our answer is technically more correct and makes the eps() function much, much faster to compute for all values since all you do is increment a register as an integer and then subtract as floats from the input value. The Matlab and C versions have crazy branches and other nonsense. Doesn't seem worth it.

[ViralBShah]

I actually think we should not have eps() for Complex, or make it throw an DoesNotMakeSenseException. It was probably implemented in Matlab for some kind of customer convenience.

One could define eps(::Complex) to be the smallest complex number to add to the given complex number, so that the result is next closest to the input. However, you can have multiple answers here.

[StefanKarpinski]

LOL @ DoesNotMakeSenseException although it's also not a terrible idea.

[StefanKarpinski]

We could also add WTFDudeException and WoahThereNellyException.

[toivoh]

I guess DomainError might be appropriate?

[ViralBShah]

DomainError does seem more appropriate.