From 186287a2a891baf0edf658220b43e3077af00904 Mon Sep 17 00:00:00 2001
From: Andreas Noack Jensen <andreasnoackjensen@gmail.com>
Date: Fri, 17 Jan 2014 12:55:28 +0100
Subject: [PATCH] Add julia implementation of the QR algorithm. This includes
 elementary reflectors. Add Float16 and BigFloat to tests and test promotion.
 Cleaup promotion in LinAlg. Avoid promotion in mutating ! functions. Make
 Symmetric, Hermitian and QRs immutable. Make thin a keyword argument in
 svdfact. Remove cond argument in sqrtm.

---
 NEWS.md                      |   7 +
 base/deprecated.jl           |   1 +
 base/exports.jl              |   1 -
 base/linalg.jl               |   1 -
 base/linalg/bunchkaufman.jl  |   8 +-
 base/linalg/dense.jl         |  81 +++---
 base/linalg/factorization.jl | 462 ++++++++++++++++++++++++-----------
 base/linalg/generic.jl       |  64 +++++
 base/linalg/lapack.jl        |  12 +-
 base/linalg/symmetric.jl     |  93 ++++---
 base/linalg/tridiag.jl       |  16 +-
 test/linalg.jl               | 338 ++++++++++++++-----------
 12 files changed, 686 insertions(+), 398 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index eed76349576d4..1a0e73dc293ea 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -110,6 +110,8 @@ Library improvements
 
       * Balancing options for eigenvector calculations for general matrices ([#5428]).
 
+      * Mutating linear algebra functions no longer promote ([#5526]).
+
     * Sparse linear algebra
 
       * Faster sparse `kron` ([#4958]).
@@ -147,6 +149,8 @@ Library improvements
 
       * LU factorization ([#5381] and [#5430])
 
+      * QR factorization ([#5526])
+
   * New function `deleteat!` deletes a specified index or indices and
     returns the updated collection
 
@@ -177,6 +181,8 @@ Deprecated or removed
 
   * `myindexes` has been renamed to `localindexes` ([#5475])
 
+  * `factorize!` is deprecated in favor of `factorize`. ([#5526])
+
 [#4042]: https://github.com/JuliaLang/julia/issues/4042
 [#5164]: https://github.com/JuliaLang/julia/issues/5164
 [#4026]: https://github.com/JuliaLang/julia/issues/4026
@@ -220,6 +226,7 @@ Deprecated or removed
 [#5025]: https://github.com/JuliaLang/julia/pull/5025
 [#4888]: https://github.com/JuliaLang/julia/pull/4888
 [#5475]: https://github.com/JuliaLang/julia/pull/5475
+[#5526]: https://github.com/JuliaLang/julia/pull/5526
 
 Julia v0.2.0 Release Notes
 ==========================
diff --git a/base/deprecated.jl b/base/deprecated.jl
index 9d962ec267483..34f03cb48e96e 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -178,6 +178,7 @@ export PipeString
 @deprecate cholpfact(A)             cholfact(A, :U, pivot=true)
 @deprecate symmetrize!(A)      Base.LinAlg.copytri!(A, 'U')
 @deprecate symmetrize!(A, uplo)      Base.LinAlg.copytri!(A, uplo)
+@deprecate factorize!(A)       factorize(A)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
diff --git a/base/exports.jl b/base/exports.jl
index 525e7477e63db..34be8e5a09000 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -607,7 +607,6 @@ export
     eigvecs,
     expm,
     eye,
-    factorize!,
     factorize,
     givens,
     hessfact!,
diff --git a/base/linalg.jl b/base/linalg.jl
index 5335eb7bdb1e5..2e20b8cf43423 100644
--- a/base/linalg.jl
+++ b/base/linalg.jl
@@ -66,7 +66,6 @@ export
     sqrtm,
     eye,
     factorize,
-    factorize!,
     givens,
     gradient,
     hessfact,
diff --git a/base/linalg/bunchkaufman.jl b/base/linalg/bunchkaufman.jl
index 2638bb3b9934a..3a4fde9b95a95 100644
--- a/base/linalg/bunchkaufman.jl
+++ b/base/linalg/bunchkaufman.jl
@@ -2,7 +2,7 @@
 ## LD for BunchKaufman, UL for CholeskyDense, LU for LUDense and
 ## define size methods for Factorization types using it.
 
-type BunchKaufman{T<:BlasFloat} <: Factorization{T}
+immutable BunchKaufman{T} <: Factorization{T}
     LD::Matrix{T}
     ipiv::Vector{BlasInt}
     uplo::Char
@@ -18,10 +18,10 @@ function bkfact!{T<:BlasComplex}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric
     LD, ipiv = (symmetric ? LAPACK.sytrf! : LAPACK.hetrf!)(string(uplo)[1] , A)
     BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
 end
-bkfact!(A::StridedMatrix, args...) = bkfact!(float(A), args...)
-bkfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = bkfact!(copy(A), args...)
-bkfact(A::StridedMatrix, args...) = bkfact(float(A), args...)
+bkfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(copy(A), uplo, symmetric)
+bkfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(convert(Matrix{promote_type(Float32,typeof(sqrt(one(T))))},A),uplo,symmetric)
 
+convert{T}(::Type{BunchKaufman{T}},B::BunchKaufman) = BunchKaufman(convert(Matrix{T},B.LD),B.ipiv,B.uplo,B.symmetric)
 size(B::BunchKaufman) = size(B.LD)
 size(B::BunchKaufman,d::Integer) = size(B.LD,d)
 issym(B::BunchKaufman) = B.symmetric
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index bae8acee6bb98..8ed042afc42e3 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -287,29 +287,21 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
     end
 end
 
-function sqrtm{T<:Real}(A::StridedMatrix{T}, cond::Bool)
-    issym(A) && return sqrtm(Symmetric(A), cond)
-    
+function sqrtm{T<:Real}(A::StridedMatrix{T})
+    issym(A) && return sqrtm(Symmetric(A))
     n = chksquare(A)
-    SchurF = schurfact!(complex(A))
+    SchurF = schurfact(complex(A))
     R = full(sqrtm(Triangular(SchurF[:T])))
     retmat = SchurF[:vectors]*R*SchurF[:vectors]'
-    retmat2= all(imag(retmat) .== 0) ? real(retmat) : retmat
-    cond ? (retmat2, norm(R)^2/norm(SchurF[:T])) : retmat2
+    all(imag(retmat) .== 0) ? real(retmat) : retmat
 end
-function sqrtm{T<:Complex}(A::StridedMatrix{T}, cond::Bool)
-    ishermitian(A) && return sqrtm(Hermitian(A), cond)
-    
+function sqrtm{T<:Complex}(A::StridedMatrix{T})
+    ishermitian(A) && return sqrtm(Hermitian(A))
     n = chksquare(A)
     SchurF = schurfact(A)
     R = full(sqrtm(Triangular(SchurF[:T])))
-    retmat = SchurF[:vectors]*R*SchurF[:vectors]'
-    cond ? (retmat, norm(R)^2/norm(SchurF[:T])) : retmat
+    SchurF[:vectors]*R*SchurF[:vectors]'
 end
-
-sqrtm{T<:Integer}(A::StridedMatrix{T}, cond::Bool) = sqrtm(float(A), cond)
-sqrtm{T<:Integer}(A::StridedMatrix{Complex{T}}, cond::Bool) = sqrtm(complex128(A), cond)
-sqrtm(A::StridedMatrix) = sqrtm(A, false)
 sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
 sqrtm(a::Complex) = sqrt(a)
 
@@ -327,7 +319,7 @@ function inv(A::Matrix)
     return inv(lufact(A))
 end
 
-function factorize!{T}(A::Matrix{T})
+function factorize{T}(A::Matrix{T})
     m, n = size(A)
     if m == n
         if m == 1 return A[1] end
@@ -377,9 +369,9 @@ function factorize!{T}(A::Matrix{T})
             end
             if utri1
                 if (herm & (T <: Complex)) | sym
-                    return ldltd!(SymTridiagonal(diag(A), diag(A, -1)))
+                    return ldltd(SymTridiagonal(diag(A), diag(A, -1)))
                 end
-                return lufact!(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
+                return lufact(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
             end
         end
         if utri
@@ -388,32 +380,36 @@ function factorize!{T}(A::Matrix{T})
         if herm
             if T <: BlasFloat
                 C, info = LAPACK.potrf!('U', copy(A))
-            elseif typeof(one(T)/one(T)) <: BlasFloat
-                C, info = LAPACK.potrf!('U', float(A))
-            else
-                error("Unable to factorize hermitian $(typeof(A)). Try converting to other element type or use explicit factorization.")
+            else 
+                S = typeof(one(T)/one(T))
+                if S <: BlasFloat
+                    C, info = LAPACK.potrf!('U', convert(Matrix{S}, A))
+                else
+                    C, info = S <: Real ? LAPACK.potrf!('U', complex128(A)) : LAPACK.potrf!('U', complex128(A))
+                end
             end
             if info == 0 return Cholesky(C, 'U') end
-            return factorize!(Hermitian(A))
+            return factorize(Hermitian(A))
         end
         if sym
             if T <: BlasFloat
                 C, info = LAPACK.potrf!('U', copy(A))
-            elseif eltype(one(T)/one(T)) <: BlasFloat
-                C, info = LAPACK.potrf!('U', float(A))
             else
-                error("Unable to factorize symmetric $(typeof(A)). Try converting to other element type or use explicit factorization.")
+                S = eltype(one(T)/one(T))
+                if S <: BlasFloat
+                    C, info = LAPACK.potrf!('U', convert(Matrix{S},A))
+                else
+                    C, info = S <: Real ? LAPACK.potrf!('U', float64(A)) : LAPACK.potrf!('U', complex(A))
+                end
             end
             if info == 0 return Cholesky(C, 'U') end
-            return factorize!(Symmetric(A))
+            return factorize(Symmetric(A))
         end
-        return lufact!(A)
+        return lufact(A)
     end
-    return qrfact!(A,pivot=true)
+    qrfact(A,pivot=T<:BlasFloat)
 end
 
-factorize(A::AbstractMatrix) = factorize!(copy(A))
-
 (\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)
 function (\)(A::StridedMatrix, B::StridedVecOrMat)
     m, n = size(A)
@@ -424,32 +420,31 @@ function (\)(A::StridedMatrix, B::StridedVecOrMat)
         istriu(A) && return \(Triangular(A, :U),B)
         return \(lufact(A),B)
     end
-    return qrfact(A,pivot=true)\B
+    return qrfact(A,pivot=eltype(A)<:BlasFloat)\B
 end
 
 ## Moore-Penrose inverse
-function pinv{T<:BlasFloat}(A::StridedMatrix{T})
-    m, n = size(A)
-    (m == 0 || n == 0) && return Array(T, n, m)
-    SVD         = svdfact(A, true)
-    Sinv        = zeros(T, length(SVD[:S]))
-    index       = SVD[:S] .> eps(real(one(T)))*max(m,n)*maximum(SVD[:S])
-    Sinv[index] = 1.0 ./ SVD[:S][index]
+function pinv{T}(A::StridedMatrix{T})
+    SVD         = svdfact(A, thin=true)
+    S           = eltype(SVD[:S])
+    m, n        = size(A)
+    (m == 0 || n == 0) && return Array(S, n, m)
+    Sinv        = zeros(S, length(SVD[:S]))
+    index       = SVD[:S] .> eps(real(float(one(T))))*max(m,n)*maximum(SVD[:S])
+    Sinv[index] = one(S) ./ SVD[:S][index]
     return SVD[:Vt]'scale(Sinv, SVD[:U]')
 end
-pinv{T<:Integer}(A::StridedMatrix{T}) = pinv(float(A))
 pinv(a::StridedVector) = pinv(reshape(a, length(a), 1))
 pinv(x::Number) = one(x)/x
 
 ## Basis for null space
-function null{T<:BlasFloat}(A::StridedMatrix{T})
+function null{T}(A::StridedMatrix{T})
     m, n = size(A)
     (m == 0 || n == 0) && return eye(T, n)
-    SVD = svdfact(A, false)
+    SVD = svdfact(A, thin=false)
     indstart = sum(SVD[:S] .> max(m,n)*maximum(SVD[:S])*eps(eltype(SVD[:S]))) + 1
     return SVD[:V][:,indstart:end]
 end
-null{T<:Integer}(A::StridedMatrix{T}) = null(float(A))
 null(a::StridedVector) = null(reshape(a, length(a), 1))
 
 function cond(A::StridedMatrix, p::Real=2) 
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 3eafd5d579573..432e65270f2e8 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -13,11 +13,11 @@ end
 ##########################
 # Cholesky Factorization #
 ##########################
-type Cholesky{T<:BlasFloat} <: Factorization{T}
+immutable Cholesky{T} <: Factorization{T}
     UL::Matrix{T}
     uplo::Char
 end
-type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
+immutable CholeskyPivoted{T} <: Factorization{T}
     UL::Matrix{T}
     uplo::Char
     piv::Vector{BlasInt}
@@ -26,7 +26,7 @@ type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
     info::BlasInt
 end
 
-function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0)
+function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
     uplochar = string(uplo)[1]
     if pivot
         A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
@@ -36,14 +36,16 @@ function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=fal
         return @assertposdef Cholesky(C, uplochar) info
     end
 end
-cholfact!(A::StridedMatrix, args...) = cholfact!(float(A), args...)
-cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol)
-cholfact(A::StridedMatrix, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(float(A), uplo, pivot=pivot, tol=tol)
+cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol)
+cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? cholfact!(convert(Matrix{S},A), uplo, pivot=pivot, tol=tol) : cholfact!(copy(A), uplo, pivot=pivot, tol=tol)) # When julia Cholesky has been implemented, the promotion should be changed.
 cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x) == 0 && real(x) > 0)
 
 chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo]
 chol(A::Union(Number, AbstractMatrix)) = triu!(cholfact(A, :U).UL)
 
+convert{T,S}(::Type{Cholesky{T}},C::Cholesky{S}) = Cholesky(convert(Matrix{T},C.UL),C.uplo)
+convert{T,S}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted{S}) = CholeskyPivoted(convert(Matrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info)
+
 size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
 size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)
 
@@ -71,12 +73,12 @@ end
 show(io::IO, C::Cholesky) = (println("$(typeof(C)) with factor:");show(io,C[symbol(C.uplo)]))
 
 A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B)
+A_ldiv_B!(C::Cholesky, B::StridedVecOrMat) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B)) : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B))
 
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T})
     chkfullrank(C)
     ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv)
 end
-
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
     chkfullrank(C)
     n = size(C, 1)
@@ -89,6 +91,8 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
     end
     B
 end
+A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)]
+A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:]
 
 function det{T}(C::Cholesky{T})
     dd = one(T)
@@ -127,7 +131,6 @@ end
 
 lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(LAPACK.getrf!(A)...)
 function lufact!{T}(A::AbstractMatrix{T})
-    typeof(one(T)/one(T)) <: BlasFloat && return lufact!(float(A))
     m, n = size(A)
     minmn = min(m,n)
     info = 0
@@ -135,7 +138,7 @@ function lufact!{T}(A::AbstractMatrix{T})
     for k = 1:minmn
         # find index max
         kp = 1
-        amax = zero(T)
+        amax = real(zero(T))
         for i = k:m
             absi = abs(A[i,k])
             if absi > amax
@@ -168,7 +171,8 @@ function lufact!{T}(A::AbstractMatrix{T})
     if minmn > 0 && A[minmn,minmn] == 0; info = minmn; end
     LU(A, ipiv, convert(BlasInt, info))
 end
-lufact(A::StridedMatrix) = lufact!(copy(A))
+lufact{T<:BlasFloat}(A::StridedMatrix{T}) = lufact!(copy(A))
+lufact{T}(A::StridedMatrix{T}) = (S = typeof(one(T)/one(T)); S != T ? lufact!(convert(Matrix{S}, A)) : lufact!(copy(A)))
 lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt))
 
 function lu(A::Union(Number, AbstractMatrix))
@@ -250,24 +254,48 @@ cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:
 # QR Factorization #
 ####################
 
-# Note. For QR factorization without pivoting, the WY representation based method introduced in LAPACK 3.4
-type QR{S<:BlasFloat} <: Factorization{S}
+immutable QR{T} <: Factorization{T}
+    factors::Matrix{T}
+    τ::Vector{T}
+end
+# Note. For QRCompactWY factorization without pivoting, the WY representation based method introduced in LAPACK 3.4
+immutable QRCompactWY{S} <: Factorization{S}
     factors::Matrix{S}
     T::Matrix{S}
 end
-QR{T<:BlasFloat}(A::StridedMatrix{T}, nb::Integer = min(minimum(size(A)), 36)) = QR(LAPACK.geqrt!(A, nb)...)
 
-type QRPivoted{T} <: Factorization{T}
+immutable QRPivoted{T} <: Factorization{T}
     factors::Matrix{T}
-    tau::Vector{T}
+    τ::Vector{T}
     jpvt::Vector{BlasInt}
 end
 
-qrfact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = pivot ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QR(A)
-qrfact!(A::StridedMatrix; pivot=false) = qrfact!(float(A), pivot=pivot)
-qrfact{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = qrfact!(copy(A), pivot=pivot)
-qrfact(A::StridedMatrix; pivot=false) = qrfact!(float(A), pivot=pivot)
-qrfact(x::Integer) = qrfact(float(x))
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = pivot ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QRCompactWY(LAPACK.geqrt!(A, min(minimum(size(A)), 36))...)
+function qrfact!{T}(A::AbstractMatrix{T}; pivot=false)
+    pivot && warn("pivoting only implemented for Float32, Float64, Complex64 and Complex128")
+    m, n = size(A)
+    τ = zeros(T, min(m,n))
+    @inbounds begin
+        for k = 1:min(m-1+(T<:Complex),n)
+            τk = elementaryLeft!(A, k, k)
+            τ[k] = τk
+            for j = k+1:n
+                νAj = A[k,j]
+                for i = k+1:m
+                    νAj += conj(A[i,k])*A[i,j]
+                end
+                νAj *= conj(τk)
+                A[k,j] -= νAj
+                for i = k+1:m
+                    A[i,j] -= A[i,k]*νAj
+                end
+            end
+        end
+    end
+    QR(A, τ)
+end
+qrfact{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = qrfact!(copy(A),pivot=pivot)
+qrfact{T}(A::StridedMatrix{T}; pivot=false) = (S = typeof(sqrt(one(T)));S != T ? qrfact!(convert(Matrix{S},A), pivot=pivot) : qrfact!(copy(A),pivot=pivot))
 qrfact(x::Number) = QR(fill(one(x), 1, 1), fill(x, 1, 1))
 
 function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true)
@@ -275,58 +303,23 @@ function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true)
     full(F[:Q], thin=thin), F[:R]
 end
 
+convert{T}(::Type{QR{T}},A::QR) = QR(convert(Matrix{T}, A.factors), convert(Vector{T}, A.τ))
+convert{T}(::Type{QRCompactWY{T}},A::QRCompactWY) = QRCompactWY(convert(Matrix{T}, A.factors), convert(Matrix{T}, A.T))
+convert{T}(::Type{QRPivoted{T}},A::QRPivoted) = QRPivoted(convert(Matrix{T}, A.factors), convert(Vector{T}, A.τ), A.jpvt)
+
 function getindex(A::QR, d::Symbol)
     d == :R && return triu(A.factors[1:minimum(size(A)),:])
-    d == :Q && return QRPackedQ(A)
+    d == :Q && return QRPackedQ(A.factors,A.τ)
     throw(KeyError(d))
 end
-
-type QRPackedQ{S} <: AbstractMatrix{S} 
-    factors::Matrix{S}                      
-    T::Matrix{S}                       
-end
-QRPackedQ(A::QR) = QRPackedQ(A.factors, A.T)
-type QRPivotedQ{T} <: AbstractMatrix{T}
-    factors::Matrix{T}                       # Householder transformations and R
-    tau::Vector{T}                      # Scalar factors of transformations
-end
-QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.factors, A.tau)
-
-size(A::Union(QR, QRPivoted), dim::Integer) = size(A.factors, dim)
-size(A::Union(QR, QRPivoted)) = size(A.factors)
-size(A::Union(QRPackedQ, QRPivotedQ), dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError())
-size(A::Union(QRPackedQ, QRPivotedQ)) = size(A, 1), size(A, 2)
-
-full{T<:BlasFloat}(A::QRPackedQ{T}; thin::Bool=true) = A * (thin ? eye(T, size(A.factors)...) : eye(T, size(A.factors,1)))
-function full{T<:BlasFloat}(A::QRPivotedQ{T}; thin::Bool=true)
-    m, n = size(A.factors)
-    Ahhpad = thin ? copy(A.factors) : [A.factors zeros(T, m, max(0, m - n))]
-    LAPACK.orgqr!(Ahhpad, A.tau)
-end
-
-print_matrix(io::IO, A::QRPackedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols)
-print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols)
-
-## Multiplication by Q from the QR decomposition
-function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
-    Bpad = size(B, 1)==size(A.factors, 2) ? (isa(B, Vector) ? [B; zeros(T, size(A.factors, 1) - size(A.factors, 2))] : [B; zeros(T, size(A.factors, 1) - size(A.factors, 2), size(B, 2))]) : copy(B)
-    LAPACK.gemqrt!('L', 'N', A.factors, A.T, Bpad)
-end
-*{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.factors, B.T, copy(A))
-Ac_mul_B{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.factors,A.T,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.factors,A.T,copy(B))
-Ac_mul_B(A::QRPackedQ, B::StridedVecOrMat) = Ac_mul_B(A, float(B))
-function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
-    Apad = size(A, 2)==size(B.factors, 2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : copy(A)     
-    LAPACK.gemqrt!('R', iseltype(B.factors,Complex) ? 'C' : 'T', B.factors, B.T, Apad)
+function getindex(A::QRCompactWY, d::Symbol)
+    d == :R && return triu(A.factors[1:minimum(size(A)),:])
+    d == :Q && return QRCompactWYQ(A.factors,A.T)
+    throw(KeyError(d))
 end
-## Least squares solution.  Should be more careful about cases with m < n
-\(A::QR, B::StridedVector) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2)]
-\(A::QR, B::StridedMatrix) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2),:]
-
-function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
+function getindex{T}(A::QRPivoted{T}, d::Symbol)
     d == :R && return triu(A.factors[1:minimum(size(A)),:])
