diff --git a/NEWS.md b/NEWS.md
index 967e03a48513d..3e4a15ceaf01b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -154,6 +154,8 @@ Library improvements
 
   * copy(DArray) will now make a copy of the DArray ([#9745])
 
+  * Split `Triangular` type into `UpperTriangular`, `LowerTriangular`, `UnitUpperTriagular` and `UnitLowerTriangular` ([#9779])
+
 Deprecated or removed
 ---------------------
 
diff --git a/base/exports.jl b/base/exports.jl
index e319785717794..3b257f3e38a94 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -62,6 +62,7 @@ export
     IOBuffer,
     IOStream,
     LocalProcess,
+    LowerTriangular,
     MathConst,
     Matrix,
     MergeSort,
@@ -108,9 +109,9 @@ export
     SymTridiagonal,
     Timer,
     TmStruct,
-    Triangular,
     Tridiagonal,
     UnitRange,
+    UpperTriangular,
     UTF16String,
     UTF32String,
     Val,
diff --git a/base/linalg.jl b/base/linalg.jl
index cfc7c74eca088..e81942bb55c1c 100644
--- a/base/linalg.jl
+++ b/base/linalg.jl
@@ -30,7 +30,8 @@ export
     SVD,
     Hermitian,
     Symmetric,
-    Triangular,
+    LowerTriangular,
+    UpperTriangular,
     Diagonal,
     UniformScaling,
 
diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl
index b804a5f03a8fe..8701c11b8bdf5 100644
--- a/base/linalg/bidiag.jl
+++ b/base/linalg/bidiag.jl
@@ -114,24 +114,23 @@ end
 /(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper)
 ==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)
 
-SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular)
+SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, TriangularUnion)
 *(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B)
 
 #Generic multiplication
 for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
     @eval begin
         ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B)
-        #($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B)
     end
 end
 
 
 #Linear solvers
-A_ldiv_B!(A::Union(Bidiagonal, Triangular), b::AbstractVector) = naivesub!(A, b)
-At_ldiv_B!(A::Union(Bidiagonal, Triangular), b::AbstractVector) = naivesub!(transpose(A), b)
-Ac_ldiv_B!(A::Union(Bidiagonal, Triangular), b::AbstractVector) = naivesub!(ctranspose(A), b)
+A_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(A, b)
+At_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(transpose(A), b)
+Ac_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(ctranspose(A), b)
 for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
-    function ($func)(A::Union(Bidiagonal, Triangular), B::AbstractMatrix)
+    function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix)
         tmp = similar(B[:,1])
         n = size(B, 1)
         for i = 1:size(B,2)
@@ -143,7 +142,7 @@ for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
     end
 end end
 for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
-    function ($func)(A::Union(Bidiagonal, Triangular), B::AbstractMatrix)
+    function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix)
         tmp = similar(B[:, 2])
         m, n = size(B)
         nm = n*m
diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl
index 4ee7e30e17c4e..1fca762bf4610 100644
--- a/base/linalg/cholesky.jl
+++ b/base/linalg/cholesky.jl
@@ -18,12 +18,35 @@ immutable CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T}
 end
 CholeskyPivoted{T}(UL::AbstractMatrix{T}, uplo::Char, piv::Vector{BlasInt}, rank::BlasInt, tol::Real, info::BlasInt) = CholeskyPivoted{T,typeof(UL)}(UL, uplo, piv, rank, tol, info)
 
-function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U)
+function chol!{T<:BlasFloat}(A::StridedMatrix{T})
+    C, info = LAPACK.potrf!('U', A)
+    return @assertposdef UpperTriangular(C) info
+end
+function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol)
     C, info = LAPACK.potrf!(char_uplo(uplo), A)
-    return @assertposdef Triangular{eltype(C),typeof(C),uplo,false}(C) info
+    return @assertposdef uplo == :U ? UpperTriangular(C) : LowerTriangular(C) info
 end
 
-function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U)
+function chol!{T}(A::AbstractMatrix{T})
+    n = chksquare(A)
+    @inbounds begin
+        for k = 1:n
+            for i = 1:k - 1
+                A[k,k] -= A[i,k]'A[i,k]
+            end
+            A[k,k] = chol!(A[k,k], uplo)
+            AkkInv = inv(A[k,k])
+            for j = k + 1:n
+                for i = 1:k - 1
+                    A[k,j] -= A[i,k]'A[i,j]
+                end
+                A[k,j] = A[k,k]'\A[k,j]
+            end
+        end
+    end
+    return UpperTriangular(A)
+end
+function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol)
     n = chksquare(A)
     @inbounds begin
         if uplo == :L
@@ -58,7 +81,7 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U)
             throw(ArgumentError("uplo must be either :U or :L but was $(uplo)"))
         end
     end
-    return Triangular(A, uplo, false)
+    return uplo == :U ? UpperTriangular(A) : LowerTriangular(A)
 end
 
 function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
@@ -118,14 +141,14 @@ size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
 size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)
 
 function getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, d::Symbol)
-    d == :U && return Triangular(UpLo == d ? C.UL : C.UL',:U)
-    d == :L && return Triangular(UpLo == d ? C.UL : C.UL',:L)
-    d == :UL && return Triangular(C.UL, UpLo)
+    d == :U && return UpperTriangular(UpLo == d ? C.UL : C.UL')
+    d == :L && return LowerTriangular(UpLo == d ? C.UL : C.UL')
+    d == :UL && return UpLo == :U ? UpperTriangular(C.UL) : LowerTriangular(C.UL)
     throw(KeyError(d))
 end
 function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol)
-    d == :U && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :U)
-    d == :L && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :L)
+    d == :U && return UpperTriangular(symbol(C.uplo) == d ? C.UL : C.UL')
+    d == :L && return LowerTriangular(symbol(C.uplo) == d ? C.UL : C.UL')
     d == :p && return C.piv
     if d == :P
         n = size(C, 1)
@@ -142,8 +165,8 @@ show{T,S<:AbstractMatrix,UpLo}(io::IO, C::Cholesky{T,S,UpLo}) = (println("$(type
 
 A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat{T}) = LAPACK.potrs!('U', C.UL, B)
 A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat{T}) = LAPACK.potrs!('L', C.UL, B)
-A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat) = Ac_ldiv_B!(Triangular(C.UL, :L, false), A_ldiv_B!(Triangular(C.UL, :L, false), B))
-A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat) = A_ldiv_B!(Triangular(C.UL, :U, false), Ac_ldiv_B!(Triangular(C.UL, :U, false), B))
+A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat) = Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B))
+A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat) = A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B))
 
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T})
     chkfullrank(C)
@@ -161,8 +184,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, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)]
-A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:]
+A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B[C.piv]))[invperm(C.piv)]
+A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(LowerTriangular(C.UL), A_ldiv_B!(LowerTriangular(C.UL), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(UpperTriangular(C.UL), Ac_ldiv_B!(UpperTriangular(C.UL), B[C.piv,:]))[invperm(C.piv),:]
 
 function det{T,S,UpLo}(C::Cholesky{T,S,UpLo})
     dd = one(T)
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index 2a1966753c5e4..5477c91a03afb 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -307,7 +307,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T})
     issym(A) && return sqrtm(Symmetric(A))
     n = chksquare(A)
     SchurF = schurfact(complex(A))
-    R = full(sqrtm(Triangular(SchurF[:T], :U, false)))
+    R = full(sqrtm(UpperTriangular(SchurF[:T])))
     retmat = SchurF[:vectors]*R*SchurF[:vectors]'
     all(imag(retmat) .== 0) ? real(retmat) : retmat
 end
@@ -315,7 +315,7 @@ 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], :U, false)))
+    R = full(sqrtm(UpperTriangular(SchurF[:T])))
     SchurF[:vectors]*R*SchurF[:vectors]'
 end
 sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
@@ -325,9 +325,9 @@ function inv{S}(A::StridedMatrix{S})
     T = typeof(one(S)/one(S))
     Ac = convert(AbstractMatrix{T}, A)
     if istriu(Ac)
-        Ai = inv(Triangular(A, :U, false))
+        Ai = inv(UpperTriangular(A))
     elseif istril(Ac)
-        Ai = inv(Triangular(A, :L, false))
+        Ai = inv(LowerTriangular(A))
     else
         Ai = inv(lufact(Ac))
     end
@@ -377,7 +377,7 @@ function factorize{T}(A::Matrix{T})
                 if utri1
                     return Bidiagonal(diag(A), diag(A, -1), false)
                 end
-                return Triangular(A, :L)
+                return LowerTriangular(A)
             end
             if utri
                 return Bidiagonal(diag(A), diag(A, 1), true)
@@ -392,7 +392,7 @@ function factorize{T}(A::Matrix{T})
             end
         end
         if utri
-            return Triangular(A, :U)
+            return UpperTriangular(A)
         end
         if herm
             try
@@ -413,9 +413,9 @@ function (\)(A::StridedMatrix, B::StridedVecOrMat)
     m, n = size(A)
     if m == n
         if istril(A)
-            return istriu(A) ? \(Diagonal(A),B) : \(Triangular(A, :L),B)
+            return istriu(A) ? \(Diagonal(A),B) : \(LowerTriangular(A),B)
         end
-        istriu(A) && return \(Triangular(A, :U),B)
+        istriu(A) && return \(UpperTriangular(A),B)
         return \(lufact(A),B)
     end
     return qrfact(A,pivot=eltype(A)<:BlasFloat)\B
diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl
index 7a621e8a73ec5..d2c52e328271b 100644
--- a/base/linalg/diagonal.jl
+++ b/base/linalg/diagonal.jl
@@ -8,7 +8,8 @@ Diagonal(A::Matrix) = Diagonal(diag(A))
 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)
-convert{T}(::Type{Triangular}, A::Diagonal{T}) = Triangular{T, Diagonal{T}, :U, false}(A)
+convert{T}(::Type{UpperTriangular}, A::Diagonal{T}) = UpperTriangular(A)
+convert{T}(::Type{LowerTriangular}, A::Diagonal{T}) = LowerTriangular(A)
 
 function similar{T}(D::Diagonal, ::Type{T}, d::(Int,Int))
     d[1] == d[2] || throw(ArgumentError("Diagonal matrix must be square"))
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 63fe761cc4fa2..9fc00ab461dae 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -274,7 +274,7 @@ function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T})
     end
     A
 end
