JuliaLang/julia — pull request #2308 by andreasnoackjensen (2013-02-15 07:23:39+00:00)
Title: WIP:New design of linalg (factorizations)

Based on the discussion in #2212 I have rewritten some of the linalg functionality. Please share your thoughts. The main points of the design are

Matlab style methods are retained for now

Removal of `factors` methods and introduction of `ref` methods for factorizations object. Symbols are used for referencing. E.g.
```julia
julia> A=randn(5,3)
5x3 Float64 Array:
-1.84808 -1.56926 -0.394121
2.5735 1.11443 0.0477358
0.751048 0.920844 0.756794
-0.987947 -1.39258 -1.59187
0.684266 1.28512 -1.37579

julia> Aqr=QRDense(copy(A));

julia> Aqr[:Q]*Aqr[:R]
5x3 Float64 Array:
-1.84808 -1.56926 -0.394121
2.5735 1.11443 0.0477358
0.751048 0.920844 0.756794
-0.987947 -1.39258 -1.59187
0.684266 1.28512 -1.37579
```

Factorizations work similarly to the matrix they were constructed from. No magic.

New `QRDensQ` and `HessenbergQ` types which are subtypes of AbstractMatrix can be used for multiplication.
```julia
julia> typeof(Aqr[:Q])
QRDenseQ{Float64}
```

For types with pivoting you can get either a permutation vector of matrix
```julia
julia> A=randn(5,5)
5x5 Float64 Array:
0.466243 1.15806 1.06016 0.479681 -0.0639797
0.603621 -0.94684 -0.0267892 0.173938 -0.458464
0.87668 -2.22612 0.23515 2.37495 0.0723865
0.849353 0.180417 1.08095 -0.45903 -0.993289
1.17308 -0.644301 0.202831 1.50354 -0.065989

julia> LU=LUDense(copy(A));

julia> LU[:P]'LU[:L]*LU[:U]
5x5 Float64 Array:
0.466243 1.15806 1.06016 0.479681 -0.0639797
0.603621 -0.94684 -0.0267892 0.173938 -0.458464
0.87668 -2.22612 0.23515 2.37495 0.0723865
0.849353 0.180417 1.08095 -0.45903 -0.993289
1.17308 -0.644301 0.202831 1.50354 -0.065989

julia> LU[:L][invperm(LU[:p]),:]*LU[:U]
5x5 Float64 Array:
0.466243 1.15806 1.06016 0.479681 -0.0639797
0.603621 -0.94684 -0.0267892 0.173938 -0.458464
0.87668 -2.22612 0.23515 2.37495 0.0723865
0.849353 0.180417 1.08095 -0.45903 -0.993289
1.17308 -0.644301 0.202831 1.50354 -0.065989
```

Removal of all fact methods. Factorizations are constructed from constructors. The fact methods were kind of constructors so I think it is more clear just to use e.g. `QRDense` instead of `qrfact`.

New Triangular and Hermitian types for dispatch. This eases the choice of algorithm and simplifies some of the functions for general matrices. E.g.
```julia
julia> B=randn(5,2)
i5x2 Float64 Array:
-1.67605 0.0792652
-2.1588 0.0286816
-0.252934 -0.082285
0.777874 1.55441
-1.9374 0.949172

julia> issym(B/(B'B)*B')
false

julia> issym(Hermitian(B/(B'B)*B'))
true

julia> eigvals(B/(B'B)*B') # non symmetric solver
5-element Float64 Array:
-1.11022e-16
1.0
0.0
1.0
-1.74743e-17

julia> eigvals(Hermitian(B/(B'B)*B')) # symmetric solver
5-element Float64 Array:
-1.96396e-16
-1.24191e-16
-7.64151e-17
1.0
1.0
```

New symmetric tridiagonal eigen solvers and methods for subsets of eigenvalues
```julia
julia> eigvals(Hermitian(B/(B'B)*B'), 4, 5)
2-element Float64 Array:
1.0
1.0

julia> eigvals(Hermitian(B/(B'B)*B'), 0.9, 1.1)
2-element Float64 Array:
1.0
1.0
```

Removal Bool argument in eig. Just use eigvals.