-    d == :Q && return QRPivotedQ(A)
+    d == :Q && return QRPackedQ(A.factors,A.τ)
     d == :p && return A.jpvt
     if d == :P
         p = A[:p]
@@ -340,9 +333,161 @@ function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
     throw(KeyError(d))
 end
 
+immutable QRPackedQ{T} <: AbstractMatrix{T}
+    factors::Matrix{T}
+    τ::Vector{T}
+end
+immutable QRCompactWYQ{S} <: AbstractMatrix{S} 
+    factors::Matrix{S}                      
+    T::Matrix{S}                       
+end
+
+convert{T,S}(::Type{QRPackedQ{T}}, Q::QRPackedQ{S}) = QRPackedQ(convert(Matrix{T}, Q.factors), convert(Vector{T}, Q.τ))
+convert{S1,S2}(::Type{QRCompactWYQ{S1}}, Q::QRCompactWYQ{S2}) = QRCompactWYQ(convert(Matrix{S1}, Q.factors), convert(Matrix{S1}, Q.T))
+
+size(A::Union(QR,QRCompactWY,QRPivoted), dim::Integer) = size(A.factors, dim)
+size(A::Union(QR,QRCompactWY,QRPivoted)) = size(A.factors)
+size(A::Union(QRPackedQ,QRCompactWYQ), dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError())
+size(A::Union(QRPackedQ,QRCompactWYQ)) = size(A, 1), size(A, 2)
+
+full{T}(A::Union(QRPackedQ{T},QRCompactWYQ{T}); thin::Bool=true) = A_mul_B!(A, thin ? eye(T, size(A.factors)...) : eye(T, size(A.factors,1)))
+
+print_matrix(io::IO, A::Union(QRPackedQ,QRCompactWYQ), rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols)
+
+## Multiplication by Q
+### QB
+A_mul_B!{T<:BlasFloat}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','N',A.factors,A.T,B)
+A_mul_B!{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','N',A.factors,A.τ,B)
+function A_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
+    mA, nA = size(A.factors)
+    mB, nB = size(B,1), size(B,2)
+    mA == mB || throw(DimensionMismatch(""))
+    Afactors = A.factors
+    @inbounds begin
+        for k = min(mA,nA):-1:1
+            for j = 1:nB
+                νBj = B[k,j]
+                for i = k+1:mB
+                    νBj += conj(Afactors[i,k])*B[i,j]
+                end
+                νBj *= conj(A.τ[k])
+                B[k,j] -= νBj
+                for i = k+1:mB
+                    B[i,j] -= Afactors[i,k]*νBj
+                end
+            end
+        end
+    end
+    B
+end
+function *{TA,Tb}(A::Union(QRPackedQ{TA},QRCompactWYQ{TA}), b::StridedVector{Tb})
+    TAb = promote_type(TA,Tb)
+    Anew = convert(typeof(A).name.primary{TAb},A)
+    bnew = size(A.factors,1) == length(b) ? (Tb == TAb ? copy(b) : convert(Vector{TAb}, b)) : (size(A.factors,2) == length(b) ? [b,zeros(TAb, size(A.factors,1)-length(b))] : throw(DimensionMismatch("")))
+    A_mul_B!(Anew,bnew)
+end
+function *{TA,TB}(A::Union(QRPackedQ{TA},QRCompactWYQ{TA}), B::StridedMatrix{TB})
+    TAB = promote_type(TA,TB)
+    Anew = convert(typeof(A).name.primary{TAB},A)
+    Bnew = size(A.factors,1) == size(B,1) ? (TB == TAB ? copy(B) : convert(Matrix{TAB}, B)) : (size(A.factors,2) == size(B,1) ? [B;zeros(TAB, size(A.factors,1)-size(B,1),size(B,2))] : throw(DimensionMismatch("")))
+    A_mul_B!(Anew,Bnew)
+end
+### QcB
+Ac_mul_B!{T<:BlasReal}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.factors,A.T,B)
+Ac_mul_B!{T<:BlasComplex}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.factors,A.T,B)
+Ac_mul_B!{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.factors,A.τ,B)
+Ac_mul_B!{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.factors,A.τ,B)
+function Ac_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
+    mA, nA = size(A.factors)
+    mB, nB = size(B,1), size(B,2)
+    mA == mB || throw(DimensionMismatch(""))
+    Afactors = A.factors
+    @inbounds begin
+        for k = 1:min(mA,nA)
+            for j = 1:nB
+                νBj = B[k,j]
+                for i = k+1:mB
+                    νBj += conj(Afactors[i,k])*B[i,j]
+                end
+                νBj *= A.τ[k]
+                B[k,j] -= νBj
+                for i = k+1:mB
+                    B[i,j] -= Afactors[i,k]*νBj
+                end
+            end
+        end
+    end
+    B
+end
+function Ac_mul_B{TQ<:Number,TB<:Number}(Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}), B::StridedVecOrMat{TB})
+    TQB = promote_type(TQ,TB)
+    Ac_mul_B!(convert(typeof(Q).name.primary{TQB}, Q), TB == TQB ? copy(B) : convert(typeof(B).name.primary{TQB}, B))
+end
+### AQ
+A_mul_B!{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','N', B.factors, B.T, A)
+A_mul_B!(A::StridedVecOrMat, B::QRPackedQ) = LAPACK.ormqr!('R', 'N', B.factors, B.τ, A)
+function A_mul_B!{T}(A::StridedMatrix{T},Q::QRPackedQ{T})
+    mQ, nQ = size(Q.factors)
+    mA, nA = size(A,1), size(A,2)
+    nA == mQ || throw(DimensionMismatch(""))
+    Qfactors = Q.factors
+    @inbounds begin
+        for k = 1:min(mQ,nQ)
+            for i = 1:mA
+                νAi = A[i,k]
+                for j = k+1:mQ
+                    νAi += Qfactors[j,k]*A[i,j]
+                end
+                νAi *= conj(Q.τ[k])
+                A[i,k] -= νAi
+                for j = k+1:nA
+                    A[i,j] -= Qfactors[j,k]*νAi
+                end
+            end
+        end
+    end
+    A
+end
+function *{TA,TQ}(A::StridedVecOrMat{TA}, Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}))
+    TAQ = promote_type(TA, TQ)
+    A_mul_B!(TA==TAQ ? copy(A) : convert(typeof(A).name.primary{TAQ}, A), convert(typeof(Q).name.primary{TAQ}, Q))
+end
+### AQc
+A_mul_Bc!{T<:BlasReal}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','T',B.factors,B.T,A)
+A_mul_Bc!{T<:BlasComplex}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','C',B.factors,B.T,A)
+A_mul_Bc!{T<:BlasReal}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.ormqr!('R','T',B.factors,B.τ,A)
+A_mul_Bc!{T<:BlasComplex}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.ormqr!('R','C',B.factors,B.τ,A)
+function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T})
+    mQ, nQ = size(Q.factors)
+    mA, nA = size(A,1), size(A,2)
+    nA == mQ || throw(DimensionMismatch(""))
+    Qfactors = Q.factors
+    @inbounds begin
+        for k = min(mQ,nQ):-1:1
+            for i = 1:mA
+                νAi = A[i,k]
+                for j = k+1:mQ
+                    νAi += Qfactors[j,k]*A[i,j]
+                end
+                νAi *= Q.τ[k]
+                A[i,k] -= νAi
+                for j = k+1:nA
+                    A[i,j] -= Qfactors[j,k]*νAi
+                end
+            end
+        end
+    end
+    A
+end
+function A_mul_Bc{TA,TB}(A::AbstractVecOrMat{TA}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB}))
+    TAB = promote_type(TA,TB)
+    A_mul_Bc!(size(A,2)==size(B.factors,1) ? copy(A) : (size(A,2)==size(B.factors,2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : throw(DimensionMismatch(""))))
+end
+
 # Julia implementation similarly to xgelsy
-function \{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
-    nr = minimum(size(A.factors))
+function A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedMatrix{T}, rcond::Real)
+    mA, nA = size(A.factors)
+    nr = min(mA,nA)
     nrhs = size(B, 2)
     if nr == 0 return zeros(0, nrhs), 0 end
     ar = abs(A.factors[1])
@@ -360,27 +505,79 @@ function \{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
         rnk += 1
         # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end
     end
-    C, tau = LAPACK.tzrzf!(A.factors[1:rnk,:])
-    X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.factors, 2) - rnk, nrhs)]
-    LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, tau, X)
-    return X[invperm(A[:p]),:], rnk
+    C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:])
+    A_ldiv_B!(Triangular(C[1:rnk,1:rnk],:U),sub(Ac_mul_B!(A[:Q],sub(B, 1:mA, 1:nrhs)),1:rnk,1:nrhs))
+    B[rnk+1:end,:] = zero(T)
+    LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, τ, sub(B,1:nA,1:nrhs))
+    return isa(A,QRPivoted) ? B[invperm(A[:p]),:] : B[1:nA,:], rnk
+end
+A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedVector, rcond::Real) = A_ldiv_B!(A,reshape(B,length(B),1), rcond)
+A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedVecOrMat{T}) = A_ldiv_B!(A, B, sqrt(eps(real(float(one(eltype(B)))))))[1]
+# A_ldiv_B!(A::QRCompactWY, B::StridedVector) = A_ldiv_B!(Triangular(A[:R], :U), Ac_mul_B(A[:Q], B))[1:size(A, 2)]
+# A_ldiv_B!(A::QRCompactWY, B::StridedMatrix) = A_ldiv_B!(Triangular(A[:R], :U), Ac_mul_B(A[:Q], B))[1:size(A, 2),:]
+function A_ldiv_B!{T}(A::QR{T},B::StridedMatrix{T})
+    m, n = size(A)
+    minmn = min(m,n)
+    mB, nB = size(B)
+    Ac_mul_B!(A[:Q],sub(B,1:m,1:nB)) # Reconsider when arrayviews are merged.
+    R = A[:R]
+    @inbounds begin
+        if n > m # minimum norm solution
+            τ = zeros(T,m)
+            for k = m:-1:1 # Trapezoid to triangular by elementary operation
+                τ[k] = elementaryRightTrapezoid!(R,k)
+                for i = 1:k-1
+                    νRi = R[i,k]
+                    for j = m+1:n
+                        νRi += R[i,j]*R[k,j]
+                    end
+                    νRi *= τ[k]
+                    R[i,k] -= νRi
+                    for j = m+1:n
+                        R[i,j] -= νRi*R[k,j]
+                    end
+                end
+            end
+        end
+        for k = 1:nB # solve triangular system. When array views are implemented, consider exporting    to function.
+            for i = minmn:-1:1
+                for j = i+1:minmn
+                    B[i,k] -= R[i,j]*B[j,k]
+                end
+                B[i,k] /= R[i,i]
+            end
+        end
+        if n > m # Apply elemenary transformation to solution
+            B[m+1:mB,1:nB] = zero(T)
+            for j = 1:nB
+                for k = 1:m
+                    νBj = B[k,j]
+                    for i = m+1:n
+                        νBj += B[i,j]*conj(R[k,i])
+                    end
+                    νBj *= τ[k]
+                    B[k,j] -= νBj
+                    for i = m+1:n
+                        B[i,j] -= R[k,i]*νBj
+                    end
+                end
+            end
+        end
+    end
+    return B[1:n,:]
 end