-A_mul_Bc(A::Triangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B)
+A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B)
 function A_mul_Bc{TA,TB}(A::AbstractArray{TA}, 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(AbstractMatrix{TAB}, A)) :
@@ -284,8 +284,8 @@ function A_mul_Bc{TA,TB}(A::AbstractArray{TA}, B::Union(QRCompactWYQ{TB},QRPacke
               convert(AbstractMatrix{TAB}, B))
 end
 
-A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, b::StridedVector{T}) = (A_ldiv_B!(Triangular(A[:R], :U), sub(Ac_mul_B!(A[:Q], b), 1:size(A, 2))); b)
-A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, B::StridedMatrix{T}) = (A_ldiv_B!(Triangular(A[:R], :U), sub(Ac_mul_B!(A[:Q], B), 1:size(A, 2), 1:size(B, 2))); B)
+A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, b::StridedVector{T}) = (A_ldiv_B!(UpperTriangular(A[:R]), sub(Ac_mul_B!(A[:Q], b), 1:size(A, 2))); b)
+A_ldiv_B!{T<:BlasFloat}(A::QRCompactWY{T}, B::StridedMatrix{T}) = (A_ldiv_B!(UpperTriangular(A[:R]), sub(Ac_mul_B!(A[:Q], B), 1:size(A, 2), 1:size(B, 2))); B)
 
 # Julia implementation similarly to xgelsy
 function A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
@@ -313,7 +313,7 @@ function A_ldiv_B!{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Re
         # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end
     end
     C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:])
-    A_ldiv_B!(Triangular(C[1:rnk,1:rnk],:U),sub(Ac_mul_B!(getq(A),sub(B, 1:mA, 1:nrhs)),1:rnk,1:nrhs))
+    A_ldiv_B!(UpperTriangular(C[1:rnk,1:rnk]),sub(Ac_mul_B!(getq(A),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))
     B[1:nA,:] = sub(B, 1:nA, :)[invperm(A[:p]::Vector{BlasInt}),:]
diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index 419f947cc2dd6..861d101afe83b 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -1,5 +1,9 @@
 ## linalg.jl: Some generic Linear Algebra definitions
 
+# Fall back arithmetic
++(A::AbstractMatrix, B::AbstractMatrix) = full(A) + full(B)
+-(A::AbstractMatrix, B::AbstractMatrix) = full(A) - full(B)
+
 scale(X::AbstractArray, s::Number) = scale!(copy(X), s)
 scale(s::Number, X::AbstractArray) = scale!(copy(X), s)
 
@@ -424,7 +428,7 @@ function elementaryRightTrapezoid!(A::AbstractMatrix, row::Integer)
 end
 
 function det(A::AbstractMatrix)
-    (istriu(A) || istril(A)) && return det(Triangular(A, :U, false))
+    (istriu(A) || istril(A)) && return det(UpperTriangular(A))
     return det(lufact(A))
 end
 det(x::Number) = x
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index 38450cb86b3b7..4fd661735eabc 100644
--- a/base/linalg/lu.jl
+++ b/base/linalg/lu.jl
@@ -114,8 +114,8 @@ function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol)
 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)),:]))
+A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, b::StridedVector) = A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), b[ipiv2perm(A.ipiv, length(b))]))
+A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, B::StridedMatrix) = A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), 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
diff --git a/base/linalg/sparse.jl b/base/linalg/sparse.jl
index 14c4c3081c923..d507bbb8e6a53 100644
--- a/base/linalg/sparse.jl
+++ b/base/linalg/sparse.jl
@@ -126,7 +126,7 @@ end
 *{TvA,TiA}(X::BitArray{2}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (AbstractMatrix, SparseMatrixCSC), X, A)
 # TODO: Tridiagonal * Sparse should be implemented more efficiently
 *{TX,TvA,TiA}(X::Tridiagonal{TX}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (Tridiagonal, AbstractMatrix), X, A)
-*{TvA,TiA}(X::Triangular, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A
+*{TvA,TiA}(X::TriangularUnion, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A
 function *{TX,TvA,TiA}(X::AbstractMatrix{TX}, A::SparseMatrixCSC{TvA,TiA})
     mX, nX = size(X)
     nX == A.m || throw(DimensionMismatch())
diff --git a/base/linalg/special.jl b/base/linalg/special.jl
index 4113e49446d95..0c0af842fe4c0 100644
--- a/base/linalg/special.jl
+++ b/base/linalg/special.jl
@@ -4,8 +4,12 @@
 convert{T}(::Type{Bidiagonal}, A::Diagonal{T})=Bidiagonal(A.diag, zeros(T, size(A.diag,1)-1), true)
 convert{T}(::Type{SymTridiagonal}, A::Diagonal{T})=SymTridiagonal(A.diag, zeros(T, size(A.diag,1)-1))
 convert{T}(::Type{Tridiagonal}, A::Diagonal{T})=Tridiagonal(zeros(T, size(A.diag,1)-1), A.diag, zeros(T, size(A.diag,1)-1))
-convert(::Type{Triangular}, A::Diagonal) = Triangular(full(A), :L)
-convert(::Type{Triangular}, A::Bidiagonal) = Triangular(full(A), A.isupper ? :U : :L)
+convert(::Type{UpperTriangular}, A::Diagonal) = UpperTriangular(full(A), :L)
+convert(::Type{UnitUpperTriangular}, A::Diagonal) = UnitUpperTriangular(full(A), :L)
+convert(::Type{LowerTriangular}, A::Diagonal) = LowerTriangular(full(A), :L)
+convert(::Type{UnitLowerTriangular}, A::Diagonal) = UnitLowerTriangular(full(A), :L)
+convert(::Type{LowerTriangular}, A::Bidiagonal) = !A.isupper ? LowerTriangular(full(A)) : throw(ArgumentError("Bidiagonal matrix must have lower off diagonal to be converted to LowerTriangular"))
+convert(::Type{UpperTriangular}, A::Bidiagonal) = A.isupper ? UpperTriangular(full(A)) : throw(ArgumentError("Bidiagonal matrix must have upper off diagonal to be converted to UpperTriangular"))
 convert(::Type{Matrix}, D::Diagonal) = diagm(D.diag)
 
 function convert(::Type{Diagonal}, A::Union(Bidiagonal, SymTridiagonal))
@@ -42,12 +46,12 @@ function convert(::Type{SymTridiagonal}, A::Tridiagonal)
     SymTridiagonal(A.d, A.dl)
 end
 
-function convert(::Type{Diagonal}, A::Triangular)
+function convert(::Type{Diagonal}, A::TriangularUnion)
     full(A) == diagm(diag(A)) || throw(ArgumentError("Matrix cannot be represented as Diagonal"))
     Diagonal(diag(A))
 end
 
-function convert(::Type{Bidiagonal}, A::Triangular)
+function convert(::Type{Bidiagonal}, A::TriangularUnion)
     fA = full(A)
     if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1)
         return Bidiagonal(diag(A), diag(fA,1), true)
@@ -58,9 +62,9 @@ function convert(::Type{Bidiagonal}, A::Triangular)
     end
 end
 
-convert(::Type{SymTridiagonal}, A::Triangular) = convert(SymTridiagonal, convert(Tridiagonal, A))
+convert(::Type{SymTridiagonal}, A::TriangularUnion) = convert(SymTridiagonal, convert(Tridiagonal, A))
 
-function convert(::Type{Tridiagonal}, A::Triangular)
+function convert(::Type{Tridiagonal}, A::TriangularUnion)
     fA = full(A)
     if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) + diagm(diag(fA, -1), -1)
         return Tridiagonal(diag(fA, -1), diag(A), diag(fA,1))
@@ -82,7 +86,7 @@ macro commutative(myexpr)
 end
 
 for op in (:+, :-)
-    SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Triangular, :Matrix]
+    SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Matrix]
     for (idx, matrixtype1) in enumerate(SpecialMatrices) #matrixtype1 is the sparser matrix type
         for matrixtype2 in SpecialMatrices[idx+1:end] #matrixtype2 is the denser matrix type
             @eval begin #TODO quite a few of these conversions are NOT defined...
@@ -93,7 +97,7 @@ for op in (:+, :-)
     end
 
     for  matrixtype1 in (:SymTridiagonal,)                      #matrixtype1 is the sparser matrix type
-        for matrixtype2 in (:Tridiagonal, :Triangular, :Matrix) #matrixtype2 is the denser matrix type
+        for matrixtype2 in (:Tridiagonal, :Matrix) #matrixtype2 is the denser matrix type
             @eval begin
                 ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                 ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
@@ -111,7 +115,7 @@ for op in (:+, :-)
     end
 end
 
-A_mul_Bc!(A::Triangular, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B)
-A_mul_Bc!(A::Triangular, B::QRPackedQ) = A_mul_Bc!(full!(A),B)
-A_mul_Bc(A::Triangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B)
+A_mul_Bc!(A::TriangularUnion, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B)
+A_mul_Bc!(A::TriangularUnion, B::QRPackedQ) = A_mul_Bc!(full!(A),B)
+A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B)
 
diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl
index b344573187ab1..b0658b1029e3f 100644
--- a/base/linalg/triangular.jl
+++ b/base/linalg/triangular.jl
@@ -1,141 +1,226 @@
 ## Triangular
-immutable Triangular{T,S<:AbstractMatrix,UpLo,IsUnit} <: AbstractMatrix{T}
-    data::S
-end
-function Triangular{T}(A::AbstractMatrix{T}, uplo::Symbol, isunit::Bool=false)
-    chksquare(A)
-    Triangular{T,typeof(A),uplo,isunit}(A)
+
+# First loop through all methods that don't need special care for upper/lower and unit diagonal
+for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular)
+    @eval begin
+        immutable $t{T,S<:AbstractMatrix} <: AbstractMatrix{T}
+            data::S
+        end
+        function $t(A::AbstractMatrix)
+            Base.LinAlg.chksquare(A)
+            return $t{eltype(A), typeof(A)}(A)
+        end
+
+        size(A::$t, args...) = size(A.data, args...)
+
+        convert{T,S}(::Type{$t{T}}, A::$t{T,S}) = A
+        convert{Tnew,Told,S}(::Type{$t{Tnew}}, A::$t{Told,S}) = (Anew = convert(AbstractMatrix{Tnew}, A.data); $t(Anew))
+        convert{Tnew,Told,S}(::Type{AbstractMatrix{Tnew}}, A::$t{Told,S}) = convert($t{Tnew}, A)
+        convert{T,S}(::Type{Matrix}, A::$t{T,S}) = convert(Matrix{T}, A)
+
+        function similar{T,S,Tnew}(A::$t{T,S}, ::Type{Tnew}, dims::Dims)
+            dims[1] == dims[2] || throw(ArgumentError("a Triangular matrix must be square"))
+            length(dims) == 2 || throw(ArgumentError("a Triangular matrix must have two dimensions"))
+            B = similar(A.data, Tnew, dims)
+            return $t(B)
+        end
+
+        copy{T,S}(A::$t{T,S}) = $t{T,S}(copy(A.data))
+
+        big(A::$t) = $t(big(A.data))
+
+        real{T<:Real}(A::$t{T}) = A
+        real{T<:Complex}(A::$t{T}) = (B = real(A.data); $t(B))
+
+    end
 end
 
