diff --git a/base/exports.jl b/base/exports.jl
index aaf51cf2881c3..c862ebfda3323 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -620,6 +620,8 @@ export
     istril,
     istriu,
     kron,
+    ldltfact,
+    ldltfact!,
     linreg,
     logdet,
     lu,
diff --git a/base/linalg.jl b/base/linalg.jl
index 1e13a4c2f0eed..7abb721544057 100644
--- a/base/linalg.jl
+++ b/base/linalg.jl
@@ -24,7 +24,7 @@ export
     Hessenberg,
     LU,
     LUTridiagonal,
-    LDLTTridiagonal,
+    LDLt,
     QR,
     QRPivoted,
     Schur,
@@ -78,8 +78,8 @@ export
     istril,
     istriu,
     kron,
-    ldltd!,
-    ldltd,
+    ldltfact!,
+    ldltfact,
     linreg,
     logdet,
     lu,
@@ -196,13 +196,14 @@ include("linalg/matmul.jl")
 include("linalg/lapack.jl")
 
 include("linalg/dense.jl")
+include("linalg/tridiag.jl")
 include("linalg/factorization.jl")
+include("linalg/lu.jl")
 
 include("linalg/bunchkaufman.jl")
 include("linalg/triangular.jl")
 include("linalg/symmetric.jl")
 include("linalg/woodbury.jl")
-include("linalg/tridiag.jl")
 include("linalg/diagonal.jl")
 include("linalg/bidiag.jl")
 include("linalg/uniformscaling.jl")
@@ -210,6 +211,7 @@ include("linalg/rectfullpacked.jl")
 include("linalg/givens.jl")
 include("linalg/special.jl")
 include("linalg/bitarray.jl")
+include("linalg/ldlt.jl")
 
 include("linalg/sparse.jl")
 include("linalg/umfpack.jl")
diff --git a/base/linalg/bunchkaufman.jl b/base/linalg/bunchkaufman.jl
index 3a4fde9b95a95..36c873ba6a010 100644
--- a/base/linalg/bunchkaufman.jl
+++ b/base/linalg/bunchkaufman.jl
@@ -22,6 +22,8 @@ bkfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym
 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)
+convert{T}(::Type{Factorization{T}}, B::BunchKaufman) = convert(BunchKaufman{T}, B)
+
 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/cholmod.jl b/base/linalg/cholmod.jl
index 28768126f963a..78bb732d48b12 100644
--- a/base/linalg/cholmod.jl
+++ b/base/linalg/cholmod.jl
@@ -18,7 +18,7 @@ import Base: (*), convert, copy, ctranspose, eltype, findnz, getindex, hcat,
 
 import ..LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
                  cholfact, cholfact!, copy, det, diag,
-                 full, logdet, norm, scale, scale!, solve, sparse
+                 full, logdet, norm, scale, scale!, sparse
 
 include("cholmod_h.jl")
 
diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl
index 7d7cc275a4675..b788c909285a2 100644
--- a/base/linalg/diagonal.jl
+++ b/base/linalg/diagonal.jl
@@ -1,11 +1,20 @@
 ## Diagonal matrices
 
-type Diagonal{T} <: AbstractMatrix{T}
+immutable Diagonal{T} <: AbstractMatrix{T}
     diag::Vector{T}
 end
 Diagonal(A::Matrix) = Diagonal(diag(A))
 
-convert{T<:Number,S<:Number}(::Type{Diagonal{T}}, A::Diagonal{S}) = T == S ? A : Diagonal(convert(Vector{T}, A.diag))
+convert{T}(::Type{Diagonal{T}}, D::Diagonal{T}) = D
+convert{T}(::Type{Diagonal{T}}, D::Diagonal) = Diagonal{T}(convert(Vector{T}, D.diag))
+convert{T}(::Type{AbstractMatrix{T}}, D::Diagonal) = convert(Diagonal{T}, D)
+
+function similar{T}(D::Diagonal, ::Type{T}, d::(Int,Int))
+    d[1] == d[2] || throw(ArgumentError("Diagonal matrix must be square"))
+    return Diagonal{T}(Array(T,d[1]))
+end
+
+copy!(D1::Diagonal, D2::Diagonal) = (copy!(D1.diag, D2.diag); D1)
 
 size(D::Diagonal) = (length(D.diag),length(D.diag))
 size(D::Diagonal,d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(D.diag) : 1)
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 8a7e656fe53df..4037fc2b54664 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -43,8 +43,10 @@ cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x)
 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)
+convert{T}(::Type{Cholesky{T}},C::Cholesky) = Cholesky(convert(Matrix{T},C.UL),C.uplo)
+convert{T}(::Type{Factorization{T}}, C::Cholesky) = convert(Cholesky{T}, C)
+convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = CholeskyPivoted(convert(Matrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info)
+convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C)
 
 size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
 size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)