-\(A::QRPivoted, B::StridedMatrix, rcond::Real) = \(A, float(B), rcond)
-\(A::QRPivoted, B::StridedMatrix) = (\)(A, B, sqrt(eps(typeof(float(real(B[1]))))))[1]
-\(A::QRPivoted, B::StridedVector) = (\)(A, reshape(B, length(B), 1))[:]
-
-## Multiplication by Q from the Pivoted QR decomposition
-function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
-    Bpad = size(A.factors, 2)==size(B, 1) ? (isa(B, Vector) ? [B; zeros(T, size(A.factors, 1) - size(A.factors, 2))] : [B; zeros(T, size(A.factors, 1) - size(A.factors, 2), size(B, 2))]) : copy(B)
-    LAPACK.ormqr!('L', 'N', A.factors, A.tau, Bpad)
+A_ldiv_B!(A::QR, B::StridedVector) = A_ldiv_B!(A, reshape(B, length(B), 1))[:]
+A_ldiv_B!(A::QRPivoted, B::StridedVector) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt)]
+A_ldiv_B!(A::QRPivoted, B::StridedMatrix) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt),:]
+function \{TA,TB}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),B::StridedVector{TB})
+    S = promote_type(TA,TB)
+    m,n = size(A)
+    n > m ? A_ldiv_B!(convert(typeof(A).name.primary{S},A),[B,zeros(S,n-m)]) : A_ldiv_B!(convert(typeof(A).name.primary{S},A), S == TB ? copy(B) : convert(Vector{S}, B))
 end