-size(A::Triangular, args...) = size(A.data, args...)
+# some cases can be handled with a union
+typealias TriangularUnion{T,S} Union(UpperTriangular{T,S}, UnitUpperTriangular{T,S}, LowerTriangular{T,S}, UnitLowerTriangular{T,S})
+
+full{T,S}(A::TriangularUnion{T,S}) = convert(Matrix, A)
 
-convert{T,S,UpLo,IsUnit}(::Type{Triangular{T}}, A::Triangular{T,S,UpLo,IsUnit}) = A
-convert{Tnew,Told,S,UpLo,IsUnit}(::Type{Triangular{Tnew}}, A::Triangular{Told,S,UpLo,IsUnit}) = (Anew = convert(AbstractMatrix{Tnew}, A.data); Triangular(Anew, UpLo, IsUnit))
-convert{Tnew,Told,S,UpLo,IsUnit}(::Type{AbstractMatrix{Tnew}}, A::Triangular{Told,S,UpLo,IsUnit}) = convert(Triangular{Tnew}, A)
-function convert{Tret,T,S,UpLo,IsUnit}(::Type{Matrix{Tret}}, A::Triangular{T,S,UpLo,IsUnit})
+fill!(A::TriangularUnion, x) = (fill!(A.data, x); A)
+
+# then handle all methods that requires specific handling of upper/lower and unit diagonal
+
+function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::LowerTriangular{T,S})
     B = Array(Tret, size(A, 1), size(A, 1))
     copy!(B, A.data)
-    (UpLo == :L ? tril! : triu!)(B)
-    if IsUnit
-        for i = 1:size(B,1)
-            B[i,i] = 1
-        end
+    tril!(B)
+    B
+end
+function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UnitLowerTriangular{T,S})
+    B = Array(Tret, size(A, 1), size(A, 1))
+    copy!(B, A.data)
+    tril!(B)
+    for i = 1:size(B,1)
+        B[i,i] = 1
+    end
+    B
+end
+function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UpperTriangular{T,S})
+    B = Array(Tret, size(A, 1), size(A, 1))
+    copy!(B, A.data)
+    triu!(B)
+    B
+end
+function convert{Tret,T,S}(::Type{Matrix{Tret}}, A::UnitUpperTriangular{T,S})
+    B = Array(Tret, size(A, 1), size(A, 1))
+    copy!(B, A.data)
+    triu!(B)
+    for i = 1:size(B,1)
+        B[i,i] = 1
     end
     B
 end
-convert{T,S,UpLo,IsUnit}(::Type{Matrix}, A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix{T}, A)
 
-function full!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
+function full!{T,S}(A::LowerTriangular{T,S})
     B = A.data
-    (UpLo == :L ? tril! : triu!)(B)
-    if IsUnit
-        for i = 1:size(A,1)
-            B[i,i] = 1
-        end
+    tril!(B)
+    B
+end
+function full!{T,S}(A::UnitLowerTriangular{T,S})
+    B = A.data
+    tril!(B)
+    for i = 1:size(A,1)
+        B[i,i] = 1
     end
     B
 end
-full{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = convert(Matrix, A)
-
-fill!(A::Triangular, x) = (fill!(A.data, x); A)
-
-function similar{T,S,UpLo,IsUnit,Tnew}(A::Triangular{T,S,UpLo,IsUnit}, ::Type{Tnew}, dims::Dims)
-    dims[1] == dims[2] || throw(ArgumentError("a Triangular matrix must be square"))
-    length(dims) == 2 || throw(ArgumentError("a Triangular matrix must have two dimensions"))
-    A = similar(A.data, Tnew, dims)
-    return Triangular{Tnew, typeof(A), UpLo, IsUnit}(A)
+function full!{T,S}(A::UpperTriangular{T,S})
+    B = A.data
+    triu!(B)
+    B
 end
-
-copy{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T,S,UpLo,IsUnit}(copy(A.data))
-
-getindex(A::Triangular, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1])
-getindex{T,S}(A::Triangular{T,S,:L,true}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j]))
-getindex{T,S}(A::Triangular{T,S,:L,false}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j])
-getindex{T,S}(A::Triangular{T,S,:U,true}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j]))
-getindex{T,S}(A::Triangular{T,S,:U,false}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(A.data[i,j])
-
-istril{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :L
-istriu{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :U
-
-transpose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(transpose(A.data))
-ctranspose{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(ctranspose(A.data))
-transpose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U'))
-ctranspose!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = Triangular{T, S, UpLo == :U ? :L : :U, IsUnit}(copytri!(A.data, UpLo == :L ? 'L' : 'U', true))
-diag{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = IsUnit ? ones(T, size(A,1)) : diag(A.data)
-function big{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
-    M = big(A.data)
-    Triangular{eltype(M),typeof(M),UpLo,IsUnit}(M)
+function full!{T,S}(A::UnitUpperTriangular{T,S})
+    B = A.data
+    triu!(B)
+    for i = 1:size(A,1)
+        B[i,i] = 1
+    end
+    B
 end
 
-real{T<:Real}(A::Triangular{T}) = A
-real{T<:Complex,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = (B = real(A.data); Triangular{eltype(B), typeof(B), UpLo, IsUnit}(B))
+getindex(A::TriangularUnion, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1])
+getindex{T,S}(A::UnitLowerTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j]))
+getindex{T,S}(A::LowerTriangular{T,S}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j])
+getindex{T,S}(A::UnitUpperTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j]))
+getindex{T,S}(A::UpperTriangular{T,S}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(A.data[i,j])
+
+istril{T,S}(A::LowerTriangular{T,S}) = true
+istril{T,S}(A::UnitLowerTriangular{T,S}) = true
+istril{T,S}(A::UpperTriangular{T,S}) = false
+istril{T,S}(A::UnitUpperTriangular{T,S}) = false
+istriu{T,S}(A::LowerTriangular{T,S}) = false
+istriu{T,S}(A::UnitLowerTriangular{T,S}) = false
+istriu{T,S}(A::UpperTriangular{T,S}) = true
+istriu{T,S}(A::UnitUpperTriangular{T,S}) = true
+
+transpose{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(transpose(A.data))
+transpose{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(transpose(A.data))
+transpose{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(transpose(A.data))
+transpose{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(transpose(A.data))
+ctranspose{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(ctranspose(A.data))
+ctranspose{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(ctranspose(A.data))
+ctranspose{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(ctranspose(A.data))
+ctranspose{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(ctranspose(A.data))
+
+transpose!{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(copytri!(A.data, 'L'))
+transpose!{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(copytri!(A.data, 'L'))
+transpose!{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(copytri!(A.data, 'U'))
+transpose!{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(copytri!(A.data, 'U'))
+ctranspose!{T,S}(A::LowerTriangular{T,S}) = UpperTriangular{T, S}(copytri!(A.data, 'L' , true))
+ctranspose!{T,S}(A::UnitLowerTriangular{T,S}) = UnitUpperTriangular{T, S}(copytri!(A.data, 'L' , true))
+ctranspose!{T,S}(A::UpperTriangular{T,S}) = LowerTriangular{T, S}(copytri!(A.data, 'U' , true))
+ctranspose!{T,S}(A::UnitUpperTriangular{T,S}) = UnitLowerTriangular{T, S}(copytri!(A.data, 'U' , true))
+
+diag(A::LowerTriangular) = diag(A.data)
+diag(A::UnitLowerTriangular) = ones(eltype(A), size(A,1))
+diag(A::UpperTriangular) = diag(A.data)
+diag(A::UnitUpperTriangular) = ones(eltype(A), size(A,1))
 
 # Unary operations
--{T, S, UpLo}(A::Triangular{T, S, UpLo, false}) = Triangular{T, S, UpLo, false}(-A.data)
-function -{T, S, UpLo}(A::Triangular{T, S, UpLo, true})
+-(A::LowerTriangular) = LowerTriangular(-A.data)
+-(A::UpperTriangular) = UpperTriangular(-A.data)
+function -(A::UnitLowerTriangular)
     Anew = -A.data
     for i = 1:size(A, 1)
         Anew[i, i] = -1
     end
-    Triangular{T, S, UpLo, false}(Anew)
+    LowerTriangular(Anew)
+end
+function -(A::UnitUpperTriangular)
+    Anew = -A.data
+    for i = 1:size(A, 1)
+        Anew[i, i] = -1
+    end
+    UpperTriangular(Anew)
 end
 
 # Binary operations
-+{TA, TB, SA, SB, uplo}(A::Triangular{TA, SA, uplo, false}, B::Triangular{TB, SB, uplo, false}) = Triangular(A.data + B.data, uplo)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, false}, B::Triangular{TB, SB, :U, true}) = Triangular(A.data + triu(B.data, 1) + I, :U)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, false}, B::Triangular{TB, SB, :L, true}) = Triangular(A.data + tril(B.data, -1) + I, :L)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, false}) = Triangular(triu(A.data, 1) + B.data + I, :U)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, false}) = Triangular(tril(A.data, -1) + B.data + I, :L)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, true}) = Triangular(triu(A.data, 1) + triu(B.data, 1) + 2I, :U)
-+{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, true}) = Triangular(tril(A.data, -1) + tril(B.data, -1) + 2I, :L)
-+{TA, SA, TB, SB, uplo1, uplo2, IsUnit1, IsUnit2}(A::Triangular{TA, SA, uplo1, IsUnit1}, B::Triangular{TB, SB, uplo2, IsUnit2}) = full(A) + full(B)
--{TA, SA, TB, SB, uplo}(A::Triangular{TA, SA, uplo, false}, B::Triangular{TB, SB, uplo, false}) = Triangular(A.data - B.data, uplo)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, false}, B::Triangular{TB, SB, :U, true}) = Triangular(A.data - triu(B.data, 1) - I, :U)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, false}, B::Triangular{TB, SB, :L, true}) = Triangular(A.data - tril(B.data, -1) - I, :L)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, false}) = Triangular(triu(A.data, 1) - B.data + I, :U)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, false}) = Triangular(tril(A.data, -1) - B.data + I, :L)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :U, true}, B::Triangular{TB, SB, :U, true}) = Triangular(triu(A.data, 1) - triu(B.data, 1), :U)
--{TA, SA, TB, SB}(A::Triangular{TA, SA, :L, true}, B::Triangular{TB, SB, :L, true}) = Triangular(tril(A.data, -1) - tril(B.data, -1), :L)
--{TA, SA, TB, SB, uplo1, uplo2, IsUnit1, IsUnit2}(A::Triangular{TA, SA, uplo1, IsUnit1}, B::Triangular{TB, SB, uplo2, IsUnit2}) = full(A) - full(B)
++(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data + B.data)
++(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data + B.data)
++(A::UpperTriangular, B::UnitUpperTriangular) = UpperTriangular(A.data + triu(B.data, 1) + I)
++(A::LowerTriangular, B::UnitLowerTriangular) = LowerTriangular(A.data + tril(B.data, -1) + I)
++(A::UnitUpperTriangular, B::UpperTriangular) = UpperTriangular(triu(A.data, 1) + B.data + I)
++(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) + B.data + I)
++(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) + triu(B.data, 1) + 2I)
++(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) + tril(B.data, -1) + 2I)
+
+-(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data - B.data)
+-(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data - B.data)
+-(A::UpperTriangular, B::UnitUpperTriangular) = UpperTriangular(A.data - triu(B.data, 1) - I)
+-(A::LowerTriangular, B::UnitLowerTriangular) = LowerTriangular(A.data - tril(B.data, -1) - I)
+-(A::UnitUpperTriangular, B::UpperTriangular) = UpperTriangular(triu(A.data, 1) - B.data + I)
+-(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) - B.data + I)
+-(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) - triu(B.data, 1))
+-(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) - tril(B.data, -1))
 
 ######################
 # BlasFloat routines #
 ######################
 