@@ -120,139 +122,6 @@ chkfullrank(C::CholeskyPivoted) = C.rank<size(C.UL, 1) && throw(RankDeficientExc
 
 rank(C::CholeskyPivoted) = C.rank
 
-####################
-# LU Factorization #
-####################
-immutable LU{T} <: Factorization{T}
-    factors::Matrix{T}
-    ipiv::Vector{BlasInt}
-    info::BlasInt
-end
-
-lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(LAPACK.getrf!(A)...)
-function lufact!{T}(A::AbstractMatrix{T})
-    m, n = size(A)
-    minmn = min(m,n)
-    info = 0
-    ipiv = Array(BlasInt, minmn)
-    for k = 1:minmn
-        # find index max
-        kp = k
-        amax = real(zero(T))
-        for i = k:m
-            absi = abs(A[i,k])
-            if absi > amax
-                kp = i
-                amax = absi
-            end
-        end
-        ipiv[k] = kp
-        if A[kp,k] != 0
-            if k != kp
-                # Interchange
-                for i = 1:n
-                    tmp = A[k,i]
-                    A[k,i] = A[kp,i]
-                    A[kp,i] = tmp
-                end
-            end
-            # Scale first column
-            Akkinv = inv(A[k,k])
-            for i = k+1:m
-                A[i,k] *= Akkinv
-            end
-        elseif info == 0
-            info = k
-        end
-        # Update the rest
-        for j = k+1:n
-            for i = k+1:m
-                A[i,j] -= A[i,k]*A[k,j]
-            end
-        end        
-    end
-    LU(A, ipiv, convert(BlasInt, info))
-end
-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))
-    F = lufact(A)
-    F[:L], F[:U], F[:p]
-end
-
-convert{T}(::Type{LU{T}}, F::LU) = LU(convert(Matrix{T}, F.factors), F.ipiv, F.info)
-
-size(A::LU) = size(A.factors)
-size(A::LU,n) = size(A.factors,n)
-
-function ipiv2perm{T}(v::AbstractVector{T}, maxi::Integer)
-    p = T[1:maxi]
-    @inbounds for i in 1:length(v)
-        p[i], p[v[i]] = p[v[i]], p[i]
-    end
-    return p
-end
-
-function getindex{T}(A::LU{T}, d::Symbol)
-    m, n = size(A)
-    d == :L && return tril(A.factors[1:m, 1:min(m,n)], -1) + eye(T, m, min(m,n))
-    d == :U && return triu(A.factors[1:min(m,n),1:n])
-    d == :p && return ipiv2perm(A.ipiv, m)
-    if d == :P
-        p = A[:p]
-        P = zeros(T, m, m)
-        for i in 1:m
-            P[i,p[i]] = one(T)
-        end
-        return P
-    end
-    throw(KeyError(d))
-end
-
-function det{T}(A::LU{T})
-    n = chksquare(A)
-    A.info > 0 && return zero(typeof(A.factors[1]))
-    return prod(diag(A.factors)) * (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
-end
-
-function logdet2{T<:Real}(A::LU{T})  # return log(abs(det)) and sign(det)
-    n = chksquare(A)
-    dg = diag(A.factors)
-    s = (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
-    sum(log(abs(dg))), s 
-end
-
-function logdet{T<:Real}(A::LU{T})
-    d,s = logdet2(A)
-    s>=0 || error("DomainError: determinant is negative")
-    d
-end
-
-function logdet{T<:Complex}(A::LU{T})
-    n = chksquare(A)
-    s = sum(log(diag(A.factors))) + (bool(sum(A.ipiv .!= 1:n) % 2) ? complex(0,pi) : 0) 
-    r, a = reim(s)
-    a = pi-mod(pi-a,2pi) #Take principal branch with argument (-pi,pi] 
-    complex(r,a)    
-end
-
-A_ldiv_B!{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info
-A_ldiv_B!(A::LU, b::StridedVector) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), b[ipiv2perm(A.ipiv, length(b))]))
-A_ldiv_B!(A::LU, B::StridedMatrix) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), B[ipiv2perm(A.ipiv, size(B, 1)),:]))
-At_ldiv_B{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, copy(B)) A.info
-Ac_ldiv_B{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, copy(B)) A.info
-At_ldiv_Bt{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info
-Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info
-
-/{T}(B::Matrix{T},A::LU{T}) = At_ldiv_Bt(A,B).'
-
-inv{T<:BlasFloat}(A::LU{T}) = @assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
-
-cond{T<:BlasFloat}(A::LU{T}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
-cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
-
 ####################
 # QR Factorization #
 ####################
@@ -307,8 +176,11 @@ function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true)
 end
 
 convert{T}(::Type{QR{T}},A::QR) = QR(convert(Matrix{T}, A.factors), convert(Vector{T}, A.τ))
+convert{T}(::Type{Factorization{T}}, A::QR) = convert(QR{T}, A)
 convert{T}(::Type{QRCompactWY{T}},A::QRCompactWY) = QRCompactWY(convert(Matrix{T}, A.factors), convert(Matrix{T}, A.T))
+convert{T}(::Type{Factorization{T}}, A::QRCompactWY) = convert(QRCompactWY{T}, A)
 convert{T}(::Type{QRPivoted{T}},A::QRPivoted) = QRPivoted(convert(Matrix{T}, A.factors), convert(Vector{T}, A.τ), A.jpvt)
+convert{T}(::Type{Factorization{T}}, A::QRPivoted) = convert(QRPivoted{T}, A)
 
 function getindex(A::QR, d::Symbol)
     d == :R && return triu(A.factors[1:minimum(size(A)),:])
@@ -345,8 +217,10 @@ immutable QRCompactWYQ{S} <: AbstractMatrix{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))
+convert{T}(::Type{QRPackedQ{T}}, Q::QRPackedQ) = QRPackedQ(convert(Matrix{T}, Q.factors), convert(Vector{T}, Q.τ))
+convert{T}(::Type{AbstractMatrix{T}}, Q::QRPackedQ) = convert(QRPackedQ{T}, Q)
+convert{S}(::Type{QRCompactWYQ{S}}, Q::QRCompactWYQ) = QRCompactWYQ(convert(Matrix{S}, Q.factors), convert(Matrix{S}, Q.T))
+convert{S}(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ) = convert(QRCompactWYQ{S}, Q)
 
 size(A::Union(QR,QRCompactWY,QRPivoted), dim::Integer) = size(A.factors, dim)
 size(A::Union(QR,QRCompactWY,QRPivoted)) = size(A.factors)
@@ -383,13 +257,13 @@ function A_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
 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)
+    Anew = convert(AbstractMatrix{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)
+    Anew = convert(AbstractMatrix{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
@@ -420,9 +294,9 @@ function Ac_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
     end
     B
 end
-function Ac_mul_B{TQ<:Number,TB<:Number}(Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}), B::StridedVecOrMat{TB})
+function Ac_mul_B{TQ<:Number,TB<:Number,N}(Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}), B::StridedArray{TB,N})
     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))
+    Ac_mul_B!(convert(AbstractMatrix{TQB}, Q), TB == TQB ? copy(B) : convert(AbstractArray{TQB,N}, B))
 end
 ### AQ
 A_mul_B!{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','N', B.factors, B.T, A)
@@ -449,9 +323,9 @@ function A_mul_B!{T}(A::StridedMatrix{T},Q::QRPackedQ{T})
     end
     A
 end
-function *{TA,TQ}(A::StridedVecOrMat{TA}, Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}))
+function *{TA,TQ,N}(A::StridedArray{TA,N}, Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}))
     TAQ = promote_type(TA, TQ)
-    A_mul_B!(TA==TAQ ? copy(A) : convert(Matrix{TAQ}, A), convert(typeof(Q).name.primary{TAQ}, Q))
+    A_mul_B!(TA==TAQ ? copy(A) : convert(AbstractArray{TAQ,N}, A), convert(AbstractMatrix{TAQ}, Q))
 end
 ### AQc
 A_mul_Bc!{T<:BlasReal}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','T',B.factors,B.T,A)
@@ -480,9 +354,9 @@ function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T})
     end
     A
 end
-function A_mul_Bc{TA,TB}(A::AbstractVecOrMat{TA}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB}))
+function A_mul_Bc{TA,TB,N}(A::AbstractArray{TA,N}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB}))
     TAB = promote_type(TA,TB)
-    A_mul_Bc!(size(A,2)==size(B.factors,1) ? (TA == TAB ? copy(A) : convert(Matrix{TAB}, A)) : (size(A,2)==size(B.factors,2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : throw(DimensionMismatch(""))),B)
+    A_mul_Bc!(size(A,2)==size(B.factors,1) ? (TA == TAB ? copy(A) : convert(AbstractMatrix{TAB,N}, A)) : (size(A,2)==size(B.factors,2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : throw(DimensionMismatch(""))),convert(AbstractMatrix{TAB}, B))
 end
 
 # Julia implementation similarly to xgelsy
@@ -572,13 +446,13 @@ function \{TA,Tb}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),b::StridedVecto
     S = promote_type(TA,Tb)
     m,n = size(A)
     m == length(b) || throw(DimensionMismatch("left hand side has $(m) rows, but right hand side has length $(length(b))"))
-    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))
+    n > m ? A_ldiv_B!(convert(Factorization{S},A),[b,zeros(S,n-m)]) : A_ldiv_B!(convert(Factorization{S},A), S == Tb ? copy(b) : convert(AbstractVector{S}, b))
 end
 function \{TA,TB}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),B::StridedMatrix{TB})
     S = promote_type(TA,TB)
     m,n = size(A)
     m == size(B,1) || throw(DimensionMismatch("left hand side has $(m) rows, but right hand side has $(size(B,1)) rows"))
-    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))
+    n > m ? A_ldiv_B!(convert(Factorization{S},A),[B;zeros(S,n-m,size(B,2))]) : A_ldiv_B!(convert(Factorization{S},A), S == TB ? copy(B) : convert(AbstractMatrix{S}, B))
 end
 
 ##TODO:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
@@ -917,18 +791,19 @@ function schur(A::AbstractMatrix, B::AbstractMatrix)
 end
 
 ### General promotion rules
+convert{T}(::Type{Factorization{T}}, F::Factorization{T}) = F
 inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1)))
-function \{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB})
+function \{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
     TFB = typeof(one(TF)/one(TB)) 
-    A_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B))
+    A_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
 end
 
