RFC: Some linalg changes.

A new general matrix factorization method. Some changes to the treatment of symmetric and hermitian matrices and some bug fixes. And some tests.

This is a follow up on the discussion in #3295. The main new thing is `factorize` which tries to guess a good factorization depending on the structure of matrix. For example

```julia
julia> A=randn(5,5);

julia> typeof(factorize(A))
LU{Float64}

julia> typeof(factorize(A + A'))
BunchKaufman{Float64}

julia> typeof(factorize(A'A))
Cholesky{Float64}

julia> typeof(factorize(A[:,1:2]))
QRPivoted{Float64}
```
For complex indifinite matrices it is necessary to distinguish between symmetric and hermitian matrices (the Bunch-Kaufman case) and therefore I have added a `Symmetric` type which also works for real matrices. Hence `Hermitian` is now a complex matrix only type. I have also added the `bkfact` method for constructing the Bunch-Kaufman decoposition directly.

To handle least squares, I had to reimplement the functionality of `xgelsy` in Julia because I wanted to be able to solve rank deficient problems. You can now reuse a QR based "predicter":
```julia
julia> mX = randn(10,2)*randn(2,3);y1 = randn(10);y2 = randn(10);

julia> xf = factorize(mX);

julia> xf\y1
3-element Float64 Array:
 -0.18619   
  0.811663  
  0.00219657

julia> xf\y2
3-element Float64 Array:
  0.0448524
 -0.197673 
  0.0139636
```
I got some random failures in the tests but I couldn't get them if I ran the tests with `include` from a Julia session.