-for (uplos, uploc) in ((:(:U), 'U'), (:(:L), 'L'))
-    for (isunitb, isunitc) in ((true, 'U'), (false, 'N'))
-        @eval begin
+A_mul_B!(A::Tridiagonal, B::TriangularUnion) = A*full!(B)
+A_mul_B!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B))
+A_mul_Bc!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B))
+
+for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
+                            (:UnitLowerTriangular, 'L', 'U'),
+                            (:UpperTriangular, 'U', 'N'),
+                            (:UnitUpperTriangular, 'U', 'U'))
+    @eval begin
         # Vector multiplication
-        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, b::StridedVector{T}) = BLAS.trmv!($uploc, 'N', $isunitc, A.data, b)
+        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, b::StridedVector{T}) = BLAS.trmv!($uploc, 'N', $isunitc, A.data, b)
 
         # Matrix multiplication
-        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B)
-        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
+        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'N', $isunitc, one(T), A.data, B)
+        A_mul_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
 
-        Ac_mul_B!{T<:BlasComplex,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B)
-        Ac_mul_B!{T<:BlasReal,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)
+        Ac_mul_B!{T<:BlasComplex,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'C', $isunitc, one(T), A.data, B)
+        Ac_mul_B!{T<:BlasReal,S<:StridedMatrix}(A::$t{T,S}, B::StridedMatrix{T}) = BLAS.trmm!('L', $uploc, 'T', $isunitc, one(T), A.data, B)
 
-        A_mul_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)
-        A_mul_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)
+        A_mul_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)
+        A_mul_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trmm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)
 
         # Left division
-        A_ldiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'N', $isunitc, A.data, B)
-        Ac_ldiv_B!{T<:BlasReal,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'T', $isunitc, A.data, B)
-        Ac_ldiv_B!{T<:BlasComplex,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'C', $isunitc, A.data, B)
+        A_ldiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'N', $isunitc, A.data, B)
+        Ac_ldiv_B!{T<:BlasReal,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'T', $isunitc, A.data, B)
+        Ac_ldiv_B!{T<:BlasComplex,S<:StridedMatrix}(A::$t{T,S}, B::StridedVecOrMat{T}) = LAPACK.trtrs!($uploc, 'C', $isunitc, A.data, B)
 
         # Right division
-        A_rdiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedVecOrMat{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
-        A_rdiv_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)
-        A_rdiv_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::Triangular{T,S,$uplos,$isunitb}) = BLAS.trsm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)
+        A_rdiv_B!{T<:BlasFloat,S<:StridedMatrix}(A::StridedVecOrMat{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'N', $isunitc, one(T), B.data, A)
+        A_rdiv_Bc!{T<:BlasReal,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'T', $isunitc, one(T), B.data, A)
+        A_rdiv_Bc!{T<:BlasComplex,S<:StridedMatrix}(A::StridedMatrix{T}, B::$t{T,S}) = BLAS.trsm!('R', $uploc, 'C', $isunitc, one(T), B.data, A)
 
         # Matrix inverse
-        inv{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}) = Triangular{T,S,$uplos,$isunitb}(LAPACK.trtri!($uploc, $isunitc, copy(A.data)))
+        inv{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}) = $t{T,S}(LAPACK.trtri!($uploc, $isunitc, copy(A.data)))
 
         # Error bounds for triangular solve
-        errorbounds{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,$uplos,$isunitb}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = LAPACK.trrfs!($uploc, 'N', $isunitc, A.data, B, X)
+        errorbounds{T<:BlasFloat,S<:StridedMatrix}(A::$t{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = LAPACK.trrfs!($uploc, 'N', $isunitc, A.data, B, X)
 
         # Condition numbers
-        function cond{T<:BlasFloat,S}(A::Triangular{T,S,$uplos,$isunitb}, p::Real=2)
+        function cond{T<:BlasFloat,S}(A::$t{T,S}, p::Real=2)
             chksquare(A)
             if p==1
                 return inv(LAPACK.trcon!('O', $uploc, $isunitc, A.data))
@@ -146,68 +231,115 @@ for (uplos, uploc) in ((:(:U), 'U'), (:(:L), 'L'))
             end
         end
     end
-    end
 end
 
-errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::Triangular{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.")
-function errorbounds{TA<:Number,S<:StridedMatrix,UpLo,IsUnit,TX<:Number,TB<:Number}(A::Triangular{TA,S,UpLo,IsUnit}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB})
+errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::TriangularUnion{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.")
+function errorbounds{TA<:Number,S<:StridedMatrix,TX<:Number,TB<:Number}(A::TriangularUnion{TA,S}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB})
     TAXB = promote_type(TA, TB, TX, Float32)
     errorbounds(convert(AbstractMatrix{TAXB}, A), convert(AbstractArray{TAXB}, X), convert(AbstractArray{TAXB}, B))
 end
 
 # Eigensystems
 ## Notice that trecv works for quasi-triangular matrices and therefore the lower sub diagonal must be zeroed before calling the subroutine
-eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:U,false}) = LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data))
-eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:U,true}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data)))
-eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:L,false}) = LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)')
-eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::Triangular{T,S,:L,true}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)'))
+eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UpperTriangular{T,S}) = LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data))
+eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UnitUpperTriangular{T,S}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('R', 'A', BlasInt[], triu!(A.data)))
+eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::LowerTriangular{T,S}) = LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)')
+eigvecs{T<:BlasFloat,S<:StridedMatrix}(A::UnitLowerTriangular{T,S}) = (for i = 1:size(A, 1); A.data[i,i] = 1;end;LAPACK.trevc!('L', 'A', BlasInt[], tril!(A.data)'))
 
 ####################
 # Generic routines #
 ####################
 
-(*){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data*x, UpLo)
-function (*){T,S,UpLo}(A::Triangular{T,S,UpLo,true}, x::Number)
+(*)(A::UpperTriangular, x::Number) = UpperTriangular(A.data*x)
+(*)(A::LowerTriangular, x::Number) = LowerTriangular(A.data*x)
+function (*)(A::UnitUpperTriangular, x::Number)
+    B = A.data*x
+    for i = 1:size(A, 1)
+        B[i,i] = x
+    end
+    UpperTriangular(B)
+end
+function (*)(A::UnitLowerTriangular, x::Number)
     B = A.data*x
     for i = 1:size(A, 1)
         B[i,i] = x
     end
-    Triangular(B, UpLo)
+    LowerTriangular(B)
+end
+(*)(x::Number, A::UpperTriangular) = UpperTriangular(x*A.data)
+(*)(x::Number, A::LowerTriangular) = LowerTriangular(x*A.data)
+function (*)(x::Number, A::UnitUpperTriangular)
+    B = x*A.data
+    for i = 1:size(A, 1)
+        B[i,i] = x
+    end
+    UpperTriangular(B)
 end
-(*){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x*A.data, UpLo)
-function (*){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,true})
+function (*)(x::Number, A::UnitLowerTriangular)
     B = x*A.data
     for i = 1:size(A, 1)
         B[i,i] = x
     end
-    Triangular(B, UpLo)
+    LowerTriangular(B)
 end
-(/){T,S,UpLo}(A::Triangular{T,S,UpLo,false}, x::Number) = Triangular(A.data/x, UpLo)
-function (/){T,S,UpLo}(A::Triangular{T,S,UpLo,true}, x::Number)
+(/)(A::UpperTriangular, x::Number) = UpperTriangular(A.data/x)
+(/)(A::LowerTriangular, x::Number) = LowerTriangular(A.data/x)
+function (/)(A::UnitUpperTriangular, x::Number)
     B = A.data*x
     invx = inv(x)
     for i = 1:size(A, 1)
         B[i,i] = x
     end
-    Triangular(B, UpLo)
+    UpperTriangular(B)
 end
-(\){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,false}) = Triangular(x\A.data, UpLo)
-function (\){T,S,UpLo}(x::Number, A::Triangular{T,S,UpLo,true})
+function (/)(A::UnitLowerTriangular, x::Number)
+    B = A.data*x
+    invx = inv(x)
+    for i = 1:size(A, 1)
+        B[i,i] = x
+    end
+    LowerTriangular(B)
+end
+(\)(x::Number, A::UpperTriangular) = UpperTriangular(x\A.data)
+(\)(x::Number, A::LowerTriangular) = LowerTriangular(x\A.data)
+function (\)(x::Number, A::UnitUpperTriangular)
     B = x\A.data
     invx = inv(x)
     for i = 1:size(A, 1)
         B[i,i] = invx
     end