-Ac_mul_B{T<:BlasReal}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.factors,A.tau,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.factors,A.tau,copy(B))
-Ac_mul_B(A::QRPivotedQ, B::StridedVecOrMat) = Ac_mul_B(A, float(B))
-*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.factors, B.tau, copy(A))
-function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
-    Apad = size(A, 2)==size(B.factors, 2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : copy(A)
-    LAPACK.ormqr!('R', iseltype(B.factors,Complex) ? 'C' : 'T', B.factors, B.tau, Apad)
+function \{TA,TB}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),B::StridedMatrix{TB})
+    S = promote_type(TA,TB)
+    m,n = size(A)
+    n > m ? A_ldiv_B!(convert(typeof(A).name.primary{S},A),[B;zeros(S,n-m,size(B,2))]) : A_ldiv_B!(convert(typeof(A).name.primary{S},A), S == TB ? copy(B) : convert(Matrix{S}, B))
 end
 
 ##TODO:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
@@ -388,27 +585,26 @@ end
 ## Lower priority: Add LQ, QL and RQ factorizations
 
 # FIXME! Should add balancing option through xgebal
-type Hessenberg{T} <: Factorization{T}
+immutable Hessenberg{T} <: Factorization{T}
     factors::Matrix{T}
-    tau::Vector{T}
-    function Hessenberg(hh::Matrix{T}, tau::Vector{T})
+    τ::Vector{T}
+    function Hessenberg(hh::Matrix{T}, τ::Vector{T})
         chksquare(hh)
-        new(hh, tau)
+        new(hh, τ)
     end
 end
-Hessenberg{T<:BlasFloat}(hh::Matrix{T}, tau::Vector{T}) = Hessenberg{T}(hh, tau)
+Hessenberg{T<:BlasFloat}(hh::Matrix{T}, τ::Vector{T}) = Hessenberg{T}(hh, τ)
 Hessenberg(A::StridedMatrix) = Hessenberg(LAPACK.gehrd!(A)...)
 
 hessfact!{T<:BlasFloat}(A::StridedMatrix{T}) = Hessenberg(A)
-hessfact!(A::StridedMatrix) = hessfact!(float(A))
 hessfact{T<:BlasFloat}(A::StridedMatrix{T}) = hessfact!(copy(A))
-hessfact(A::StridedMatrix) = hessfact!(float(A))
+hessfact{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? hessfact!(convert(Matrix{S},A)) : hessfact!(copy(A)))
 
-type HessenbergQ{T} <: AbstractMatrix{T}
+immutable HessenbergQ{T} <: AbstractMatrix{T}
     hh::Matrix{T}
-    tau::Vector{T}
+    τ::Vector{T}
 end
-HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.tau)
+HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.τ)
 size(A::HessenbergQ, args...) = size(A.factors, args...)
 getindex(A::HessenbergQ, i::Real) = getindex(full(A), i)
 getindex(A::HessenbergQ, i::AbstractArray) = getindex(full(A), i)
@@ -420,20 +616,20 @@ function getindex(A::Hessenberg, d::Symbol)
     throw(KeyError(d))
 end
 
-full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.tau)
+full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.τ)
 
 #######################
 # Eigendecompositions #
 #######################
 
 # Eigenvalues
-type Eigen{T,V} <: Factorization{T}
+immutable Eigen{T,V} <: Factorization{T}
     values::Vector{V}
     vectors::Matrix{T}
 end
 
 # Generalized eigenvalue problem.
-type GeneralizedEigen{T,V} <: Factorization{T}
+immutable GeneralizedEigen{T,V} <: Factorization{T}
     values::Vector{V}
     vectors::Matrix{T}
 end
@@ -450,7 +646,6 @@ function eigfact!{T<:BlasReal}(A::StridedMatrix{T}; balance::Symbol=:balance)
     n = size(A, 2)
     n==0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
     issym(A) && return eigfact!(Symmetric(A))
-
     A, WR, WI, VL, VR, _ = LAPACK.geevx!(balance == :balance ? 'B' : (balance == :diagonal ? 'S' : (balance == :permute ? 'P' : (balance == :nobalance ? 'N' : throw(ArgumentError("balance must be either ':balance', ':diagonal', ':permute' or ':nobalance'"))))), 'N', 'V', 'N', A)
     all(WI .== 0.) && return Eigen(WR, VR)
     evec = zeros(Complex{T}, n, n)
@@ -474,9 +669,8 @@ function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}; balance::Symbol=:balance)
     ishermitian(A) && return eigfact!(Hermitian(A)) 
     return Eigen(LAPACK.geevx!(balance == :balance ? 'B' : (balance == :diagonal ? 'S' : (balance == :permute ? 'P' : (balance == :nobalance ? 'N' : throw(ArgumentError("balance must be either ':balance', ':diagonal', ':permute' or ':nobalance'"))))), 'N', 'V', 'N', A)[[2,4]]...)
 end