-function Ac_ldiv_B{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB})
+function Ac_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
     TFB = typeof(one(TF)/one(TB)) 
-    Ac_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B))
+    Ac_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
 end
 
-function At_ldiv_B{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB})
+function At_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
     TFB = typeof(one(TF)/one(TB)) 
-    At_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B))
+    At_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
 end
diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index 1b2fbbe606f44..8f8268bd2c9ac 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -213,9 +213,9 @@ trace(x::Number) = x
 inv(a::AbstractVector) = error("argument must be a square matrix")
 inv{T}(A::AbstractMatrix{T}) = A_ldiv_B!(A,eye(T, chksquare(A)))
 
-function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
+function \{TA,TB,N}(A::AbstractMatrix{TA}, B::AbstractArray{TB,N})
     TC = typeof(one(TA)/one(TB))
-    A_ldiv_B!(convert(typeof(A).name.primary{TC}, A), TB == TC ? copy(B) : convert(typeof(B).name.primary{TC}, B))
+    A_ldiv_B!(TA == TC ? copy(A) : convert(AbstractMatrix{TC}, A), TB == TC ? copy(B) : convert(AbstractArray{TC,N}, B))
 end
 \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
 /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'
diff --git a/base/linalg/ldlt.jl b/base/linalg/ldlt.jl
new file mode 100644
index 0000000000000..1bae952fca5f0
--- /dev/null
+++ b/base/linalg/ldlt.jl
@@ -0,0 +1,54 @@
+immutable LDLt{T,S<:AbstractMatrix{T}} <: Factorization{T}
+    data::S
+end
+
+size(S::LDLt) = size(S.data)
+size(S::LDLt, i::Integer) = size(S.data, i)
+
+convert{T,S}(::Type{LDLt{T,S}}, F::LDLt) = LDLt{T,S}(convert(S, F.data))
+convert{T,S,U}(::Type{LDLt{T}}, F::LDLt{S,U}) = convert(LDLt{T,U}, F)
+convert{T,S,U}(::Type{Factorization{T}}, F::LDLt{S,U}) = convert(LDLt{T,U}, F)
+
+# SymTridiagonal
+function ldltfact!{T<:Real}(S::SymTridiagonal{T})
+	n = size(S,1)
+	d = S.dv
+	e = S.ev
+	@inbounds for i = 1:n-1
+		e[i] /= d[i]
+		d[i+1] -= abs2(e[i])*d[i]
+		d[i+1] > 0 || throw(PosDefException(i+1))
+	end
+	return LDLt{T,SymTridiagonal{T}}(S)
+end
+function ldltfact{T}(M::SymTridiagonal{T})
+	S = typeof(one(T)/one(T))
+	return S == T ? ldltfact!(copy(M)) : ldltfact!(convert(SymTridiagonal{S}, M))
+end
+function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
+	n, nrhs = size(B, 1), size(B, 2)
+	size(S,1) == n || throw(DimensionMismatch(""))
+	d = S.data.dv
+	l = S.data.ev
+	@inbounds begin
+		for i = 2:n
+			li1 = l[i-1]
+			for j = 1:nrhs
+				B[i,j] -= li1*B[i-1,j]
+			end
+		end
+		dn = d[n]
+		for j = 1:nrhs
+			B[n,j] /= dn
+		end
+		for i = n-1:-1:1
+			di = d[i]
+			li = l[i]
+			for j = 1:nrhs
+				B[i,j] /= di
+				B[i,j] -= li*B[i+1,j]
+			end
+		end
+	end
+	return B
+end
\ No newline at end of file
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
new file mode 100644
index 0000000000000..44c4c5a7365fb
--- /dev/null
+++ b/base/linalg/lu.jl
@@ -0,0 +1,313 @@
+####################
+# LU Factorization #
+####################
+immutable LU{T,S<:AbstractMatrix{T}} <: Factorization{T}
+    factors::S
+    ipiv::Vector{BlasInt}
+    info::BlasInt
+end
+
+# StridedMatrix
+function lufact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot = true)
+    !pivot && return generic_lufact!(A, pivot=pivot)
+    lpt = LAPACK.getrf!(A)
+    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
+end
+lufact!(A::StridedMatrix; pivot = true) = generic_lufact!(A, pivot=pivot)
+function generic_lufact!{T}(A::StridedMatrix{T}; pivot = true)
+    m, n = size(A)
+    minmn = min(m,n)
+    info = 0
+    ipiv = Array(BlasInt, minmn)
+    @inbounds begin
+        for k = 1:minmn
+            # find index max
+            kp = k
+            if pivot
+                amax = real(zero(T))
+                for i = k:m
+                    absi = abs(A[i,k])
+                    if absi > amax
+                        kp = i
+                        amax = absi
+                    end
+                end
+            end
+            ipiv[k] = kp
+            if A[kp,k] != 0
+                if k != kp
+                    # Interchange
+                    for i = 1:n
+                        tmp = A[k,i]
+                        A[k,i] = A[kp,i]
+                        A[kp,i] = tmp
+                    end
+                end
+                # Scale first column
+                Akkinv = inv(A[k,k])
+                for i = k+1:m
+                    A[i,k] *= Akkinv
+                end
+            elseif info == 0
+                info = k
+            end
+            # Update the rest
+            for j = k+1:n
+                for i = k+1:m
+                    A[i,j] -= A[i,k]*A[k,j]
+                end
+            end        
+        end
+    end
+    LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
+end
+lufact{T<:BlasFloat}(A::AbstractMatrix{T}; pivot = true) = lufact!(copy(A), pivot=pivot)
+lufact{T}(A::AbstractMatrix{T}; pivot = true) = (S = typeof(one(T)/one(T)); S != T ? lufact!(convert(AbstractMatrix{S}, A), pivot=pivot) : lufact!(copy(A), pivot=pivot))
+lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt))
+
+function lu(A::Union(Number, AbstractMatrix))
+    F = lufact(A)
+    F[:L], F[:U], F[:p]
+end
+
+function convert{T}(::Type{LU{T}}, F::LU)
+    M = convert(AbstractMatrix{T}, F.factors)
+    LU{T,typeof(M)}(M, F.ipiv, F.info)
+end    
+convert{T,S}(::Type{LU{T,S}}, F::LU) = LU{T,S}(convert(S, F.factors), F.ipiv, F.info)
+convert{T}(::Type{Factorization{T}}, F::LU) = convert(LU{T}, F)
+
+
+size(A::LU) = size(A.factors)
+size(A::LU,n) = size(A.factors,n)
+
+function ipiv2perm{T}(v::AbstractVector{T}, maxi::Integer)
+    p = T[1:maxi]
+    @inbounds for i in 1:length(v)
+        p[i], p[v[i]] = p[v[i]], p[i]
+    end
+    return p
+end
+
+function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol)
+    m, n = size(A)
+    d == :L && return tril(A.factors[1:m, 1:min(m,n)], -1) + eye(T, m, min(m,n))
+    d == :U && return triu(A.factors[1:min(m,n),1:n])
+    d == :p && return ipiv2perm(A.ipiv, m)
+    if d == :P
+        p = A[:p]
+        P = zeros(T, m, m)
+        for i in 1:m
+            P[i,p[i]] = one(T)
+        end
+        return P
+    end
+    throw(KeyError(d))
+end
+
+A_ldiv_B!{T<:BlasFloat, S<:StridedMatrix}(A::LU{T, S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info
+A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, b::StridedVector) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), b[ipiv2perm(A.ipiv, length(b))]))
+A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, B::StridedMatrix) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), B[ipiv2perm(A.ipiv, size(B, 1)),:]))
+At_ldiv_B{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, copy(B)) A.info
+Ac_ldiv_B{T<:BlasComplex,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, copy(B)) A.info
+At_ldiv_Bt{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info
+Ac_ldiv_Bc{T<:BlasComplex,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info
+
+function det{T,S}(A::LU{T,S})
+    n = chksquare(A)
+    A.info > 0 && return zero(typeof(A.factors[1]))
+    return prod(diag(A.factors)) * (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
+end
+
+function logdet2{T<:Real,S}(A::LU{T,S})  # return log(abs(det)) and sign(det)
+    n = chksquare(A)
+    dg = diag(A.factors)
+    s = (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
+    sum(log(abs(dg))), s 
+end
+
+function logdet{T<:Real,S}(A::LU{T,S})
+    d,s = logdet2(A)
+    s>=0 || error("DomainError: determinant is negative")
+    d
+end
+
+function logdet{T<:Complex,S}(A::LU{T,S})
+    n = chksquare(A)
+    s = sum(log(diag(A.factors))) + (bool(sum(A.ipiv .!= 1:n) % 2) ? complex(0,pi) : 0) 
+    r, a = reim(s)
+    a = pi-mod(pi-a,2pi) #Take principal branch with argument (-pi,pi] 
+    complex(r,a)    
+end
+
+inv{T<:BlasFloat}(A::LU{T,StridedMatrix{T}}) = @assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
+
+cond{T<:BlasFloat}(A::LU{T,StridedMatrix{T}}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
+
+# Tridiagonal
+
+# See dgttrf.f
+function lufact!{T}(A::Tridiagonal{T}; pivot = true)
+    n = size(A, 1)
+    info = 0
+    ipiv = Array(BlasInt, n)
+    dl = A.dl
+    d = A.d
+    du = A.du
+    du2 = A.du2
+
+    @inbounds begin
+        for i = 1:n
+            ipiv[i] = i
+        end
+        for i = 1:n-2
+            # pivot or not?
+            if !pivot || abs(d[i]) >= abs(dl[i])
+                # No interchange
+                if d[i] != 0
+                    fact = dl[i]/d[i]
+                    dl[i] = fact
+                    d[i+1] -= fact*du[i]
+                    du2[i] = 0
+                end
+            else
+                # Interchange
+                fact = d[i]/dl[i]
+                d[i] = dl[i]
+                dl[i] = fact
+                tmp = du[i]
+                du[i] = d[i+1]
+                d[i+1] = tmp - fact*d[i+1]
+                du2[i] = du[i+1]
+                du[i+1] = -fact*du[i+1]
+                ipiv[i] = i+1
+            end
+        end
+        if n > 1
+            i = n-1
+            if !pivot || abs(d[i]) >= abs(dl[i])
+                if d[i] != 0
+                    fact = dl[i]/d[i]
+                    dl[i] = fact
+                    d[i+1] -= fact*du[i]
+                end
+            else
+                fact = d[i]/dl[i]
+                d[i] = dl[i]
+                dl[i] = fact
+                tmp = du[i]
+                du[i] = d[i+1]
+                d[i+1] = tmp - fact*d[i+1]
+                ipiv[i] = i+1
+            end
+        end
+        # check for a zero on the diagonal of U
+        for i = 1:n
+            if d[i] == 0
+                info = i
+                break
+            end
+        end
+    end
+    LU{T,Tridiagonal{T}}(A, ipiv, convert(BlasInt, info))
+end
+factorize(A::Tridiagonal) = lufact(A)
+
+# See dgtts2.f
+function A_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
+    n = size(A,1)
+    n == size(B,1) || throw(DimentionsMismatch(""))
+    nrhs = size(B,2)
+    dl = A.factors.dl
+    d = A.factors.d
+    du = A.factors.du
+    du2 = A.factors.du2
+    ipiv = A.ipiv
+    @inbounds begin
+        for j = 1:nrhs
+            for i = 1:n-1
+                ip = ipiv[i]
+                tmp = B[i+1-ip+i,j] - dl[i]*B[ip,j]
+                B[i,j] = B[ip,j]
+                B[i+1,j] = tmp
+            end
+            B[n,j] /= d[n]
+            if n > 1
+                B[n-1,j] = (B[n-1,j] - du[n-1]*B[n,j])/d[n-1]
+            end
+            for i = n-2:-1:1
+                B[i,j] = (B[i,j] - du[i]*B[i+1,j] - du2[i]*B[i+2,j])/d[i]
+            end
+        end
+    end
+    return B
+end
+
+function At_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
+    n = size(A,1)
+    n == size(B,1) || throw(DimentionsMismatch(""))
+    nrhs = size(B,2)
+    dl = A.factors.dl
+    d = A.factors.d
+    du = A.factors.du
+    du2 = A.factors.du2
+    ipiv = A.ipiv
+    @inbounds begin
+        for j = 1:nrhs
+            B[1,j] /= d[1]
+            if n > 1
+                B[2,j] = (B[2,j] - du[1]*B[1,j])/d[2]
+            end
+            for i = 3:n
+                B[i,j] = (B[i,j] - du[i-1]*B[i-1,j] - du2[i-2]*B[i-2,j])/d[i]
+            end
+            for i = n-1:-1:1
+                if ipiv[i] == i
+                    B[i,j] = B[i,j] - dl[i]*B[i+1,j]
+                else
+                    tmp = B[i+1,j]
+                    B[i+1,j] = B[i,j] - dl[i]*tmp
+                    B[i,j] = tmp
+                end
+            end
+        end
+    end
+    return B
+end
+
+# Ac_ldiv_B!{T<:Real}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) = At_ldiv_B!(A,B)
+function Ac_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
+    n = size(A,1)
+    n == size(B,1) || throw(DimentionsMismatch(""))
+    nrhs = size(B,2)
+    dl = A.factors.dl
+    d = A.factors.d
+    du = A.factors.du
+    du2 = A.factors.du2
+    ipiv = A.ipiv
+    @inbounds begin
+        for j = 1:nrhs
+            B[1,j] /= conj(d[1])
+            if n > 1
+                B[2,j] = (B[2,j] - conj(du[1])*B[1,j])/conj(d[2])
+            end
+            for i = 3:n
+                B[i,j] = (B[i,j] - conj(du[i-1])*B[i-1,j] - conj(du2[i-2])*B[i-2,j])/conj(d[i])
+            end
+            for i = n-1:-1:1
+                if ipiv[i] == i
+                    B[i,j] = B[i,j] - conj(dl[i])*B[i+1,j]
+                else
+                    tmp = B[i+1,j]
+                    B[i+1,j] = B[i,j] - conj(dl[i])*tmp
+                    B[i,j] = tmp
+                end
+            end
+        end
+    end
+    return B
+end
+
+/(B::AbstractMatrix,A::LU) = At_ldiv_Bt(A,B).'
+
diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
index d1ba5ab940242..2e073c37eb0cb 100644
--- a/base/linalg/symmetric.jl
+++ b/base/linalg/symmetric.jl
@@ -16,7 +16,9 @@ getindex(A::HermOrSym, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? ge
 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{AbstractMatrix{T}}, A::Symmetric) = convert(Symmetric{T}, A)
 convert{T}(::Type{Hermitian{T}},A::Hermitian) = Hermitian(convert(Matrix{T},A.S),A.uplo)
+convert{T}(::Type{AbstractMatrix{T}}, A::Hermitian) = convert(Hermitian{T}, A)
 copy(A::Symmetric) = Symmetric(copy(A.S),A.uplo)
 copy(A::Hermitian) = Hermitian(copy(A.S),A.uplo)
 ishermitian(A::Hermitian) = true
@@ -38,13 +40,13 @@ factorize(A::HermOrSym) = bkfact(A.S, symbol(A.uplo), issym(A))
 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(one(T)/norm(one(T)))); S != T ? eigfact!(convert(typeof(A).name.primary{S}, A)) : eigfact!(copy(A)))
+eigfact{T}(A::HermOrSym{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(AbstractMatrix{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(one(T)/norm(one(T)))); S != T ? eigvals!(convert(typeof(A).name.primary{S}, A, l, h)) : eigvals!(copy(A), l, h))
+eigvals{T}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigvals!(convert(AbstractMatrix{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]
 
diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl
index 04f8e466c3770..beac0fe8c55c5 100644
--- a/base/linalg/tridiag.jl
+++ b/base/linalg/tridiag.jl
@@ -1,7 +1,7 @@
 #### Specialized matrix types ####
 
 ## Hermitian tridiagonal matrices
-type SymTridiagonal{T} <: AbstractMatrix{T}
+immutable SymTridiagonal{T} <: AbstractMatrix{T}
     dv::Vector{T}                        # diagonal
     ev::Vector{T}                        # subdiagonal
     function SymTridiagonal(dv::Vector{T}, ev::Vector{T})
@@ -10,7 +10,7 @@ type SymTridiagonal{T} <: AbstractMatrix{T}
     end
 end
 
-SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(copy(dv), copy(ev))
+SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(dv, ev)
 
 function SymTridiagonal{Td,Te}(dv::Vector{Td}, ev::Vector{Te})
     T = promote_type(Td,Te)
@@ -19,6 +19,8 @@ end
 
 SymTridiagonal(A::AbstractMatrix) = diag(A,1)==diag(A,-1)?SymTridiagonal(diag(A), diag(A,1)):throw(DimensionMismatch("matrix is not symmetric; cannot convert to SymTridiagonal"))
 full{T}(M::SymTridiagonal{T}) = convert(Matrix{T}, M)
+convert{T}(::Type{SymTridiagonal{T}}, S::SymTridiagonal) = SymTridiagonal(convert(Vector{T}, S.dv), convert(Vector{T}, S.ev))
+convert{T}(::Type{AbstractMatrix{T}}, S::SymTridiagonal) = SymTridiagonal(convert(Vector{T}, S.dv), convert(Vector{T}, S.ev))
 convert{T}(::Type{Matrix{T}}, M::SymTridiagonal{T})=diagm(M.dv)+diagm(M.ev,-1)+conj(diagm(M.ev,1))
 
 size(m::SymTridiagonal) = (length(m.dv), length(m.dv))
@@ -136,32 +138,21 @@ function getindex{T}(A::SymTridiagonal{T}, i::Integer, j::Integer)
 end
 
 ## Tridiagonal matrices ##
-type Tridiagonal{T} <: AbstractMatrix{T}
+immutable Tridiagonal{T} <: AbstractMatrix{T}
     dl::Vector{T}    # sub-diagonal
     d::Vector{T}     # diagonal
     du::Vector{T}    # sup-diagonal
-    dutmp::Vector{T} # scratch space for vector RHS solver, sup-diagonal
-    rhstmp::Vector{T}# scratch space, rhs
-
-    function Tridiagonal(N::Integer)
-        dutmp = Array(T, N-1)
-        rhstmp = Array(T, N)
-        new(dutmp, rhstmp, similar(dutmp), similar(dutmp), similar(rhstmp))
-    end
-
-    function Tridiagonal(dl::Vector{T}, d::Vector{T}, du::Vector{T})
-        N = length(d)
-        if (length(dl) != N-1 || length(du) != N-1)
-            error(string("Cannot make Tridiagonal from incompatible lengths of subdiagonal, diagonal and superdiagonal: (", length(dl), ", ", length(d), ", ", length(du),")"))
-        end
-        new(copy(dl), copy(d), copy(du), Array(T,N-1), Array(T,N))
+    du2::Vector{T}   # supsup-diagonal for pivoting
+end
+function Tridiagonal{T}(dl::Vector{T}, d::Vector{T}, du::Vector{T})
+    n = length(d)
+    if (length(dl) != n-1 || length(du) != n-1)
+        error(string("Cannot make Tridiagonal from incompatible lengths of subdiagonal, diagonal and superdiagonal: (", length(dl), ", ", length(d), ", ", length(du),")"))
     end
+    Tridiagonal(dl, d, du, zeros(T,n-2))
 end
-
-Tridiagonal{T}(dl::Vector{T}, d::Vector{T}, du::Vector{T}) = Tridiagonal{T}(dl, d, du)
-
 function Tridiagonal{Tl, Td, Tu}(dl::Vector{Tl}, d::Vector{Td}, du::Vector{Tu})
-    Tridiagonal(map(v->copy(convert(Vector{promote_type(Tl,Td,Tu)}, v)), (dl, d, du))...)
+    Tridiagonal(map(v->convert(Vector{promote_type(Tl,Td,Tu)}, v), (dl, d, du))...)
 end
 
 size(M::Tridiagonal) = (length(M.d), length(M.d))
@@ -183,16 +174,16 @@ function similar(M::Tridiagonal, T, dims::Dims)
     if length(dims) != 2 || dims[1] != dims[2]
         throw(DimensionMismatch("Tridiagonal matrices must be square"))
     end
-    Tridiagonal{T}(dims[1])
+    Tridiagonal{T}(similar(M.dl), similar(M.d), similar(M.du), similar(M.du2))
 end
 
 # Operations on Tridiagonal matrices
-copy!(dest::Tridiagonal, src::Tridiagonal) = Tridiagonal(copy!(dest.dl, src.dl), copy!(dest.d, src.d), copy!(dest.du, src.du))
+copy!(dest::Tridiagonal, src::Tridiagonal) = Tridiagonal(copy!(dest.dl, src.dl), copy!(dest.d, src.d), copy!(dest.du, src.du), copy!(dest.du2, src.du2))
 
 #Elementary operations
 for func in (:copy, :round, :iround, :conj) 
     @eval begin
-        ($func)(M::Tridiagonal) = Tridiagonal(map(($func), (M.dl, M.d, M.du))...)
+        ($func)(M::Tridiagonal) = Tridiagonal(map(($func), (M.dl, M.d, M.du, M.du2))...)
     end
 end
 
@@ -200,13 +191,11 @@ transpose(M::Tridiagonal) = Tridiagonal(M.du, M.d, M.dl)
 ctranspose(M::Tridiagonal) = conj(transpose(M))
 
 diag{T}(M::Tridiagonal{T}, n::Integer=0) = n==0 ? M.d : n==-1 ? M.dl : n==1 ? M.du : abs(n)<size(M,1) ? zeros(T,size(M,1)-abs(n)) : throw(BoundsError()) 
-
 function getindex{T}(A::Tridiagonal{T}, i::Integer, j::Integer)
     (1<=i<=size(A,2) && 1<=j<=size(A,2)) || throw(BoundsError())
     i==j ? A.d[i] : i==j+1 ? A.dl[j] : i+1==j ? A.du[i] : zero(T)
 end
 
-
 ###################
 # Generic methods #
 ###################
@@ -231,181 +220,32 @@ 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{AbstractMatrix{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::Tridiagonal) = Tridiagonal(convert(Vector{T}, M.dl), convert(Vector{T}, M.d), convert(Vector{T}, M.du), convert(Vector{T}, M.du2))
+convert{T}(::Type{AbstractMatrix{T}},M::Tridiagonal) = convert(Tridiagonal{T}, M)
 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
-
-#### Tridiagonal matrix routines ####
-function \{T<:BlasFloat}(M::Tridiagonal{T}, rhs::StridedVecOrMat{T})
-    if stride(rhs, 1) == 1
-        return LAPACK.gtsv!(copy(M.dl), copy(M.d), copy(M.du), copy(rhs))
-    end
-    solve(M, rhs)  # use the Julia "fallback"
-end
-
-# This is definitely not going to work
-#eig(M::Tridiagonal) = LAPACK.stev!('V', copy(M))
-
-# Allocation-free variants
-# Note that solve is non-aliasing, so you can use the same array for
-# input and output
-function solve!{T<:BlasFloat}(x::AbstractArray{T}, xrng::Range{Int}, M::Tridiagonal{T}, rhs::AbstractArray{T}, rhsrng::Range{Int})
-    d = M.d
-    N = length(d)
-    if length(xrng) != N || length(rhsrng) != N
-        throw(DimensionMismatch(""))
-    end
-    dl = M.dl
-    du = M.du
-    dutmp = M.dutmp
-    rhstmp = M.rhstmp
-    xstart = first(xrng)
-    xstride = step(xrng)
-    rhsstart = first(rhsrng)
-    rhsstride = step(rhsrng)
-    # Forward sweep
-    denom = d[1]
-    dulast = du[1] / denom
-    dutmp[1] = dulast
-    rhslast = rhs[rhsstart] / denom
-    rhstmp[1] = rhslast
-    irhs = rhsstart+rhsstride
-    for i in 2:N-1
-        dltmp = dl[i-1]
-        denom = d[i] - dltmp*dulast
-        dulast = du[i] / denom
-        dutmp[i] = dulast
-        rhslast = (rhs[irhs] - dltmp*rhslast)/denom
-        rhstmp[i] = rhslast
-        irhs += rhsstride
-    end
-    dltmp = dl[N-1]
-    denom = d[N] - dltmp*dulast
-    xlast = (rhs[irhs] - dltmp*rhslast)/denom
-    # Backward sweep
-    ix = xstart + (N-2)*xstride
-    x[ix+xstride] = xlast
-    for i in N-1:-1:1
-        xlast = rhstmp[i] - dutmp[i]*xlast
-        x[ix] = xlast
-        ix -= xstride
-    end
-    nothing
-end
-
-function solve!(x::StridedVector, M::Tridiagonal, rhs::StridedVector)
-    solve!(x, 1:length(x), M, rhs, 1:length(rhs))
-    x
-end
-solve{TM<:BlasFloat,TB<:BlasFloat}(M::Tridiagonal{TM}, B::StridedVecOrMat{TB}) = solve!(zeros(typeof(one(TM)/one(TB)), size(B)), M, B)
-solve(M::Tridiagonal, B::StridedVecOrMat) = solve(float(M), float(B))
-function solve!(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
-    size(B, 1) == size(M, 1) || throw(DimensionMismatch(""))
-    size(X) == size(B) || throw(DimensionMismatch(""))
-    m, n = size(B)
-    for j = 1:n
-        r = ((j-1)*m+1):((j-1)*m+m)
-        solve!(X, r, M, B, r)
-    end
-    X
-end
-
-# User-friendly solver
-\(M::Tridiagonal, rhs::StridedVecOrMat) = solve(M, rhs)
-
-# Tridiagonal multiplication
-function mult(x::AbstractArray, xrng::Range{Int}, M::Tridiagonal, v::AbstractArray, vrng::Range{Int})
-    dl = M.dl
-    d = M.d
-    du = M.du
-    N = length(d)
-    xi = first(xrng)
-    xstride = step(xrng)
-    vi = first(vrng)
-    vstride = step(vrng)
-    x[xi] = d[1]*v[vi] + du[1]*v[vi+vstride]
-    xi += xstride
-    for i = 2:N-1
-        x[xi] = dl[i-1]*v[vi] + d[i]*v[vi+vstride] + du[i]*v[vi+2*vstride]
-        xi += xstride
-        vi += vstride
-    end
-    x[xi] = dl[N-1]*v[vi] + d[N]*v[vi+vstride]
-    x
-end
-
-mult(x::StridedVector, M::Tridiagonal, v::StridedVector) = mult(x, 1:length(x), M, v, 1:length(v))
-
-function mult(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
-    size(B, 1) == size(M, 1) || throw(DimensionMismatch(""))
-    size(X) == size(B) || throw(DimensionMismatch(""))
-    m, n = size(B)
-    for j = 1:n
-        r = ((j-1)*m+1):((j-1)*m+m)
-        mult(X, r, M, B, r)
-    end
-    X
-end
-
-mult(X::StridedMatrix, M1::Tridiagonal, M2::Tridiagonal) = mult(X, M1, full(M2))
-
-*(M::Tridiagonal, B::Union(StridedVector,StridedMatrix)) = mult(similar(B), M, B)
-*(A::Tridiagonal, B::Tridiagonal) = A*full(B)
-
-#### Factorizations for Tridiagonal ####
-type LDLTTridiagonal{T<:BlasFloat,S<:BlasFloat} <: Factorization{T}
-    D::Vector{S}
-    E::Vector{T}
-    function LDLTTridiagonal(D::Vector{S}, E::Vector{T})
-        typeof(real(E[1])) == eltype(D) ? new(D, E) : error("element types do not match")
-        new(D, E)
+function A_mul_B!(C::AbstractVecOrMat, A::Tridiagonal, B::AbstractVecOrMat)
+    size(C,1) == size(B,1) == (nA = size(A,1)) || throw(DimensionMismatch(""))
+    size(C,2) == (nB = size(B,2)) || throw(DimensionMismatch(""))
+    l = A.dl
+    d = A.d
+    u = A.du
+    @inbounds begin
+        for j = 1:nB
+            C[1,j] = d[1]*B[1,j] + u[1]*B[2,j]
+            for i = 2:nA-1
+                C[i,j] = l[i-1]*B[i-1,j] + d[i]*B[i,j] + u[i]*B[i+1,j]
+            end
+            C[nA,j] = l[nA-1]*B[nA-1,j] + d[nA]*B[nA,j]
+        end
     end
+    C
 end
+*(A::Tridiagonal, B::AbstractVecOrMat) = A_mul_B!(similar(B), A, B)
 
-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<: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)
-A_ldiv_B!(C::LDLTTridiagonal, B::StridedVecOrMat) = A_ldiv_B!(C, float(B))
-
-type LUTridiagonal{T} <: Factorization{T}
-    dl::Vector{T}
-    d::Vector{T}
-    du::Vector{T}
-    du2::Vector{T}
-    ipiv::Vector{BlasInt}
-    # function LUTridiagonal(dl::Vector{T}, d::Vector{T}, du::Vector{T},
-    #                        du2::Vector{T}, ipiv::Vector{BlasInt})
-    #     n = length(d)
-    #     if length(dl) != n - 1 || length(du) != n - 1 || length(ipiv) != n || length(du2) != n-2
-    #         throw(DimensionMismatch("LUTridiagonal")
-    #     end
-    #     new(dl, d, du, du2, ipiv)
-    # end
-end
-lufact!{T<:BlasFloat}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
-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})
-    n = length(lu.d)
-    prod(lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
-end
-
-det{T<:BlasFloat}(A::Tridiagonal{T}) = det(lufact(A))
-
-A_ldiv_B!{T<:BlasFloat}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
-    LAPACK.gttrs!('N', lu.dl, lu.d, lu.du, lu.du2, lu.ipiv, B)
-A_ldiv_B!(lu::LUTridiagonal, B::StridedVecOrMat) = A_ldiv_B!(lu, float(B))
+A_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = A_ldiv_B!(lufact!(A), B)
+At_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = At_ldiv_B!(lufact!(A), B)
+Ac_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = Ac_ldiv_B!(lufact!(A), B)
diff --git a/base/linalg/umfpack.jl b/base/linalg/umfpack.jl
index ed4e5b45e9283..05140f7decaa0 100644
--- a/base/linalg/umfpack.jl
+++ b/base/linalg/umfpack.jl
@@ -8,7 +8,7 @@ export UmfpackLU,
 
 import Base: (\), Ac_ldiv_B, At_ldiv_B, findnz, getindex, show, size
 
-import ..LinAlg: A_ldiv_B!, Ac_ldiv_B!, At_ldiv_B!, Factorization, det, lufact, lufact!, solve
+import ..LinAlg: A_ldiv_B!, Ac_ldiv_B!, At_ldiv_B!, Factorization, det, lufact, lufact!
 
 include("umfpack_h.jl")
 type MatrixIllConditionedException <: Exception
diff --git a/base/linalg/woodbury.jl b/base/linalg/woodbury.jl
index 0a9b29c4a747b..99a51ec9493e6 100644
--- a/base/linalg/woodbury.jl
+++ b/base/linalg/woodbury.jl
@@ -76,38 +76,35 @@ end
 
 det(W::Woodbury)=det(W.A)*det(W.C)/det(W.Cp)
 
-# Allocation-free solver for arbitrary strides (requires that W.A has a
-# non-aliasing "solve" routine, e.g., is Tridiagonal)
-function solve!(x::AbstractArray, xrng::Range{Int}, W::Woodbury, rhs::AbstractArray, rhsrng::Range{Int})
-    solve!(W.tmpN1, 1:length(W.tmpN1), W.A, rhs, rhsrng)
+function A_ldiv_B!(W::Woodbury, B::AbstractArray)
+    copy!(W.tmpN1, B)
+    Alu = lufact(W.A) # Note. This makes an allocation. Alternative is to destroy W.A.
+    A_ldiv_B!(Alu, W.tmpN1)
     A_mul_B!(W.tmpk1, W.V, W.tmpN1)
     A_mul_B!(W.tmpk2, W.Cp, W.tmpk1)
     A_mul_B!(W.tmpN2, W.U, W.tmpk2)
-    solve!(W.tmpN2, W.A, W.tmpN2)
-    indx = first(xrng)
-    xinc = step(xrng)
+    A_ldiv_B!(Alu, W.tmpN2)
     for i = 1:length(W.tmpN2)
-        x[indx] = W.tmpN1[i] - W.tmpN2[i]
-        indx += xinc
+        B[i] = W.tmpN1[i] - W.tmpN2[i]
     end
-    nothing
+    B
 end
 
-function solve!(x::AbstractVector, W::Woodbury, rhs::AbstractVector)
-    solve!(x, 1:length(x), W, rhs, 1:length(rhs))
-    x
-end
-solve(W::Woodbury, rhs::AbstractVector) = solve!(similar(rhs), W, rhs)
-solve(W::Woodbury, B::StridedMatrix)=solve!(similar(B), W, B)
-function solve!(X::StridedMatrix, W::Woodbury, B::StridedMatrix)
-    size(B, 1) == size(W, 1) || throw(DimensionMismatch(""))
-    size(X) == size(B) || throw(DimensionMismatch(""))
-    m, n = size(B)
-    r = 1:m
-    for j = 1:n
-        r.start = (j-1)*m+1
-        solve!(X, r, W, B, r)
-    end
-    X
-end
+# function solve!(x::AbstractVector, W::Woodbury, rhs::AbstractVector)
+#     solve!(x, 1:length(x), W, rhs, 1:length(rhs))
+#     x
+# end
+# solve(W::Woodbury, rhs::AbstractVector) = solve!(similar(rhs), W, rhs)
+# solve(W::Woodbury, B::StridedMatrix)=solve!(similar(B), W, B)
+# function solve!(X::StridedMatrix, W::Woodbury, B::StridedMatrix)
+#     size(B, 1) == size(W, 1) || throw(DimensionMismatch(""))
+#     size(X) == size(B) || throw(DimensionMismatch(""))
+#     m, n = size(B)
+#     r = 1:m
+#     for j = 1:n
+#         r.start = (j-1)*m+1
+#         solve!(X, r, W, B, r)
+#     end
+#     X
+# end
 
diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl
index f13c507787315..088d16732c61d 100644
--- a/base/sparse/sparsematrix.jl
+++ b/base/sparse/sparsematrix.jl
@@ -162,7 +162,7 @@ function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, M::Matrix)
                              convert(Vector{Tv},V), 
                              m, n)
 end
-
+convert{T}(::Type{AbstractMatrix{T}}, A::SparseMatrixCSC) = convert(SparseMatrixCSC{T}, A)
 convert(::Type{Matrix}, S::SparseMatrixCSC) = full(S)
 
 function full{Tv}(S::SparseMatrixCSC{Tv})
diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst
index ed21df40112d7..63808f1ca7352 100644
--- a/doc/stdlib/linalg.rst
+++ b/doc/stdlib/linalg.rst
@@ -43,41 +43,41 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Compute the LU factorization of ``A``, such that ``A[p,:] = L*U``.
 
-.. function:: lufact(A) -> F
+.. function:: lufact(A, [pivot=true]) -> F
 
-   Compute the LU factorization of ``A``. The return type of ``F`` depends on the type of ``A``.
+   Compute the LU factorization of ``A``. The return type of ``F`` depends on the type of ``A``. In most cases, if ``A`` is a subtype ``S`` of AbstractMatrix with an element type ``T``` supporting ``+``, ``-``, ``*`` and ``/`` the return type is ``LU{T,S{T}}``. If pivoting is chosen (default) the element type should also support ``abs`` and ``<``. When ``A`` is sparse and have element of type ``Float32``, ``Float64``, ``Complex{Float32}``, or ``Complex{Float64}`` the return type is ``UmfpackLU``. Some examples are shown in the table below.
 
-      ======================= ==================== ========================================
-      Type of input ``A``     Type of output ``F`` Relationship between ``F`` and ``A``
-      ----------------------- -------------------- ----------------------------------------
-      :func:`Matrix`           ``LU``              ``F[:L]*F[:U] == A[F[:p], :]``
-      :func:`Tridiagonal`      ``LUTridiagonal``     N/A
-      :func:`SparseMatrixCSC`  ``UmfpackLU``       ``F[:L]*F[:U] == Rs .* A[F[:p], F[:q]]``
-      ======================= ==================== ========================================
+      ======================= ======================== ========================================
+      Type of input ``A``     Type of output ``F``     Relationship between ``F`` and ``A``
+      ----------------------- ------------------------ ----------------------------------------
+      :func:`Matrix`           ``LU``                  ``F[:L]*F[:U] == A[F[:p], :]``
+      :func:`Tridiagonal`      ``LU{T,Tridiagonal{T}}``  N/A
+      :func:`SparseMatrixCSC`  ``UmfpackLU``           ``F[:L]*F[:U] == Rs .* A[F[:p], F[:q]]``
+      ======================= ======================== ========================================
 
    The individual components of the factorization ``F`` can be accessed by indexing:
 
-      =========== ======================================= ====== ================= =============
-      Component   Description                             ``LU`` ``LUTridiagonal`` ``UmfpackLU``
-      ----------- --------------------------------------- ------ ----------------- -------------
-      ``F[:L]``   ``L`` (lower triangular) part of ``LU``    ✓                      ✓
-      ``F[:U]``   ``U`` (upper triangular) part of ``LU``    ✓                      ✓
-      ``F[:p]``   (right) permutation ``Vector``             ✓                      ✓
-      ``F[:P]``   (right) permutation ``Matrix``             ✓
-      ``F[:q]``   left permutation ``Vector``                                       ✓
-      ``F[:Rs]``  ``Vector`` of scaling factors                                     ✓
-      ``F[:(:)]`` ``(L,U,p,q,Rs)`` components                                       ✓
-      =========== ======================================= ====== ================= =============
-
-      ================== ====== ================= =============
-      Supported function ``LU`` ``LUTridiagonal`` ``UmfpackLU``
-      ------------------ ------ ----------------- -------------
+      =========== ======================================= ====== ======================= ==============
+      Component   Description                             ``LU`` ``LU{T,Tridiagonal{T}}`` ``UmfpackLU``
+      ----------- --------------------------------------- ------ ----------------------- -------------
+      ``F[:L]``   ``L`` (lower triangular) part of ``LU``    ✓                                     ✓
+      ``F[:U]``   ``U`` (upper triangular) part of ``LU``    ✓                                     ✓
+      ``F[:p]``   (right) permutation ``Vector``             ✓                                     ✓
+      ``F[:P]``   (right) permutation ``Matrix``             ✓              
+      ``F[:q]``   left permutation ``Vector``                                                      ✓
+      ``F[:Rs]``  ``Vector`` of scaling factors                                                    ✓
+      ``F[:(:)]`` ``(L,U,p,q,Rs)`` components                                                      ✓
+      =========== ======================================= ====== ======================= =============
+
+      ================== ====== ======================= =============
+      Supported function ``LU`` ``LU{T,Tridiagonal{T}}`` ``UmfpackLU``
+      ------------------ ------ ----------------------- -------------
            ``/``            ✓
-           ``\``            ✓        ✓             ✓
-           ``cond``         ✓                      ✓
-           ``det``          ✓                      ✓
-           ``size``         ✓
-      ================== ====== ================= =============
+           ``\``            ✓                       ✓             ✓
+           ``cond``         ✓                                     ✓
+           ``det``          ✓                       ✓             ✓
+           ``size``         ✓                       ✓
+      ================== ====== ======================= =============
 
 .. function:: lufact!(A) -> LU
 
@@ -93,12 +93,16 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: cholfact(A, [ll]) -> CholmodFactor
 
-   Compute the sparse Cholesky factorization of a sparse matrix ``A``.  If ``A`` is Hermitian its Cholesky factor is determined.  If ``A`` is not Hermitian the Cholesky factor of ``A*A'`` is determined. A fill-reducing permutation is used.  Methods for ``size``, ``solve``, ``\``, ``findn_nzs``, ``diag``, ``det`` and ``logdet``.  One of the solve methods includes an integer argument that can be used to solve systems involving parts of the factorization only.  The optional boolean argument, ``ll`` determines whether the factorization returned is of the ``A[p,p] = L*L'`` form, where ``L`` is lower triangular or ``A[p,p] = scale(L,D)*L'`` form where ``L`` is unit lower triangular and ``D`` is a non-negative vector.  The default is LDL.
+   Compute the sparse Cholesky factorization of a sparse matrix ``A``.  If ``A`` is Hermitian its Cholesky factor is determined.  If ``A`` is not Hermitian the Cholesky factor of ``A*A'`` is determined. A fill-reducing permutation is used.  Methods for ``size``, ``solve``, ``\``, ``findn_nzs``, ``diag``, ``det`` and ``logdet``.  One of the solve methods includes an integer argument that can be used to solve systems involving parts of the factorization only.  The optional boolean argument, ``ll`` determines whether the factorization returned is of the ``A[p,p] = L*L'`` form, where ``L`` is lower triangular or ``A[p,p] = L*Diagonal(D)*L'`` form where ``L`` is unit lower triangular and ``D`` is a non-negative vector.  The default is LDL.
 
 .. function:: cholfact!(A, [LU,][pivot=false,][tol=-1.0]) -> Cholesky
 
    ``cholfact!`` is the same as :func:`cholfact`, but saves space by overwriting the input ``A``, instead of creating a copy.
 
+.. function:: ldltfact(A) -> LDLtFactorization
+
+   Compute a factorization of a positive definite matrix ``A`` such that ``A=L*Diagonal(d)*L'`` where ``L`` is a unit lower triangular matrix and ``d`` is a vector with non-negative elements.
+
 .. function:: qr(A, [pivot=false,][thin=true]) -> Q, R, [p]
 
    Compute the (pivoted) QR factorization of ``A`` such that either ``A = Q*R`` or ``A[:,p] = Q*R``. Also see ``qrfact``. The default is to compute a thin factorization. Note that ``R`` is not extended with zeros when the full ``Q`` is requested. 
diff --git a/test/linalg2.jl b/test/linalg2.jl
index dc8ce5e028dc7..46198fe43fe16 100644
--- a/test/linalg2.jl
+++ b/test/linalg2.jl
@@ -68,20 +68,22 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
     @test_approx_eq T*v F*v
     invFv = F\v
     @test_approx_eq T\v invFv
-    @test_approx_eq Base.solve(T,v) invFv
-    @test_approx_eq Base.solve(T, B) F\B
+    # @test_approx_eq Base.solve(T,v) invFv
+    # @test_approx_eq Base.solve(T, B) F\B
     Tlu = factorize(T)
     x = Tlu\v
     @test_approx_eq x invFv
     @test_approx_eq det(T) det(F)
 
     # symmetric tridiagonal
-    Ts = SymTridiagonal(d, dl)
-    Fs = full(Ts)
-    invFsv = Fs\v
-    Tldlt = Base.ldltd(Ts)
-    x = Tldlt\v
-    @test_approx_eq x invFsv
+    if elty <: Real
+        Ts = SymTridiagonal(d, dl)
+        Fs = full(Ts)
+        invFsv = Fs\v
+        Tldlt = ldltfact(Ts)
+        x = Tldlt\v
+        @test_approx_eq x invFsv
+    end
 
     # eigenvalues/eigenvectors of symmetric tridiagonal
     if elty === Float32 || elty === Float64
@@ -100,7 +102,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
     @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
     iWv = similar(iFv)
     if elty != Int
-        Base.LinAlg.solve!(iWv, W, v)
+        iWv = A_ldiv_B!(W, copy(v))
         @test_approx_eq iWv iFv
     end
 
@@ -477,7 +479,7 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
             A = convert(Matrix{elty}, complex(randn(10,nn),randn(10,nn)))
         end    ## LU (only equal for real because LAPACK uses different absolute value when choosing permutations)
         if elty <: Real
-            FJulia  = invoke(lufact!, (AbstractMatrix,), copy(A)) 
+            FJulia  = Base.LinAlg.generic_lufact!(copy(A))
             FLAPACK = Base.LinAlg.LAPACK.getrf!(copy(A))
             @test_approx_eq FJulia.factors FLAPACK[1]
             @test_approx_eq FJulia.ipiv FLAPACK[2]
diff --git a/test/suitesparse.jl b/test/suitesparse.jl
index 1fe3d5f1c3e1d..c2b8d3a8760bf 100644
--- a/test/suitesparse.jl
+++ b/test/suitesparse.jl
@@ -149,7 +149,7 @@ afiro = CholmodSparse!(int32([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,
                         -1.0,2.279,1.4,-1.0,1.0,-1.0,1.0,1.0,1.0], 27, 51, 0)
 chmaf = cholfact(afiro)
 y = afiro'*ones(size(afiro,1))
-sol = Base.solve(chmaf, afiro*y) # least squares solution
+sol = chmaf\(afiro*y) # least squares solution
 @test isvalid(sol)
 pred = afiro'*sol
 @test norm(afiro * (y.mat - pred.mat)) < 1e-8