-    Triangular(B, UpLo)
+    UpperTriangular(B)
+end
+function (\)(x::Number, A::UnitLowerTriangular)
+    B = x\A.data
+    invx = inv(x)
+    for i = 1:size(A, 1)
+        B[i,i] = invx
+    end
+    LowerTriangular(B)
 end
 
 ## Generic triangular multiplication
-function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::StridedVecOrMat)
+function A_mul_B!(A::UpperTriangular, B::StridedVecOrMat)
+    m, n = size(B, 1), size(B, 2)
+    m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for j = 1:n
+        for i = 1:m
+            Bij = A.data[i,i]*B[i,j]
+            for k = i + 1:m
+                Bij += A.data[i,k]*B[k,j]
+            end
+            B[i,j] = Bij
+        end
+    end
+    B
+end
+function A_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for j = 1:n
         for i = 1:m
-            Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j])
+            Bij = B[i,j]
             for k = i + 1:m
                 Bij += A.data[i,k]*B[k,j]
             end
@@ -217,12 +349,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::S
     B
 end
 
-function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::StridedVecOrMat)
+function A_mul_B!(A::LowerTriangular, B::StridedVecOrMat)
+    m, n = size(B, 1), size(B, 2)
+    m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for j = 1:n
+        for i = m:-1:1
+            Bij = A.data[i,i]*B[i,j]
+            for k = 1:i - 1
+                Bij += A.data[i,k]*B[k,j]
+            end
+            B[i,j] = Bij
+        end
+    end
+    B
+end
+function A_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for j = 1:n
         for i = m:-1:1
-            Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j])
+            Bij = B[i,j]
             for k = 1:i - 1
                 Bij += A.data[i,k]*B[k,j]
             end
@@ -232,12 +378,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::S
     B
 end
 
-function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::StridedVecOrMat)
+function Ac_mul_B!(A::UpperTriangular, B::StridedVecOrMat)
+    m, n = size(B, 1), size(B, 2)
+    m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for j = 1:n
+        for i = m:-1:1
+            Bij = A.data[i,i]*B[i,j]
+            for k = 1:i - 1
+                Bij += A.data[k,i]'B[k,j]
+            end
+            B[i,j] = Bij
+        end
+    end
+    B
+end
+function Ac_mul_B!(A::UnitUpperTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for j = 1:n
         for i = m:-1:1
-            Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j])
+            Bij = B[i,j]
             for k = 1:i - 1
                 Bij += A.data[k,i]'B[k,j]
             end
@@ -247,12 +407,26 @@ function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:U,IsUnit}, B::
     B
 end
 
-function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::StridedVecOrMat)
+function Ac_mul_B!(A::LowerTriangular, B::StridedVecOrMat)
+    m, n = size(B, 1), size(B, 2)
+    m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for j = 1:n
+        for i = 1:m
+            Bij = A.data[i,i]*B[i,j]
+            for k = i + 1:m
+                Bij += A.data[k,i]'B[k,j]
+            end
+            B[i,j] = Bij
+        end
+    end
+    B
+end
+function Ac_mul_B!(A::UnitLowerTriangular, B::StridedVecOrMat)
     m, n = size(B, 1), size(B, 2)
     m == size(A, 1) || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for j = 1:n
         for i = 1:m
-            Bij = ifelse(IsUnit, B[i,j], A.data[i,i]*B[i,j])
+            Bij = B[i,j]
             for k = i + 1:m
                 Bij += A.data[k,i]'B[k,j]
             end
@@ -262,12 +436,26 @@ function Ac_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::Triangular{T,S,:L,IsUnit}, B::
     B
 end
 
-function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
+function A_mul_B!(A::StridedMatrix, B::UpperTriangular)
+    m, n = size(A)
+    size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = n:-1:1
+            Aij = A[i,j]*B[j,j]
+            for k = 1:j - 1
+                Aij += A[i,k]*B.data[k,j]
+            end
+            A[i,j] = Aij
+        end
+    end
+    A
+end
+function A_mul_B!(A::StridedMatrix, B::UnitUpperTriangular)
     m, n = size(A)
     size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
         for j = n:-1:1
-            Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j])
+            Aij = A[i,j]
             for k = 1:j - 1
                 Aij += A[i,k]*B.data[k,j]
             end
@@ -277,12 +465,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,
     A
 end
 
-function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
+function A_mul_B!(A::StridedMatrix, B::LowerTriangular)
+    m, n = size(A)
+    size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = 1:n
+            Aij = A[i,j]*B[j,j]
+            for k = j + 1:n
+                Aij += A[i,k]*B.data[k,j]
+            end
+            A[i,j] = Aij
+        end
+    end
+    A
+end
+function A_mul_B!(A::StridedMatrix, B::UnitLowerTriangular)
     m, n = size(A)
     size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
         for j = 1:n
-            Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j])
+            Aij = A[i,j]
             for k = j + 1:n
                 Aij += A[i,k]*B.data[k,j]
             end
@@ -292,12 +494,26 @@ function A_mul_B!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,
     A
 end
 
-function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
+function A_mul_Bc!(A::StridedMatrix, B::UpperTriangular)
     m, n = size(A)
     size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
         for j = 1:n
-            Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j])
+            Aij = A[i,j]*B[j,j]
+            for k = j + 1:n
+                Aij += A[i,k]*B.data[j,k]'
+            end
+            A[i,j] = Aij
+        end
+    end
+    A
+end
+function A_mul_Bc!(A::StridedMatrix, B::UnitUpperTriangular)
+    m, n = size(A)
+    size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = 1:n
+            Aij = A[i,j]
             for k = j + 1:n
                 Aij += A[i,k]*B.data[j,k]'
             end
@@ -307,12 +523,26 @@ function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T
     A
 end
 
-function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
+function A_mul_Bc!(A::StridedMatrix, B::LowerTriangular)
     m, n = size(A)
     size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
         for j = n:-1:1
-            Aij = ifelse(IsUnit, A[i,j], A[i,j]*B[j,j])
+            Aij = A[i,j]*B[j,j]
+            for k = 1:j - 1
+                Aij += A[i,k]*B.data[j,k]'
+            end
+            A[i,j] = Aij
+        end
+    end
+    A
+end
+function A_mul_Bc!(A::StridedMatrix, B::UnitLowerTriangular)
+    m, n = size(A)
+    size(B, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = n:-1:1
+            Aij = A[i,j]
             for k = 1:j - 1
                 Aij += A[i,k]*B.data[j,k]'
             end
@@ -321,43 +551,70 @@ function A_mul_Bc!{T,S<:AbstractMatrix,IsUnit}(A::StridedMatrix, B::Triangular{T
     end
     A
 end
-
-A_mul_B!(A::Tridiagonal, B::Triangular) = A*full!(B)
-A_mul_B!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B))
-A_mul_Bc!(C::AbstractVecOrMat, A::Triangular, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B))
 
 #Generic solver using naive substitution
-function naivesub!{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, b::AbstractVector, x::AbstractVector=b)
+function naivesub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector=b)
     N = size(A, 2)
-    N==length(b)==length(x) || throw(DimensionMismatch())
-
-    if UpLo == :L #do forward substitution
-        for j = 1:N
-            x[j] = b[j]
-            for k = 1:j-1
-                x[j] -= A[j,k] * x[k]
-            end
-            if !IsUnit
-                x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
-            end
+    N == length(b) == length(x) || throw(DimensionMismatch())
+    for j = N:-1:1
+        x[j] = b[j]
+        for k = j+1:1:N
+            x[j] -= A[j,k] * x[k]
         end
-    elseif UpLo == :U #do backward substitution
-        for j = N:-1:1
-            x[j] = b[j]
-            for k = j+1:1:N
-                x[j] -= A[j,k] * x[k]
-            end
-            if !IsUnit
-                x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
-            end
+        x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
+    end
+    x
+end
+function naivesub!(A::UnitUpperTriangular, b::AbstractVector, x::AbstractVector=b)
+    N = size(A, 2)
+    N == length(b) == length(x) || throw(DimensionMismatch())
+    for j = N:-1:1
+        x[j] = b[j]
+        for k = j+1:1:N
+            x[j] -= A[j,k] * x[k]
+        end
+    end
+    x
+end
+function naivesub!(A::LowerTriangular, b::AbstractVector, x::AbstractVector=b)
+    N = size(A, 2)
+    N == length(b) == length(x) || throw(DimensionMismatch())
+    for j = 1:N
+        x[j] = b[j]
+        for k = 1:j-1
+            x[j] -= A[j,k] * x[k]
+        end
+        x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
+    end
+    x
+end
+function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector=b)
+    N = size(A, 2)
+    N == length(b) == length(x) || throw(DimensionMismatch())
+    for j = 1:N
+        x[j] = b[j]
+        for k = 1:j-1
+            x[j] -= A[j,k] * x[k]
         end
-    else
-        throw(ArgumentError("Unknown UpLo=$(UpLo)"))
     end
     x
 end
 
-function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
+function A_rdiv_B!(A::StridedMatrix, B::UpperTriangular)
+    m, n = size(A)
+    size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = 1:n
+            Aij = A[i,j]
+            for k = 1:j - 1
+                Aij -= A[i,k]*B.data[k,j]
+            end
+            A[i,j] = Aij/B[j,j]
+        end
+    end
+    A
+end
+function A_rdiv_B!(A::StridedMatrix, B::UnitUpperTriangular)
     m, n = size(A)
     size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
@@ -366,13 +623,13 @@ function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
             for k = 1:j - 1
                 Aij -= A[i,k]*B.data[k,j]
             end
-            A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j])
+            A[i,j] = Aij
         end
     end
     A
 end
 
-function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
+function A_rdiv_B!(A::StridedMatrix, B::LowerTriangular)
     m, n = size(A)
     size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
@@ -381,13 +638,41 @@ function A_rdiv_B!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
             for k = j + 1:n
                 Aij -= A[i,k]*B.data[k,j]
             end