-eigfact!(A::StridedMatrix; balance::Symbol=:balance) = eigfact!(float(A), balance=balance)
-eigfact{T<:BlasFloat}(x::StridedMatrix{T}; balance::Symbol=:balance) = eigfact!(copy(x), balance=balance)
-eigfact(A::StridedMatrix; balance::Symbol=:balance) = eigfact!(float(A), balance=balance)
+eigfact{T<:BlasFloat}(A::StridedMatrix{T}; balance::Symbol=:balance) = eigfact!(copy(A), balance=balance)
+eigfact{T}(A::StridedMatrix{T}; balance::Symbol=:balance) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? eigfact!(convert(Matrix{S}, A), balance=balance) : eigfact!(copy(A), balance=balance))
 eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1))
 
 function eig(A::Union(Number, AbstractMatrix), args...)
@@ -496,9 +690,8 @@ function eigvals!{T<:BlasComplex}(A::StridedMatrix{T})
     ishermitian(A) && return eigvals(Hermitian(A))
     return LAPACK.geev!('N', 'N', A)[1]
 end
-eigvals!(A::AbstractMatrix, args...) = eigvals!(float(A), args...)
-eigvals{T<:BlasFloat}(A::AbstractMatrix{T}) = eigvals!(copy(A))
-eigvals(A::AbstractMatrix, args...) = eigvals!(float(A), args...)
+eigvals{T<:BlasFloat}(A::StridedMatrix{T}) = eigvals!(copy(A))
+eigvals{T}(A::AbstractMatrix{T}) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? eigvals!(convert(Matrix{S}, A)) : eigvals!(copy(A)))
 
 eigvals(x::Number) = [one(x)]
 
@@ -542,9 +735,8 @@ function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
     alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
     return GeneralizedEigen(alpha./beta, vr)
 end
-eigfact!(A::StridedMatrix, B::StridedMatrix) = eigfact!(float(A), float(B))
-eigfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = eigfact!(copy(A), copy(B))
-eigfact(A::StridedMatrix, B::StridedMatrix) = eigfact!(float(A), float(B))
+eigfact{T<:BlasFloat}(A::AbstractMatrix{T}, B::StridedMatrix{T}) = eigfact!(copy(A),copy(B))
+eigfact{TA,TB}(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) = (S = promote_type(Float32,typeof(sqrt(one(TA))),TB); eigfact!(S != TA ? convert(Matrix{S},A) : copy(A), S != TB ? convert(Matrix{S},B) : copy(B)))
 
 function eig(A::AbstractMatrix, B::AbstractMatrix)
     F = eigfact(A, B)
@@ -561,17 +753,16 @@ function eigvals!{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
     alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
     alpha./beta
 end
-eigvals!(A::AbstractMatrix, B::AbstractMatrix) = eigvals!(float(A), float(B))
-eigvals{T<:BlasFloat}(A::AbstractMatrix{T}, B::AbstractMatrix{T}) = eigvals!(copy(A), copy(B))
-eigvals(A::AbstractMatrix, B::AbstractMatrix) = eigvals!(float(A), float(B))
+eigvals{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = eigvals!(copy(A),copy(B))
+eigvals{TA,TB}(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) = (S = promote_type(Float32,typeof(sqrt(one(TA))),TB); eigvals!(S != TA ? convert(Matrix{S},A) : copy(A), S != TB ? convert(Matrix{S},B) : copy(B)))
 
 # SVD
-type SVD{T<:BlasFloat,Tr} <: Factorization{T}
+immutable SVD{T<:BlasFloat,Tr} <: Factorization{T}
     U::Matrix{T}
     S::Vector{Tr}
     Vt::Matrix{T}
 end
-function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, thin::Bool)
+function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}; thin::Bool=true)
     m,n = size(A)
     if m == 0 || n == 0
         u,s,vt = (eye(T, m, thin ? n : m), real(zeros(T,0)), eye(T,n,n))
@@ -580,15 +771,13 @@ function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, thin::Bool)
     end
     SVD(u,s,vt)
 end
-svdfact!(A::StridedVecOrMat, args...) = svdfact!(float(A), args...)
-svdfact!{T<:BlasFloat}(a::Vector{T}, thin::Bool) = svdfact!(reshape(a, length(a), 1), thin)
-svdfact!{T<:BlasFloat}(A::StridedVecOrMat{T}) = svdfact!(A, true)
-svdfact(A::StridedVecOrMat, args...) = svdfact!(eltype(A)<:BlasFloat ? copy(A) : float(A), args...)
-svdfact(x::Number, thin::Bool=true) = SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
-svdfact(x::Integer, thin::Bool=true) = svdfact(float(x), thin)
+svdfact{T<:BlasFloat}(A::StridedMatrix{T};thin=true) = svdfact!(copy(A),thin=thin)
+svdfact{T}(A::StridedVecOrMat{T};thin=true) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? svdfact!(convert(Matrix{S},A),thin=thin) : svdfact!(copy(A),thin=thin))
+svdfact(x::Number; thin::Bool=true) = SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
+svdfact(x::Integer; thin::Bool=true) = svdfact(float(x), thin=thin)
 
-function svd(A::Union(Number, AbstractArray), thin::Bool=true)
-    F = svdfact(A, thin)
+function svd(A::Union(Number, AbstractArray); thin::Bool=true)
+    F = svdfact(A, thin=thin)
     F.U, F.S, F.Vt'
 end
 
@@ -601,9 +790,8 @@ function getindex(F::SVD, d::Symbol)
 end
 
 svdvals!{T<:BlasFloat}(A::StridedMatrix{T}) = any([size(A)...].==0) ? zeros(T, 0) : LAPACK.gesdd!('N', A)[2]
-svdvals!(A::StridedMatrix) = svdvals!(float(A))
 svdvals{T<:BlasFloat}(A::StridedMatrix{T}) = svdvals!(copy(A))
-svdvals(A::StridedMatrix) = svdvals!(float(A))
+svdvals{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? svdvals!(convert(Matrix{S}, A)) : svdvals!(copy(A)))
 svdvals(x::Number) = [abs(x)]
 
 # SVD least squares
@@ -615,7 +803,7 @@ function \{T<:BlasFloat}(A::SVD{T}, B::StridedVecOrMat{T})
 end
 
 # Generalized svd
-type GeneralizedSVD{T} <: Factorization{T}
+immutable GeneralizedSVD{T} <: Factorization{T}
     U::Matrix{T}
     V::Matrix{T}
     Q::Matrix{T}
@@ -630,9 +818,8 @@ function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
     U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', A, B)
     GeneralizedSVD(U, V, Q, a, b, int(k), int(l), R)
 end
-svdfact!(A::StridedMatrix, B::StridedMatrix) = svdfact!(float(A), float(B))
-svdfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = svdfact!(copy(A), copy(B))
-svdfact(A::StridedMatrix, B::StridedMatrix) = svdfact!(float(A), float(B))
+svdfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = svdfact!(copy(A),copy(B))
+svdfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(sqrt(one(TA))),TB); svdfact!(S != TA ? convert(Matrix{S},A) : copy(A), S != TB ? convert(Matrix{S},B) : copy(B)))
 
 function svd(A::AbstractMatrix, B::AbstractMatrix)
     F = svdfact(A, B)
@@ -675,20 +862,18 @@ function svdvals!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
     _, _, _, a, b, k, l, _ = LAPACK.ggsvd!('N', 'N', 'N', A, B)
     a[1:k + l] ./ b[1:k + l]
 end
-svdvals!(A::StridedMatrix, B::StridedMatrix) = svdvals!(float(A), float(B))
-svdvals{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = svdvals!(copy(A), copy(B))
-svdvals(A::StridedMatrix, B::StridedMatrix) = svdvals!(float(A), float(B))
+svdvals{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = svdvals!(copy(A),copy(B))
+svdvals{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(sqrt(one(TA))),TB); svdvals!(S != TA ? convert(Matrix{S}, A) : copy(A), S != TB ? convert(Matrix{S}, B) : copy(B)))
 
-type Schur{Ty<:BlasFloat} <: Factorization{Ty}
+immutable Schur{Ty<:BlasFloat} <: Factorization{Ty}
     T::Matrix{Ty}
     Z::Matrix{Ty}
     values::Vector
 end
 
 schurfact!{T<:BlasFloat}(A::StridedMatrix{T}) = Schur(LinAlg.LAPACK.gees!('V', A)...)
-schurfact!(A::StridedMatrix) = schurfact!(float(A))
 schurfact{T<:BlasFloat}(A::StridedMatrix{T}) = schurfact!(copy(A))
-schurfact(A::StridedMatrix) = schurfact!(float(A))
+schurfact{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? schurfact!(convert(Matrix{S},A)) : schurfact!(copy(A)))
 
 function getindex(F::Schur, d::Symbol)
     (d == :T || d == :Schur) && return F.T
@@ -702,7 +887,7 @@ function schur(A::AbstractMatrix)
     SchurF[:T], SchurF[:Z], SchurF[:values]
 end
 
-type GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
+immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
     S::Matrix{Ty}
     T::Matrix{Ty}
     alpha::Vector
@@ -712,9 +897,8 @@ type GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
 end
 
 schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = GeneralizedSchur(LinAlg.LAPACK.gges!('V', 'V', A, B)...)
-schurfact!(A::StridedMatrix, B::StridedMatrix) = schurfact!(float(A), float(B))
-schurfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = schurfact!(copy(A), copy(B))
-schurfact(A::StridedMatrix, B::StridedMatrix) = schurfact!(float(A), float(B))
+schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = schurfact!(copy(A),copy(B))
+schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(sqrt(one(TA))),TB); schurfact!(S != TA ? convert(Matrix{S},A) : copy(A), S != TB ? convert(Matrix{S},B) : copy(B)))
 
 function getindex(F::GeneralizedSchur, d::Symbol)
     d == :S && return F.S
diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index 11c22a481ca0d..298e6ca8c849f 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -240,3 +240,67 @@ function axpy!{Ti<:Integer,Tj<:Integer}(alpha, x::AbstractArray, rx::AbstractArr
     y
 end
 
+# Elementary reflection similar to LAPACK. The reflector is not Hermitian but ensures that tridiagonalization of Hermitian matrices become real. See lawn72
+function elementaryLeft!(A::AbstractMatrix, row::Integer, col::Integer)
+    m, n = size(A)
+    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
+    1 <= col <= n || throw(BoundsError("col cannot be less than one or larger than $(size(A,2))"))
+    @inbounds begin
+        ξ1 = A[row,col]
+        normu = abs2(ξ1)
+        for i = row+1:m
+            normu += abs2(A[i,col])
+        end
+        normu = sqrt(normu)
+        ν = copysign(normu,real(ξ1))
+        A[row,col] += ν
+        ξ1 += ν
+        A[row,col] = -ν
+        for i = row+1:m
+            A[i,col] /= ξ1
+        end
+    end
+    ξ1/ν
+end
+function elementaryRight!(A::AbstractMatrix, row::Integer, col::Integer)
+    m, n = size(A)
+    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
+    1 <= col <= n || throw(BoundsError("col cannot be less than one or larger than $(size(A,2))"))
+    row <= col || error("col cannot be larger than row")
+    @inbounds begin
+        ξ1 = A[row,col]
+        normu = abs2(ξ1)
+        for i = col+1:n
+            normu += abs2(A[row,i])
+        end
+        normu = sqrt(normu)
+        ν = copysign(normu,real(ξ1))
+        A[row,col] += ν
+        ξ1 += ν
+        A[row,col] = -ν
+        for i = col+1:n
+            A[row,i] /= ξ1
+        end
+    end
+    conj(ξ1/ν)
+end
+function elementaryRightTrapezoid!(A::AbstractMatrix, row::Integer)
+    m, n = size(A)
+    1 <= row <= m || throw(BoundsError("row cannot be less than one or larger than $(size(A,1))"))
+    @inbounds begin
+        ξ1 = A[row,row]
+        normu = abs2(A[row,row])
+        for i = m+1:n
+            normu += abs2(A[row,i])
+        end
+        normu = sqrt(normu)
+        ν = copysign(normu,real(ξ1))
+        A[row,row] += ν
+        ξ1 += ν
+        A[row,row] = -ν
+        for i = m+1:n
+            A[row,i] /= ξ1
+        end
+    end
+    conj(ξ1/ν)
+end
\ No newline at end of file
diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl
index ad1f37af7ff20..dd7e03f9bbe20 100644
--- a/base/linalg/lapack.jl
+++ b/base/linalg/lapack.jl
@@ -286,7 +286,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty
             minmn = min(m, n)
             nb = size(T, 1)
             nb <= minmn || throw(ArgumentError("Block size $nb > $minmn too large"))
-            lda = max(1, m)
+            lda = max(1, stride(A,2))
             work = Array($elty, nb*n)
             if n > 0
                 info = Array(BlasInt, 1)
@@ -295,7 +295,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty
                      Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
                      Ptr{BlasInt}),
                      &m, &n, &nb, A, 
-                     &lda, T, &nb, work,
+                     &lda, T, &max(1,stride(T,2)), work,
                      info)
                 @lapackerror
             end
@@ -455,8 +455,8 @@ for (tzrzf, ormrz, elty) in
             m, n = size(C)
             k = length(tau)
             l = size(A, 2) - size(A, 1)
-            lda = max(1, k)
-            ldc = max(1, m)
+            lda = max(1, stride(A,2))
+            ldc = max(1, stride(C,2))
             work = Array($elty, 1)
             lwork = -1
             info = Array(BlasInt, 1)
@@ -1661,7 +1661,7 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in
                 wss = m*k
             end
             1 <= nb <= k || throw(DimensionMismatch("Wrong value for nb"))
-            ldc = max(1, m)
+            ldc = max(1, stride(C,2))
             work = Array($elty, wss)
             info = Array(BlasInt, 1)
             ccall(($(string(gemqrt)), liblapack), Void,
@@ -1671,7 +1671,7 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in
                  Ptr{$elty}, Ptr{BlasInt}),
                 &side, &trans, &m, &n,
                 &k, &nb, V, &ldv,
-                T, &nb, C, &ldc,
+                T, &max(1,stride(T,2)), C, &ldc,
                 work, info)
             @lapackerror
             return C
diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
index e6756290faa58..f34f9841bd01b 100644
--- a/base/linalg/symmetric.jl
+++ b/base/linalg/symmetric.jl
@@ -1,31 +1,24 @@
 #Symmetric and Hermitian matrices
-
-for ty in (:Hermitian, :Symmetric)
-    @eval begin
-        type $ty{T} <: AbstractMatrix{T}
-            S::Matrix{T}
-            uplo::Char
-        end
-        function $ty(S::Matrix, uplo::Symbol=:U)
-            chksquare(S)
-            $ty(S, string(uplo)[1])
-        end
-
-        function copy!(A::$ty, B::$ty)
-            copy!(A.S, B.S)
-            A.uplo = B.uplo
-            A
-        end
-        similar(A::$ty, args...) = $ty(similar(A.S, args...), A.uplo)
-    end
+immutable Symmetric{T} <: AbstractMatrix{T}
+    S::Matrix{T}
+    uplo::Char
 end
-
-typealias HermOrSym Union(Hermitian, Symmetric)
+Symmetric(A::Matrix, uplo::Symbol=:U) = (chksquare(A);Symmetric(A, string(uplo)[1]))
+immutable Hermitian{T} <: AbstractMatrix{T}
+    S::Matrix{T}
+    uplo::Char
+end
+Hermitian(A::Matrix, uplo::Symbol=:U) = (chksquare(A);Hermitian(A, string(uplo)[1]))
+typealias HermOrSym{T} Union(Hermitian{T}, Symmetric{T})
 
 size(A::HermOrSym, args...) = size(A.S, args...)
 getindex(A::HermOrSym, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : conj(getindex(A.S, j, i))
-full(A::Hermitian) = copytri!(A.S, A.uplo, true)
 full(A::Symmetric) = copytri!(A.S, A.uplo)
+full(A::Hermitian) = copytri!(A.S, A.uplo, true)
+convert{T}(::Type{Symmetric{T}},A::Symmetric) = Symmetric(convert(Matrix{T},A.S),A.uplo)
+convert{T}(::Type{Hermitian{T}},A::Hermitian) = Hermitian(convert(Matrix{T},A.S),A.uplo)
+copy(A::Symmetric) = Symmetric(copy(A.S),A.uplo)
+copy(A::Hermitian) = Hermitian(copy(A.S),A.uplo)
 ishermitian(A::Hermitian) = true
 ishermitian{T<:Real}(A::Symmetric{T}) = true
 ishermitian{T<:Complex}(A::Symmetric{T}) = all(imag(A.S) .== 0)
@@ -39,41 +32,37 @@ ctranspose(A::Hermitian) = A
 *(A::HermOrSym, B::StridedMatrix) = full(A)*B
 *(A::StridedMatrix, B::HermOrSym) = A*full(B)
 
-factorize!(A::HermOrSym) = bkfact!(A.S, symbol(A.uplo), issym(A))
+factorize(A::HermOrSym) = bkfact(A.S, symbol(A.uplo), issym(A))
 \(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.S, symbol(A.uplo), issym(A)), B)
 