-            A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j])
+            A[i,j] = Aij/B[j,j]
+        end
+    end
+    A
+end
+function A_rdiv_B!(A::StridedMatrix, B::UnitLowerTriangular)
+    m, n = size(A)
+    size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = n:-1:1
+            Aij = A[i,j]
+            for k = j + 1:n
+                Aij -= A[i,k]*B.data[k,j]
+            end
+            A[i,j] = Aij
         end
     end
     A
 end
 
-function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
+function A_rdiv_Bc!(A::StridedMatrix, B::UpperTriangular)
+    m, n = size(A)
+    size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = n:-1:1
+            Aij = A[i,j]
+            for k = j + 1:n
+                Aij -= A[i,k]*B.data[j,k]'
+            end
+            A[i,j] = Aij/B[j,j]
+        end
+    end
+    A
+end
+function A_rdiv_Bc!(A::StridedMatrix, B::UnitUpperTriangular)
     m, n = size(A)
     size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
@@ -396,13 +681,27 @@ function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:U,IsUnit})
             for k = j + 1:n
                 Aij -= A[i,k]*B.data[j,k]'
             end
-            A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j])
+            A[i,j] = Aij
         end
     end
     A
 end
 
-function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
+function A_rdiv_Bc!(A::StridedMatrix, B::LowerTriangular)
+    m, n = size(A)
+    size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
+    for i = 1:m
+        for j = 1:n
+            Aij = A[i,j]
+            for k = 1:j - 1
+                Aij -= A[i,k]*B.data[j,k]'
+            end
+            A[i,j] = Aij/B[j,j]
+        end
+    end
+    A
+end
+function A_rdiv_Bc!(A::StridedMatrix, B::UnitLowerTriangular)
     m, n = size(A)
     size(A, 1) == n || throw(DimensionMismatch("left and right hand side doesn't fit"))
     for i = 1:m
@@ -411,25 +710,30 @@ function A_rdiv_Bc!{T,S,IsUnit}(A::StridedMatrix, B::Triangular{T,S,:L,IsUnit})
             for k = 1:j - 1
                 Aij -= A[i,k]*B.data[j,k]'
             end
-            A[i,j] = ifelse(IsUnit, Aij, Aij/B[j,j])
+            A[i,j] = Aij
         end
     end
     A
 end
 
 # Promotion
-## Promotion methods in matmul doesn't apply to triangular multiplication since it is inplace. Hence we have to make very similar definitions, but without allocation of a result array. For multiplication and unit diagonal division the element type doesn't have to be stable under division whereas that is necessary in the general triangular solve problem.
+## Promotion methods in matmul don't apply to triangular multiplication since it is inplace. Hence we have to make very similar definitions, but without allocation of a result array. For multiplication and unit diagonal division the element type doesn't have to be stable under division whereas that is necessary in the general triangular solve problem.
 
 ## Some Triangular-Triangular cases. We might want to write taylored methods for these cases, but I'm not sure it is worth it.
-*(A::Tridiagonal, B::Triangular) = A_mul_B!(full(A), B)
+for t in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
+    @eval begin
+        *(A::Tridiagonal, B::$t) = A_mul_B!(full(A), B)
+    end
+end
+
 for f in (:*, :Ac_mul_B, :At_mul_B, :A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :\, :Ac_ldiv_B, :At_ldiv_B)
     @eval begin
-        ($f)(A::Triangular, B::Triangular) = ($f)(A, full(B))
+        ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(A, full(B))
     end
 end
 for f in (:A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :/, :A_rdiv_Bc, :A_rdiv_Bt)
     @eval begin
-        ($f)(A::Triangular, B::Triangular) = ($f)(full(A), B)
+        ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(full(A), B)
     end
 end
 
@@ -437,7 +741,7 @@ end
 ### Multiplication with triangle to the left and hence rhs cannot be transposed.
 for (f, g) in ((:*, :A_mul_B!), (:Ac_mul_B, :Ac_mul_B!), (:At_mul_B, :At_mul_B!))
     @eval begin
-        function ($f){TA,TB}(A::Triangular{TA}, B::StridedVecOrMat{TB})
+        function ($f){TA,TB}(A::TriangularUnion{TA}, B::StridedVecOrMat{TB})
             TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
@@ -446,7 +750,15 @@ end
 ### Left division with triangle to the left hence rhs cannot be transposed. No quotients.
 for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
     @eval begin
-        function ($f){TA,TB,S,UpLo}(A::Triangular{TA,S,UpLo,true}, B::StridedVecOrMat{TB})
+        function ($f){TA,TB,S}(A::UnitUpperTriangular{TA,S}, B::StridedVecOrMat{TB})
+            TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
+            ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
+        end
+    end
+end
+for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
+    @eval begin
+        function ($f){TA,TB,S}(A::UnitLowerTriangular{TA,S}, B::StridedVecOrMat{TB})
             TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
@@ -455,7 +767,15 @@ end
 ### Left division with triangle to the left hence rhs cannot be transposed. Quotients.
 for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
     @eval begin
-        function ($f){TA,TB,S,UpLo}(A::Triangular{TA,S,UpLo,false}, B::StridedVecOrMat{TB})
+        function ($f){TA,TB,S}(A::UpperTriangular{TA,S}, B::StridedVecOrMat{TB})
+            TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
+            ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
+        end
+    end
+end
+for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
+    @eval begin
+        function ($f){TA,TB,S}(A::LowerTriangular{TA,S}, B::StridedVecOrMat{TB})
             TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
@@ -464,7 +784,7 @@ end
 ### Multiplication with triangle to the rigth and hence lhs cannot be transposed.
 for (f, g) in ((:*, :A_mul_B!), (:A_mul_Bc, :A_mul_Bc!), (:A_mul_Bt, :A_mul_Bt!))
     @eval begin
-        function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::Triangular{TB})
+        function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::TriangularUnion{TB})
             TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
@@ -473,7 +793,15 @@ end
 ### Right division with triangle to the right hence lhs cannot be transposed. No quotients.
 for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
     @eval begin
-        function ($f){TA,TB,S,UpLo}(A::StridedVecOrMat{TA}, B::Triangular{TB,S,UpLo,true})
+        function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UnitUpperTriangular{TB,S})
+            TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
+            ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
+        end
+    end
+end
+for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
+    @eval begin
+        function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UnitLowerTriangular{TB,S})
             TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
@@ -482,52 +810,76 @@ end
 ### Right division with triangle to the right hence lhs cannot be transposed. Quotients.
 for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
     @eval begin
-        function ($f){TA,TB,S,UpLo}(A::StridedVecOrMat{TA}, B::Triangular{TB,S,UpLo,false})
+        function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::UpperTriangular{TB,S})
+            TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
+            ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
+        end
+    end
+end
+for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
+    @eval begin
+        function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::LowerTriangular{TB,S})
             TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
             ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B))
         end
     end
 end
 
-function sqrtm{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit})
+function sqrtm{T}(A::UpperTriangular{T})
     n = size(A, 1)
     TT = typeof(sqrt(zero(T)))
     R = zeros(TT, n, n)
-    if UpLo == :U
-        for j = 1:n
-            (T<:Complex || A[j,j]>=0) ? (R[j,j] = IsUnit ? one(T) : sqrt(A[j,j])) : throw(SingularException(j))
-            for i = j-1:-1:1
-                r = A[i,j]
-                for k = i+1:j-1
-                    r -= R[i,k]*R[k,j]
-                end
-                r==0 || (R[i,j] = r / (R[i,i] + R[j,j]))
+    for j = 1:n
+        R[j,j] = sqrt(A[j,j])
+        for i = j-1:-1:1
+            r = A[i,j]
+            for k = i+1:j-1
+                r -= R[i,k]*R[k,j]
+            end
+            r==0 || (R[i,j] = r / (R[i,i] + R[j,j]))
+        end
+    end
+    return UpperTriangular(R)
+end
+function sqrtm{T}(A::UnitUpperTriangular{T})
+    n = size(A, 1)
+    TT = typeof(sqrt(zero(T)))
+    R = zeros(TT, n, n)
+    for j = 1:n
+        R[j,j] = one(T)
+        for i = j-1:-1:1
+            r = A[i,j]
+            for k = i+1:j-1
+                r -= R[i,k]*R[k,j]
             end
+            r==0 || (R[i,j] = r / (R[i,i] + R[j,j]))
         end