-# eigvals!(A::HermOrSym, args...) = eigvals!(float(A), args...)
-eigvals!(A::HermOrSym) = eigvals!(A, 1, size(A, 1))
+eigfact!{T<:BlasReal}(A::Symmetric{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
+eigfact!{T<:BlasComplex}(A::Hermitian{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
+eigfact{T<:BlasFloat}(A::HermOrSym{T}) = eigfact!(copy(A))
+eigfact{T}(A::HermOrSym{T}) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? eigfact!(convert(typeof(A).name.primary{S}, A)) : eigfact!(copy(A)))
+eigvals!{T<:BlasReal}(A::Symmetric{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
+eigvals!{T<:BlasReal}(A::Symmetric{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
+eigvals!{T<:BlasComplex}(A::Hermitian{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
+eigvals!{T<:BlasComplex}(A::Hermitian{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
+eigvals{T<:BlasFloat}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = eigvals!(copy(A),l,h)
+eigvals{T}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = (S = promote_type(Float32,typeof(sqrt(one(T)))); S != T ? eigvals!(convert(typeof(A).name.primary{S}, A, l, h)) : eigvals!(copy(A), l, h))
 eigmax(A::HermOrSym) = eigvals(A, size(A, 1), size(A, 1))[1]
 eigmin(A::HermOrSym) = eigvals(A, 1, 1)[1]
-eigfact(A::HermOrSym) = eigfact!(copy(A))
 
-for (elty, ty) in ((:BlasFloat, :Hermitian), (:BlasReal, :Symmetric))
-    @eval begin
-        eigfact!{T<:$elty }(A::$ty{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
-        eigvals!{T<:$elty}(A::$ty{T}, il::Int, ih::Int) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
-        eigvals!{T<:$elty}(A::$ty{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
-        function eigfact!{T<:$elty}(A::$ty{T}, B::$ty{T})
-            vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
-            GeneralizedEigen(vals, vecs)
-        end
-        eigfact(A::$ty, B::$ty) = eigfact!(copy(A), copy(B))
-        eigvals!{T<:$elty}(A::$ty{T}, B::$ty{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
-    end
+function eigfact!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T})
+    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
+    GeneralizedEigen(vals, vecs)
+end
+function eigfact!{T<:BlasComplex}(A::Hermitian{T}, B::Hermitian{T})
+    vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')
+    GeneralizedEigen(vals, vecs)
 end
+eigvals!{T<:BlasReal}(A::Symmetric{T}, B::Symmetric{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
+eigvals!{T<:BlasComplex}(A::Hermitian{T}, B::Hermitian{T}) = LAPACK.sygvd!(1, 'N', A.uplo, A.S, B.uplo == A.uplo ? B.S : B.S')[1]
 
 #Matrix-valued functions
-for (elty, ty) in ((:Any, :Hermitian), (:Real, :Symmetric))
-    @eval begin
-        function expm{T<:$elty}(A::$ty{T})
-            F = eigfact(A)
-            F.vectors * Diagonal(exp(F.values)) * F.vectors'
-        end
-
-        function sqrtm{T<:$elty}(A::$ty{T}, cond::Bool=false)
-            F = eigfact(A)
-            retmat = F.vectors*Diagonal((isposdef(F) ? sqrt : x->sqrt(complex(x)))(F.values))*F.vectors'
-            return cond ? (retmat, norm(vsqrt, Inf)^2/norm(F[:values], Inf)) : retmat
-        end
-    end
-end
\ No newline at end of file
+expm(A::HermOrSym) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
+function sqrtm(A::HermOrSym)
+    F = eigfact(A)
+    isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'
+    return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors'
+end
diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl
index f29052f3b5d6e..5659dab97454e 100644
--- a/base/linalg/tridiag.jl
+++ b/base/linalg/tridiag.jl
@@ -13,7 +13,7 @@ end
 SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(copy(dv), copy(ev))
 
 function SymTridiagonal{Td,Te}(dv::Vector{Td}, ev::Vector{Te})
-    T = promote(Td,Te)
+    T = promote_type(Td,Te)
     SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
 end
 
@@ -231,9 +231,11 @@ convert(::Type{Tridiagonal}, A::SymTridiagonal) = Tridiagonal(A.ev, A.dv, A.ev)
 -(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl-B.ev, A.d-B.dv, A.du-B.ev)
 -(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev-B.dl, A.dv-B.d, A.ev-B.du)
 
+convert{T}(::Type{Tridiagonal{T}},M::Tridiagonal) = Tridiagonal(convert(Vector{T}, M.dl), convert(Vector{T}, M.d), convert(Vector{T}, M.du))
 convert{T}(::Type{Tridiagonal{T}}, M::SymTridiagonal{T}) = Tridiagonal(M)
 convert{T}(::Type{SymTridiagonal{T}}, M::Tridiagonal) = M.dl==M.du ? (SymTridiagonal(M.dl, M.d)) :
     error("Tridiagonal is not symmetric, cannot convert to SymTridiagonal")
+convert{T}(::Type{SymTridiagonal{T}},M::SymTridiagonal) = SymTridiagonal(convert(Vector{T}, M.dv), convert(Vector{T}, M.ev))
 
 ## Solvers
 
@@ -367,9 +369,9 @@ end
 LDLTTridiagonal{S<:BlasFloat,T<:BlasFloat}(D::Vector{S}, E::Vector{T}) = LDLTTridiagonal{T,S}(D, E)
 
 ldltd!{T<:BlasFloat}(A::SymTridiagonal{T}) = LDLTTridiagonal(LAPACK.pttrf!(real(A.dv),A.ev)...)
-ldltd!{T<:Integer}(A::SymTridiagonal{T}) = ldltd!(SymTridiagonal(float(A.dv),float(A.ev)))
-ldltd(A::SymTridiagonal) = ldltd!(copy(A))
-factorize!(A::SymTridiagonal) = ldltd(A)
+ldltd{T<:BlasFloat}(A::SymTridiagonal{T}) = ldltd!(copy(A))
+ldltd{T}(A::SymTridiagonal{T}) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? ldltd!(convert(SymTridiagonal{S},A)) : ldltd!(copy(A)))
+factorize(A::SymTridiagonal) = ldltd(A)
 
 A_ldiv_B!{T<:BlasReal}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) = LAPACK.pttrs!(C.D, C.E, B)
 A_ldiv_B!{T<:BlasComplex}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) = LAPACK.pttrs!('L', C.D, C.E, B)
@@ -391,9 +393,9 @@ type LUTridiagonal{T} <: Factorization{T}
     # end
 end
 lufact!{T<:BlasFloat}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
-lufact!{T<:Union(Rational,Integer)}(A::Tridiagonal{T}) = lufact!(float(A))
-lufact(A::Tridiagonal) = lufact!(copy(A))
-factorize!(A::Tridiagonal) = lufact!(A)
+lufact{T<:BlasFloat}(A::Tridiagonal{T}) = lufact!(copy(A))
+lufact{T}(A::Tridiagonal{T}) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? lufact!(convert(Tridiagonal{S},A)) : lufact!(copy(A)))
+factorize(A::Tridiagonal) = lufact(A)
 #show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))
 
 function det{T}(lu::LUTridiagonal{T})
diff --git a/test/linalg.jl b/test/linalg.jl
index 45c4d7ddd3e01..64baeacc788fc 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -1,3 +1,5 @@
+debug = false
+
 import Base.LinAlg
 import Base.LinAlg: BlasComplex, BlasFloat, BlasReal
 
@@ -14,180 +16,202 @@ for elty in (Float32, Float64, Complex64, Complex128)
     @test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01
 end
 
-areal = randn(n,n)
-aimg  = randn(n,n) 
-breal = randn(n,2)
-bimg  = randn(n,2)
-for elty in (Float32, Float64, Complex64, Complex128, Int)
-    if elty == Complex64 || elty == Complex128
-        a = complex(areal, aimg)
-        b = complex(breal, bimg)
-    else
-        a = areal
-        b = breal
+areal = randn(n,n)/2
+aimg  = randn(n,n)/2
+breal = randn(n,2)/2
+bimg  = randn(n,2)/2
+for eltya in (Float16, Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
+    for eltyb in (Float16, Float32, Float64, Complex32, Complex64, Complex128, Int)
+        a = eltya == Int ? rand(1:5, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
+        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
+        asym = a'+a                  # symmetric indefinite
+        apd  = a'*a                 # symmetric positive-definite
+
+        ε = max(eps(abs(float(one(eltya)))),eps(abs(float(one(eltyb)))))
+
+debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
+
+debug && println("(Automatic) upper Cholesky factor")
+    if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement cholesky decomposition in julia
+        capd  = factorize(apd)
+        r     = capd[:U]
+        @test_approx_eq r'*r apd
+        @test_approx_eq_eps b apd * (capd\b) 6500ε
+        @test_approx_eq apd * inv(capd) eye(n)
+        @test_approx_eq_eps a*(capd\(a'*b)) b 800ε
+        @test_approx_eq det(capd) det(apd)
+        @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow
+
+debug && println("lower Cholesky factor")
+        l = cholfact(apd, :L)[:L] 
+        @test_approx_eq l*l' apd
+    end
+
+debug && println("pivoted Choleksy decomposition")
+    if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted cholesky decomposition in julia
+        cpapd = cholfact(apd, pivot=true)
+        @test rank(cpapd) == n
+        @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
+        @test_approx_eq_eps b apd * (cpapd\b) 15000ε
+        if isreal(apd)
+            @test_approx_eq apd * inv(cpapd) eye(n)
+        end
     end
-    if elty <: BlasFloat
-        a = convert(Matrix{elty}, a)
-        b = convert(Matrix{elty}, b)
-    else
-        a = rand(1:10, n, n)
-        b = rand(1:10, n, 2)
-    end
-    asym = a' + a                  # symmetric indefinite
-    apd  = a'*a                    # symmetric positive-definite
-
-    # (Automatic) upper Cholesky factor
-    capd  = factorize(apd)
-    r     = capd[:U]
-    @test_approx_eq r'*r apd
-    @test_approx_eq b apd * (capd\b)
-    @test_approx_eq apd * inv(capd) eye(elty, n)
-    @test_approx_eq a*(capd\(a'*b)) b           # least squares soln for square a
-    @test_approx_eq det(capd) det(apd)
-    @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow
-
-    # lower Cholesky factor
-    l = cholfact(apd, :L)[:L] 
-    @test_approx_eq l*l' apd
-
-    # pivoted Choleksy decomposition
-    cpapd = cholfact(apd, pivot=true)
-    @test rank(cpapd) == n
-    @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
-    @test_approx_eq b apd * (cpapd\b)
-    if isreal(apd)
-        @test_approx_eq apd * inv(cpapd) eye(elty, n)
-    end
-
-    # (Automatic) Bunch-Kaufman factor of indefinite matrix
-    bc1 = factorize(asym) 
-    @test_approx_eq inv(bc1) * asym eye(elty, n)
-    @test_approx_eq asym * (bc1\b) b
-
-    # Bunch-Kaufman factors of a pos-def matrix
-    bc2 = bkfact(apd)    
-    @test_approx_eq inv(bc2) * apd eye(elty, n)
-    @test_approx_eq apd * (bc2\b) b
-
-    # (Automatic) Square LU decomposition
+
+debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix")
+    if eltya != BigFloat && eltyb != BigFloat # Not implemented for BigFloat and I don't think it will.
+        bc1 = factorize(asym)
+        @test_approx_eq inv(bc1) * asym eye(n)
+        @test_approx_eq_eps asym * (bc1\b) b 600ε
+
+debug && println("Bunch-Kaufman factors of a pos-def matrix")
+        bc2 = bkfact(apd)
+        @test_approx_eq inv(bc2) * apd eye(n)
+        @test_approx_eq_eps apd * (bc2\b) b 6000ε
+    end
+
+debug && println("(Automatic) Square LU decomposition")
     lua   = factorize(a)
     l,u,p = lua[:L], lua[:U], lua[:p]
     @test_approx_eq l*u a[p,:]
     @test_approx_eq l[invperm(p),:]*u a
-    @test_approx_eq a * inv(lua) eye(elty, n)
-    @test_approx_eq a*(lua\b) b
+    @test_approx_eq a * inv(lua) eye(n)
+    @test_approx_eq_eps a*(lua\b) b 80ε
 
-    # Thin LU
+debug && println("Thin LU")
     lua   = lufact(a[:,1:5])
     @test_approx_eq lua[:L]*lua[:U] lua[:P]*a[:,1:5]
 
-    # Fat LU
+debug && println("Fat LU")
     lua   = lufact(a[1:5,:])
     @test_approx_eq lua[:L]*lua[:U] lua[:P]*a[1:5,:]
 
-    # QR decomposition
-    qra   = qrfact(a)
+debug && println("QR decomposition (without pivoting)")
+    qra   = qrfact(a, pivot=false)
     q,r   = qra[:Q], qra[:R]
-    @test_approx_eq q'*full(q, thin=false) eye(elty, n)
-    @test_approx_eq q*full(q, thin=false)' eye(elty, n)
+    @test_approx_eq q'*full(q, thin=false) eye(n)
+    @test_approx_eq q*full(q, thin=false)' eye(n)
     @test_approx_eq q*r a
-    @test_approx_eq a*(qra\b) b
+    @test_approx_eq_eps a*(qra\b) b 3000ε
 
-    # (Automatic) Fat pivoted QR decomposition
+debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
     qrpa  = factorize(a[1:5,:])
-    q,r,p = qrpa[:Q], qrpa[:R], qrpa[:p]
-    @test_approx_eq q'*full(q, thin=false) eye(elty, 5)
-    @test_approx_eq q*full(q, thin=false)' eye(elty, 5)
-    @test_approx_eq q*r a[1:5,p]
-    @test_approx_eq q*r[:,invperm(p)] a[1:5,:]
-    @test_approx_eq a[1:5,:]*(qrpa\b[1:5]) b[1:5]
-
-    # (Automatic) Thin pivoted QR decomposition
+    q,r = qrpa[:Q], qrpa[:R]
+    if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
+    @test_approx_eq q'*full(q, thin=false) eye(5)
+    @test_approx_eq q*full(q, thin=false)' eye(5)
+    @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:5,p] : a[1:5,:]
+    @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:5,:]
+    @test_approx_eq_eps a[1:5,:]*(qrpa\b[1:5]) b[1:5] 5000ε
+
+debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
     qrpa  = factorize(a[:,1:5])
-    q,r,p = qrpa[:Q], qrpa[:R], qrpa[:p]
-    @test_approx_eq q'*full(q, thin=false) eye(elty, n)
-    @test_approx_eq q*full(q, thin=false)' eye(elty, n)
-    @test_approx_eq q*r a[:,p]
-    @test_approx_eq q*r[:,invperm(p)] a[:,1:5]
-
-    # symmetric eigen-decomposition
-    d,v   = eig(asym)
-    @test_approx_eq asym*v[:,1] d[1]*v[:,1]
-    @test_approx_eq v*scale(d,v') asym
-
-    # non-symmetric eigen decomposition
-    d,v   = eig(a)
-    for i in 1:size(a,2) @test_approx_eq a*v[:,i] d[i]*v[:,i] end
-
-    # symmetric generalized eigenproblem
-    a610 = a[:,6:10]
-    f = eigfact(asym[1:5,1:5], a610'a610)
-    @test_approx_eq asym[1:5,1:5]*f[:vectors] scale(a610'a610*f[:vectors], f[:values])
-    @test_approx_eq f[:values] eigvals(asym[1:5,1:5], a610'a610)
-    @test_approx_eq prod(f[:values]) prod(eigvals(asym[1:5,1:5]/(a610'a610)))
+    q,r = qrpa[:Q], qrpa[:R]
+    if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
+    @test_approx_eq q'*full(q, thin=false) eye(n)
+    @test_approx_eq q*full(q, thin=false)' eye(n)
+    @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:5]
+    @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:5]
+
+debug && println("symmetric eigen-decomposition")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        d,v   = eig(asym)
+        @test_approx_eq asym*v[:,1] d[1]*v[:,1]
+        @test_approx_eq v*scale(d,v') asym
+    end
+
+debug && println("non-symmetric eigen decomposition")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        d,v   = eig(a)
+        for i in 1:size(a,2) @test_approx_eq a*v[:,i] d[i]*v[:,i] end
+    end
+
+debug && println("symmetric generalized eigenproblem")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        a610 = a[:,6:10]
+        f = eigfact(asym[1:5,1:5], a610'a610)
+        @test_approx_eq asym[1:5,1:5]*f[:vectors] scale(a610'a610*f[:vectors], f[:values])
+        @test_approx_eq f[:values] eigvals(asym[1:5,1:5], a610'a610)
+        @test_approx_eq_eps prod(f[:values]) prod(eigvals(asym[1:5,1:5]/(a610'a610))) 15ε
+    end
  
-    # Non-symmetric generalized eigenproblem
-    f = eigfact(a[1:5,1:5], a[6:10,6:10])
-    @test_approx_eq a[1:5,1:5]*f[:vectors] scale(a[6:10,6:10]*f[:vectors], f[:values])
-    @test_approx_eq f[:values] eigvals(a[1:5,1:5], a[6:10,6:10])
-    @test_approx_eq prod(f[:values]) prod(eigvals(a[1:5,1:5]/a[6:10,6:10]))
-
-    # Schur
-    f = schurfact(a)
-    @test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a
-    @test_approx_eq sort(real(f[:values])) sort(real(d))
-    @test_approx_eq sort(imag(f[:values])) sort(imag(d))
-    @test istriu(f[:Schur]) || iseltype(a,Real)
-
-    # Generalized Schur
-    f = schurfact(a[1:5,1:5], a[6:10,6:10]) 
-    @test_approx_eq f[:Q]*f[:S]*f[:Z]' a[1:5,1:5]
-    @test_approx_eq f[:Q]*f[:T]*f[:Z]' a[6:10,6:10]
-    @test istriu(f[:S]) || iseltype(a,Real)
-    @test istriu(f[:T]) || iseltype(a,Real)
-
-    # singular value decomposition
-    usv = svdfact(a)                
-    @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a
-  
-    # Generalized svd
-    gsvd = svdfact(a,a[1:5,:])
-    @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a
-    @test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a[1:5,:]
+debug && println("Non-symmetric generalized eigenproblem")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        f = eigfact(a[1:5,1:5], a[6:10,6:10])
+        @test_approx_eq a[1:5,1:5]*f[:vectors] scale(a[6:10,6:10]*f[:vectors], f[:values])
+        @test_approx_eq f[:values] eigvals(a[1:5,1:5], a[6:10,6:10])
+        @test_approx_eq_eps prod(f[:values]) prod(eigvals(a[1:5,1:5]/a[6:10,6:10])) 50000ε
+    end
 
+debug && println("Schur")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        f = schurfact(a)
+        @test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a
+        @test_approx_eq sort(real(f[:values])) sort(real(d))
+        @test_approx_eq sort(imag(f[:values])) sort(imag(d))
+        @test istriu(f[:Schur]) || iseltype(a,Real)
+    end
+
+debug && println("Generalized Schur")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        f = schurfact(a[1:5,1:5], a[6:10,6:10]) 
+        @test_approx_eq f[:Q]*f[:S]*f[:Z]' a[1:5,1:5]
+        @test_approx_eq f[:Q]*f[:T]*f[:Z]' a[6:10,6:10]
+        @test istriu(f[:S]) || iseltype(a,Real)
+        @test istriu(f[:T]) || iseltype(a,Real)
+    end
+
+debug && println("singular value decomposition")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        usv = svdfact(a)
+        @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a
+    end
+
+debug && println("Generalized svd")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        gsvd = svdfact(a,a[1:5,:])
+        @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a
+        @test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a[1:5,:]
+    end
+
+debug && println("Solve square general system of equations")
     x = a \ b
-    @test_approx_eq a*x b
-    
+    @test_approx_eq_eps a*x b 80ε
+
+debug && println("Solve upper trianguler system")
     x = triu(a) \ b
-    @test_approx_eq triu(a)*x b
-    
+    @test_approx_eq_eps triu(a)*x b 15000ε
+
+debug && println("Solve lower triangular system")
     x = tril(a)\b
-    @test_approx_eq tril(a)*x b
-   
-    # Test null
-    a15null = null(a[:,1:5]')
-    @test rank([a[:,1:5] a15null]) == 10
-    @test_approx_eq_eps norm(a[:,1:5]'a15null) zero(elty) elty <: Union(Float32, Complex64) ? 1f-5 : 1e-13
-    @test_approx_eq_eps norm(a15null'a[:,1:5]) zero(elty) elty <: Union(Float32, Complex64) ? 1f-5 : 1e-13
-    @test size(null(b), 2) == 0
-
-    # Test pinv
-    pinva15 = pinv(a[:,1:5])
-    @test_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
-    @test_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
-    
-    # Complex vector rhs
-    x = a\complex(b)
-    @test_approx_eq a*x complex(b)
+    @test_approx_eq_eps tril(a)*x b 1000ε
+
+debug && println("Test null")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        a15null = null(a[:,1:5]')
+        @test rank([a[:,1:5] a15null]) == 10
+        @test_approx_eq_eps norm(a[:,1:5]'a15null, Inf) zero(eltya) 30ε
+        @test_approx_eq_eps norm(a15null'a[:,1:5], Inf) zero(eltya) 40ε
+        @test size(null(b), 2) == 0
+    end
+
+debug && println("Test pinv")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
+        pinva15 = pinv(a[:,1:5])
+        @test_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
+        @test_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
+    end
 
-    if isreal(a)
-        # Matrix square root
+    # if isreal(a)
+debug && println("Matrix square root")
+    if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
         asq = sqrtm(a)
         @test_approx_eq asq*asq a
         asymsq = sqrtm(asym)
         @test_approx_eq asymsq*asymsq asym
     end
 end
+end
 
 ## Least squares solutions
 a = [ones(20) 1:20 1:20]
@@ -827,6 +851,30 @@ for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128)
     @test_approx_eq gradient(x) g
 end
 
+# Test our own linear algebra functionaly against LAPACK
+for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
+    for nn in (5,10,15)
+        if elty <: Real
+            A = convert(Matrix{elty}, randn(10,nn))
+        else
+            A = convert(Matrix{elty}, complex(randn(10,nn),randn(10,nn)))
+        end    ## LU (only equal for real because LAPACK uses difference absolute value when choosing permutations)
+        if elty <: Real
+            FJulia  = invoke(lufact!, (AbstractMatrix,), copy(A)) 
+            FLAPACK = Base.LinAlg.LAPACK.getrf!(copy(A))
+            @test_approx_eq FJulia.factors FLAPACK[1]
+            @test_approx_eq FJulia.ipiv FLAPACK[2]
+            @test_approx_eq FJulia.info FLAPACK[3]
+        end
+        
+        ## QR
+        FJulia  = invoke(qrfact!, (AbstractMatrix,), copy(A)) 
+        FLAPACK = Base.LinAlg.LAPACK.geqrf!(copy(A))
+        @test_approx_eq FJulia.factors FLAPACK[1]
+        @test_approx_eq FJulia.τ FLAPACK[2]
+    end
+end
+
 # Test rational matrices
 ## Integrate in general tests when more linear algebra is implemented in julia
 a = convert(Matrix{Rational{BigInt}}, rand(1:10//1,n,n))/n