-        return Triangular(R, :U, IsUnit)
-    else # UpLo == :L #Not the usual case
-        return sqrtm(A.').'
     end
+    return UnitUpperTriangular(R)
 end
+sqrtm(A::LowerTriangular) = sqrtm(A.').'
+sqrtm(A::UnitLowerTriangular) = sqrtm(A.').'
 
 #Generic eigensystems
-eigvals(A::Triangular) = diag(A)
-function eigvecs{T}(A::Triangular{T})
+eigvals(A::TriangularUnion) = diag(A)
+function eigvecs{T}(A::TriangularUnion{T})
     TT = promote_type(T, Float32)
     TT <: BlasFloat && return eigvecs(convert(AbstractMatrix{TT}, A))
     error("type not handled yet. Please submit a pull request.")
 end
-det{T,S,UpLo}(A::Triangular{T,S,UpLo,true}) = one(T)*one(T)
-det{T,S,UpLo}(A::Triangular{T,S,UpLo,false}) = prod(diag(A.data))
+det{T}(A::UnitUpperTriangular{T}) = one(T)*one(T)
+det{T}(A::UnitLowerTriangular{T}) = one(T)*one(T)
+det{T}(A::UpperTriangular{T}) = prod(diag(A.data))
+det{T}(A::LowerTriangular{T}) = prod(diag(A.data))
 
-eigfact(A::Triangular) = Eigen(eigvals(A), eigvecs(A))
+eigfact(A::TriangularUnion) = Eigen(eigvals(A), eigvecs(A))
 
 #Generic singular systems
 for func in (:svd, :svdfact, :svdfact!, :svdvals)
     @eval begin
-        ($func)(A::Triangular) = ($func)(full(A))
+        ($func)(A::TriangularUnion) = ($func)(full(A))
     end
 end
 
-factorize(A::Triangular) = A
+factorize(A::TriangularUnion) = A
 
diff --git a/base/linalg/uniformscaling.jl b/base/linalg/uniformscaling.jl
index a13a2f0d42d51..4bf9627540511 100644
--- a/base/linalg/uniformscaling.jl
+++ b/base/linalg/uniformscaling.jl
@@ -84,7 +84,7 @@ inv(J::UniformScaling) = UniformScaling(inv(J.λ))
 /(J::UniformScaling, x::Number) = UniformScaling(J.λ/x)
 
 \(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ)
-\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), J::UniformScaling) = inv(A)*J.λ
+\{T<:Number}(A::Union(Bidiagonal{T},TriangularUnion{T}), J::UniformScaling) = inv(A)*J.λ
 \(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
 \(A::AbstractMatrix, J::UniformScaling) = inv(A)*J.λ
 
diff --git a/doc/manual/linear-algebra.rst b/doc/manual/linear-algebra.rst
index 9e6dd72f4bec0..733dc96eca888 100644
--- a/doc/manual/linear-algebra.rst
+++ b/doc/manual/linear-algebra.rst
@@ -46,7 +46,8 @@ for them in LAPACK are available.
 
 ======================= ==================================================================================
 :class:`Hermitian`      `Hermitian matrix <http://en.wikipedia.org/wiki/Hermitian_matrix>`_
-:class:`Triangular`     Upper/lower `triangular matrix <http://en.wikipedia.org/wiki/Triangular_matrix>`_
+:class:`UpperTriangular Upper `triangular matrix <http://en.wikipedia.org/wiki/Triangular_matrix>`_
+:class:`LowerTriangular Lower `triangular matrix <http://en.wikipedia.org/wiki/Triangular_matrix>`_
 :class:`Tridiagonal`    `Tridiagonal matrix <http://en.wikipedia.org/wiki/Tridiagonal_matrix>`_
 :class:`SymTridiagonal` Symmetric tridiagonal matrix
 :class:`Bidiagonal`     Upper/lower `bidiagonal matrix <http://en.wikipedia.org/wiki/Bidiagonal_matrix>`_
@@ -64,7 +65,9 @@ Elementary operations
 | :class:`Hermitian`      |       |       |       |   MV  | :func:`inv`,                   |
 |                         |       |       |       |       | :func:`sqrtm`, :func:`expm`    |
 +-------------------------+-------+-------+-------+-------+--------------------------------+
-| :class:`Triangular`     |       |       |  MV   |   MV  | :func:`inv`, :func:`det`       |
+| :class:`UpperTriangular`|       |       |  MV   |   MV  | :func:`inv`, :func:`det`       |
++-------------------------+-------+-------+-------+-------+--------------------------------+
+| :class:`LowerTriangular`|       |       |  MV   |   MV  | :func:`inv`, :func:`det`       |
 +-------------------------+-------+-------+-------+-------+--------------------------------+
 | :class:`SymTridiagonal` |   M   |   M   |  MS   |   MV  | :func:`eigmax`, :func:`eigmin` |
 +-------------------------+-------+-------+-------+-------+--------------------------------+
@@ -96,7 +99,9 @@ Matrix factorizations
 +=========================+========+=============+=================+=================+=============+=================+
 | :class:`Hermitian`      |   HE   |             |       ARI       |                 |             |                 |
 +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+
-| :class:`Triangular`     |   TR   |             |                 |                 |             |                 |
+| :class:`UpperTriangular`|   TR   |      A      |        A        |       A         |             |                 |
++-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+
+| :class:`LowerTriangular`|   TR   |      A      |        A        |       A         |             |                 |
 +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+
 | :class:`SymTridiagonal` |   ST   |      A      |       ARI       |       AV        |             |                 |
 +-------------------------+--------+-------------+-----------------+-----------------+-------------+-----------------+
diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst
index 4f1b83831787a..e118da5312267 100644
--- a/doc/stdlib/linalg.rst
+++ b/doc/stdlib/linalg.rst
@@ -39,7 +39,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: factorize(A)
 
-   Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, Triangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example: ``A=factorize(A); x=A\\b; y=A\\C``.
+   Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, LowerTriangular, UpperTriangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example: ``A=factorize(A); x=A\\b; y=A\\C``.
 
 .. function:: factorize!(A)
 
diff --git a/test/linalg/triangular.jl b/test/linalg/triangular.jl
index bbb8be43c1734..785a6c3dbb16a 100644
--- a/test/linalg/triangular.jl
+++ b/test/linalg/triangular.jl
@@ -1,6 +1,6 @@
 debug = false
 using Base.Test
-using Base.LinAlg: BlasFloat, errorbounds, full!, naivesub!, transpose!
+using Base.LinAlg: BlasFloat, errorbounds, full!, naivesub!, transpose!, UnitUpperTriangular, UnitLowerTriangular
 
 debug && println("Triangular matrices")
 
@@ -8,145 +8,147 @@ n = 9
 srand(123)
 
 debug && println("Test basic type functionality")
-@test_throws DimensionMismatch Triangular(randn(5, 4), :L)
-@test Triangular(randn(3, 3), :L) |> t -> [size(t, i) for i = 1:3] == [size(full(t), i) for i = 1:3]
+@test_throws DimensionMismatch LowerTriangular(randn(5, 4))
+@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(full(t), i) for i = 1:3]
 
 # The following test block tries to call all methods in base/linalg/triangular.jl in order for a combination of input element types. Keep the ordering when adding code.
 for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
     for elty2 in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
-        for uplo1 in (:U, :L)
-            for uplo2 in (:U, :L)
-                for isunit1 in (true, false)
-                    for isunit2 in (true, false)
-                        A1 = Triangular(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo1)), uplo1, isunit1)
-                        A2 = Triangular(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo2)), uplo2, isunit2)
-
-                        for eltyB in (Float32, Float64, Complex64, Complex128)
-                            B = convert(Matrix{eltyB}, elty1 <: Complex ? real(A1)*ones(n, n) : A1*ones(n, n))
-
-                            debug && println("A1: $elty1, $uplo1, $isunit1.A2: $elty2, $uplo2, $isunit2. B: $eltyB")
-
-                            # Convert
-                            @test convert(Triangular{elty1}, A1) == A1
-                            if elty1 <: Real && !(elty2 <: Integer)
-                                @test convert(Triangular{elty2}, A1) == Triangular(convert(Matrix{elty2}, A1.data), uplo1, isunit1)
-                            elseif elty2 <: Real && !(elty1 <: Integer)
-                                @test_throws InexactError convert(Triangular{elty2}, A1) == Triangular(convert(Matrix{elty2}, A1.data), uplo1, isunit1)
-                            end
-                            @test convert(Matrix, A1) == full(A1)
-
-                            # full!
-                            @test full!(copy(A1)) == full(A1)
-
-                            # fill!
-                            @test full!(fill!(copy(A1), 1)) == full(Triangular(ones(size(A1)...), uplo1, isunit1))
-
-                            # similar
-                            @test isa(similar(A1), Triangular{elty1, Matrix{elty1}, uplo1, isunit1})
-
-                            # Linear indexing
-                            for i = 1:length(A1)
-                                @test A1[i] == full(A1)[i]
-                            end
-
-                            # istril/istriu
-                            if uplo1 == :L
-                                @test istril(A1)
-                                @test !istriu(A1)
-                            else
-                                @test istriu(A1)
-                                @test !istril(A1)
-                            end
-
-                            # (c)transpose
-                            @test full(A1') == full(A1)'
-                            @test full(A1.') == full(A1).'
-                            @test transpose!(copy(A1)) == A1.'
-
-                            # diag
-                            @test diag(A1) == diag(full(A1))
-
-                            # real
-                            @test full(real(A1)) == real(full(A1))
-
-                            # Unary operations
-                            @test full(-A1) == -full(A1)
-
-                            # Binary operations
-                            @test full(A1 + A2) == full(A1) + full(A2)
-                            @test full(A1 - A2) == full(A1) - full(A2)
-                            @test A1*0.5 == full(A1*0.5)
-                            @test 0.5*A1 == full(0.5*A1)
-                            @test A1/0.5 == full(A1/0.5)
-                            @test 0.5\A1 == full(0.5\A1)
-
-                            # Triangular-Triangular multiplication and division
-                            elty1 != BigFloat && elty2 != BigFloat && @test_approx_eq full(A1*A2) full(A1)*full(A2)
-                            @test_approx_eq full(A1'A2) full(A1)'full(A2)
-                            @test_approx_eq full(A1*A2') full(A1)*full(A2)'
-                            @test_approx_eq full(A1'A2') full(A1)'full(A2)'
-                            @test_approx_eq full(A1/A2) full(A1)/full(A2)
-
-                            # Triangular-dense Matrix/vector multiplication
-                            @test_approx_eq A1*B[:,1] full(A1)*B[:,1]
-                            @test_approx_eq A1*B full(A1)*B
-                            @test_approx_eq A1'B[:,1] full(A1)'B[:,1]
-                            @test_approx_eq A1'B full(A1)'B
-                            @test_approx_eq A1*B' full(A1)*B'
-                            @test_approx_eq B*A1 B*full(A1)
-                            @test_approx_eq B[:,1]'A1 B[:,1]'full(A1)
-                            @test_approx_eq B'A1 B'full(A1)
-                            @test_approx_eq B*A1' B*full(A1)'
-                            @test_approx_eq B[:,1]'A1' B[:,1]'full(A1)'
-                            @test_approx_eq B'A1' B'full(A1)'
-
-                            # ... and division
-                            @test_approx_eq A1\B[:,1] full(A1)\B[:,1]
-                            @test_approx_eq A1\B full(A1)\B
-                            @test_approx_eq A1'\B[:,1] full(A1)'\B[:,1]
-                            @test_approx_eq A1'\B full(A1)'\B
-                            @test_approx_eq A1\B' full(A1)\B'
-                            @test_approx_eq A1'\B' full(A1)'\B'
-                            @test_approx_eq B/A1 B/full(A1)
-                            @test_approx_eq B/A1' B/full(A1)'
-                            @test_approx_eq B'/A1 B'/full(A1)
-                            @test_approx_eq B'/A1' B'/full(A1)'
-
-                            # inversion
-                            @test_approx_eq inv(A1) inv(lufact(full(A1)))
-
-                            # Error bounds
-                            elty1 != BigFloat && errorbounds(A1, A1\B, B)
-
-                            # Determinant
-                            @test_approx_eq det(A1) det(lufact(full(A1)))
-
-                            # Matri square root
-                            @test_approx_eq sqrtm(A1) |> t->t*t A1
-
-                            # eigenproblems
-                            if elty1 != BigFloat # Not handled yet
-                                vals, vecs = eig(A1)
-                                if !isunit1 && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues.
-                                    @test_approx_eq_eps vecs*diagm(vals)/vecs full(A1) sqrt(eps(float(real(one(vals[1])))))*(norm(A1, Inf)*n)^2
-                                end
-                            end
-
-                            # Condition number tests - can be VERY approximate
-                            if elty1 <:BlasFloat
-                                for p in (1.0, Inf)
-                                    @test_approx_eq_eps cond(A1, p) cond(A1, p) (cond(A1, p) + cond(A1, p))
-                                end
-                            end
-
-                            if elty1 != BigFloat # Not implemented yet
-                                svd(A1)
-                                svdfact(A1)
-                                elty1 <: BlasFloat && svdfact!(copy(A1))
-                                svdvals(A1)
-                            end
+        for (t1, uplo1) in ((UpperTriangular, :U),
+                            (UnitUpperTriangular, :U),
+                            (LowerTriangular, :L),
+                            (UnitLowerTriangular, :L))
+            for (t2, uplo2) in ((UpperTriangular, :U),
+                                (UnitUpperTriangular, :U),
+                                (LowerTriangular, :L),
+                                (UnitLowerTriangular, :L))
+                A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't, uplo1)))
+                A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, uplo2)))
+
+                for eltyB in (Float32, Float64, Complex64, Complex128)
+                    B = convert(Matrix{eltyB}, elty1 <: Complex ? real(A1)*ones(n, n) : A1*ones(n, n))
+
+                    debug && println("elty1: $elty1, A1: $t1, elty2: $elty2, A2: $t2, B: $eltyB")
+
+                    # Convert
+                    @test convert(AbstractMatrix{elty1}, A1) == A1
+                    if elty1 <: Real && !(elty2 <: Integer)
+                        @test convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
+                    elseif elty2 <: Real && !(elty1 <: Integer)
+                        @test_throws InexactError convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
+                    end
+                    @test convert(Matrix, A1) == full(A1)
+
+                    # full!
+                    @test full!(copy(A1)) == full(A1)
+
+                    # fill!
+                    @test full!(fill!(copy(A1), 1)) == full(t1(ones(size(A1)...)))
+
+                    # similar
+                    @test isa(similar(A1), t1)
+
+                    # Linear indexing
+                    for i = 1:length(A1)
+                        @test A1[i] == full(A1)[i]
+                    end
+
+                    # istril/istriu
+                    if uplo1 == :L
+                        @test istril(A1)
+                        @test !istriu(A1)
+                    else
+                        @test istriu(A1)
+                        @test !istril(A1)
+                    end
+
+                    # (c)transpose
+                    @test full(A1') == full(A1)'
+                    @test full(A1.') == full(A1).'
+                    @test transpose!(copy(A1)) == A1.'
+
+                    # diag
+                    @test diag(A1) == diag(full(A1))
+
+                    # real
+                    @test full(real(A1)) == real(full(A1))
+
+                    # Unary operations
+                    @test full(-A1) == -full(A1)
+
+                    # Binary operations
+                    @test full(A1 + A2) == full(A1) + full(A2)
+                    @test full(A1 - A2) == full(A1) - full(A2)
+                    @test A1*0.5 == full(A1*0.5)
+                    @test 0.5*A1 == full(0.5*A1)
+                    @test A1/0.5 == full(A1/0.5)
+                    @test 0.5\A1 == full(0.5\A1)
+
+                    # Triangular-Triangular multiplication and division
+                    elty1 != BigFloat && elty2 != BigFloat && @test_approx_eq full(A1*A2) full(A1)*full(A2)
+                    @test_approx_eq full(A1'A2) full(A1)'full(A2)
+                    @test_approx_eq full(A1*A2') full(A1)*full(A2)'
+                    @test_approx_eq full(A1'A2') full(A1)'full(A2)'
+                    @test_approx_eq full(A1/A2) full(A1)/full(A2)
+
+                    # Triangular-dense Matrix/vector multiplication
+                    @test_approx_eq A1*B[:,1] full(A1)*B[:,1]
+                    @test_approx_eq A1*B full(A1)*B
+                    @test_approx_eq A1'B[:,1] full(A1)'B[:,1]
+                    @test_approx_eq A1'B full(A1)'B
+                    @test_approx_eq A1*B' full(A1)*B'
+                    @test_approx_eq B*A1 B*full(A1)
+                    @test_approx_eq B[:,1]'A1 B[:,1]'full(A1)
+                    @test_approx_eq B'A1 B'full(A1)
+                    @test_approx_eq B*A1' B*full(A1)'
+                    @test_approx_eq B[:,1]'A1' B[:,1]'full(A1)'
+                    @test_approx_eq B'A1' B'full(A1)'
+
+                    # ... and division
+                    @test_approx_eq A1\B[:,1] full(A1)\B[:,1]
+                    @test_approx_eq A1\B full(A1)\B
+                    @test_approx_eq A1'\B[:,1] full(A1)'\B[:,1]
+                    @test_approx_eq A1'\B full(A1)'\B
+                    @test_approx_eq A1\B' full(A1)\B'
+                    @test_approx_eq A1'\B' full(A1)'\B'
+                    @test_approx_eq B/A1 B/full(A1)
+                    @test_approx_eq B/A1' B/full(A1)'
+                    @test_approx_eq B'/A1 B'/full(A1)
+                    @test_approx_eq B'/A1' B'/full(A1)'
+
+                    # inversion
+                    @test_approx_eq inv(A1) inv(lufact(full(A1)))
+
+                    # Error bounds
+                    elty1 != BigFloat && errorbounds(A1, A1\B, B)
+
+                    # Determinant
+                    @test_approx_eq_eps det(A1) det(lufact(full(A1))) sqrt(eps(real(float(one(elty1)))))*n*n
+
+                    # Matrix square root
+                    @test_approx_eq sqrtm(A1) |> t->t*t A1
+
+                    # eigenproblems
+                    if elty1 != BigFloat # Not handled yet
+                        vals, vecs = eig(A1)
+                        if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues.
+                            @test_approx_eq_eps vecs*diagm(vals)/vecs full(A1) sqrt(eps(float(real(one(vals[1])))))*(norm(A1, Inf)*n)^2
                         end
                     end
+
+                    # Condition number tests - can be VERY approximate
+                    if elty1 <:BlasFloat
+                        for p in (1.0, Inf)
+                            @test_approx_eq_eps cond(A1, p) cond(A1, p) (cond(A1, p) + cond(A1, p))
+                        end
+                    end
+
+                    if elty1 != BigFloat # Not implemented yet
+                        svd(A1)
+                        svdfact(A1)
+                        elty1 <: BlasFloat && svdfact!(copy(A1))
+                        svdvals(A1)
+                    end
                 end
             end
         end
@@ -170,14 +172,14 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
         debug && println("\ntype of A: ", eltya, " type of b: ", eltyb, "\n")
 
         debug && println("Solve upper triangular system")
-        Atri = Triangular(lufact(A)[:U], :U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
+        Atri = UpperTriangular(lufact(A)[:U]) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
         b = convert(Matrix{eltyb}, eltya <: Complex ? full(Atri)*ones(n, 2) : full(Atri)*ones(n, 2))
         x = full(Atri) \ b
 
         debug && println("Test error estimates")
         if eltya != BigFloat && eltyb != BigFloat
             for i = 1:2
-                @test  norm(x[:,1] .- 1) <= errorbounds(Triangular(A, :U), x, b)[1][i]
+                @test  norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
             end
         end
         debug && println("Test forward error [JIN 5705] if this is not a BigFloat")
@@ -198,14 +200,14 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
         end
 
         debug && println("Solve lower triangular system")
-        Atri = Triangular(lufact(A)[:U], :U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
+        Atri = UpperTriangular(lufact(A)[:U]) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
         b = convert(Matrix{eltyb}, eltya <: Complex ? full(Atri)*ones(n, 2) : full(Atri)*ones(n, 2))
         x = full(Atri)\b
 
         debug && println("Test error estimates")
         if eltya != BigFloat && eltyb != BigFloat
             for i = 1:2
-                @test  norm(x[:,1] .- 1) <= errorbounds(Triangular(A, :U), x, b)[1][i]
+                @test  norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
             end
         end
 
diff --git a/test/linalg2.jl b/test/linalg2.jl
index 3d11e4fd8bc9a..2e7a19bf2dd2a 100644
--- a/test/linalg2.jl
+++ b/test/linalg2.jl
@@ -403,7 +403,7 @@ S = sprandn(3,3,0.5)
 @test A/I == A
 @test I/λ === UniformScaling(1/λ)
 @test I\J === J
-T = Triangular(randn(3,3),:L)
+T = LowerTriangular(randn(3,3))
 @test T\I == inv(T)
 @test I\A == A
 @test A\I == inv(A)
diff --git a/test/linalg4.jl b/test/linalg4.jl
index 29187946f5aa8..acab808152d18 100644
--- a/test/linalg4.jl
+++ b/test/linalg4.jl
@@ -154,7 +154,7 @@ end
 debug && println("Test interconversion between special matrix types")
 a=[1.0:n]
 A=Diagonal(a)
-for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix]
+for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, LowerTriangular, UpperTriangular, Matrix]
     debug && println("newtype is $(newtype)")
     @test full(convert(newtype, A)) == full(A)
 end
@@ -162,29 +162,29 @@ end
 for isupper in (true, false)
     debug && println("isupper is $(isupper)")
     A=Bidiagonal(a, [1.0:n-1], isupper)
-    for newtype in [Bidiagonal, Tridiagonal, Triangular, Matrix]
+    for newtype in [Bidiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix]
         debug && println("newtype is $(newtype)")
         @test full(convert(newtype, A)) == full(A)
     end
     A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal
-    for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix]
+    for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix]
         debug && println("newtype is $(newtype)")
         @test full(convert(newtype, A)) == full(A)
     end
 end
 
-A=SymTridiagonal(a, [1.0:n-1])
+A = SymTridiagonal(a, [1.0:n-1])
 for newtype in [Tridiagonal, Matrix]
     @test full(convert(newtype, A)) == full(A)
 end
 
-A=Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal
+A = Tridiagonal(zeros(n-1), [1.0:n], zeros(n-1)) #morally Diagonal
 for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix]
     @test full(convert(newtype, A)) == full(A)
 end
 
-A=Triangular(full(Diagonal(a)), :L) #morally Diagonal
-for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix]
+A = LowerTriangular(full(Diagonal(a))) #morally Diagonal
+for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix]
     @test full(convert(newtype, A)) == full(A)
 end