diff --git a/base/export.jl b/base/export.jl
index 7858a639739dd..c81a82ca00f65 100644
--- a/base/export.jl
+++ b/base/export.jl
@@ -81,15 +81,20 @@ export
     SubOrDArray,
     SubString,
     TransformedString,
+    Tridiagonal,
     VecOrMat,
     Vector,
     VersionNumber,
     WeakKeyDict,
+    Woodbury,
     Zip,
     Stat,
     Factorization,
     Cholesky,
     LU,
+    LUTridiagonal,
+    LDLT,
+    LDLTTridiagonal,
     QR,
     QRP,
     
@@ -634,6 +639,7 @@ export
     randsym,
     rank,
     rref,
+    solve,
     svd,
     svdvals,
     trace,
diff --git a/base/factorizations.jl b/base/factorizations.jl
index 2ae7b7c10ec68..6311c5e35fbd2 100644
--- a/base/factorizations.jl
+++ b/base/factorizations.jl
@@ -281,3 +281,48 @@ end
 ##ToDo:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
 ##       Add rcond methods for Cholesky, LU, QR and QRP types
 ## Lower priority: Add LQ, QL and RQ factorizations
+
+
+#### Factorizations for Tridiagonal ####
+type LDLTTridiagonal{T} <: Factorization{T}
+    D::Vector{T}
+    E::Vector{T}
+end
+function LDLTTridiagonal{T<:LapackScalar}(A::Tridiagonal{T})
+    D = copy(A.d)
+    E = copy(A.dl)
+    _jl_lapack_pttrf(D, E)
+    LDLTTridiagonal(D, E)
+end
+LDLT(A::Tridiagonal) = LDLTTridiagonal(A)
+
+(\){T<:LapackScalar}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) =
+    _jl_lapack_pttrs(C.D, C.E, copy(B))
+
+type LUTridiagonal{T} <: Factorization{T}
+    lu::Tridiagonal{T}
+    ipiv::Vector{Int32}
+    function LUTridiagonal(lu::Tridiagonal{T}, ipiv::Vector{Int32})
+        m, n = size(lu)
+        m == numel(ipiv) ? new(lu, ipiv) : error("LU: dimension mismatch")
+    end
+end
+show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))
+
+function LU{T<:LapackScalar}(A::Tridiagonal{T})
+    lu, ipiv = _jl_lapack_gttrf(copy(A))
+    LUTridiagonal{T}(lu, ipiv)
+end
+
+function lu(A::Tridiagonal)
+    error("lu(A) is not defined when A is Tridiagonal. Use LU(A) instead.")
+end
+
+function det(lu::LUTridiagonal)
+    prod(lu.lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -1 : 1)
+end
+
+det(A::Tridiagonal) = det(LU(A))
+
+(\){T<:LapackScalar}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
+    _jl_lapack_gttrs('N', lu.lu, lu.ipiv, copy(B))
diff --git a/base/linalg.jl b/base/linalg.jl
index ab4c6cf51164f..c68fefd0f799d 100644
--- a/base/linalg.jl
+++ b/base/linalg.jl
@@ -1,4 +1,6 @@
-## linalg.jl: Basic Linear Algebra interface specifications ##
+## linalg.jl: Basic Linear Algebra interface specifications and
+## specialized matrix types
+
 #
 # This file mostly contains commented functions which are supposed
 # to be defined in type-specific linalg_<type>.jl files.
diff --git a/base/linalg_lapack.jl b/base/linalg_lapack.jl
index f0f0bac74683d..db655ef84dace 100644
--- a/base/linalg_lapack.jl
+++ b/base/linalg_lapack.jl
@@ -885,3 +885,187 @@ end
 expm{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) = expm!(copy(A))
 expm{T<:Integer}(A::StridedMatrix{T}) = expm!(float(A))
 
+
+#### Tridiagonal matrix routines ####
+function \{T<:LapackScalar}(M::Tridiagonal{T}, rhs::StridedVecOrMat{T})
+    if stride(rhs, 1) == 1
+        x = copy(rhs)
+        Mc = copy(M)
+        Mlu, x = _jl_lapack_gtsv(Mc, x)
+        return x
+    end
+    solve(M, rhs)  # use the Julia "fallback"
+end
+
+eig(M::Tridiagonal) = _jl_lapack_stev('V', copy(M))
+
+# Decompositions
+for (gttrf, pttrf, elty) in
+    ((:dgttrf_,:dpttrf_,:Float64),
+     (:sgttrf_,:spttrf_,:Float32),
+     (:zgttrf_,:zpttrf_,:Complex128),
+     (:cgttrf_,:cpttrf_,:Complex64))
+    @eval begin
+        function _jl_lapack_gttrf(M::Tridiagonal{$elty})
+            info = zero(Int32)
+            n    = int32(length(M.d))
+            ipiv = Array(Int32, n)
+            ccall(dlsym(_jl_liblapack, $string(gttrf)),
+                  Void,
+                  (Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{Int32}),
+                  &n, M.dl, M.d, M.du, M.dutmp, ipiv, &info)
+            if info != 0 throw(LapackException(info)) end
+            M, ipiv
+        end
+        function _jl_lapack_pttrf(D::Vector{$elty}, E::Vector{$elty})
+            info = zero(Int32)
+            n    = int32(length(D))
+            if length(E) != n-1
+                error("subdiagonal must be one element shorter than diagonal")
+            end
+            ccall(dlsym(_jl_liblapack, $string(pttrf)),
+                  Void,
+                  (Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32}),
+                  &n, D, E, &info)
+            if info != 0 throw(LapackException(info)) end
+            D, E
+        end
+    end
+end
+# Direct solvers
+for (gtsv, ptsv, elty) in
+    ((:dgtsv_,:dptsv_,:Float64),
+     (:sgtsv_,:sptsv,:Float32),
+     (:zgtsv_,:zptsv,:Complex128),
+     (:cgtsv_,:cptsv,:Complex64))
+    @eval begin
+        function _jl_lapack_gtsv(M::Tridiagonal{$elty}, B::StridedVecOrMat{$elty})
+            if stride(B,1) != 1
+                error("_jl_lapack_gtsv: matrix columns must have contiguous elements");
+            end
+            info = zero(Int32)
+            n    = int32(length(M.d))
+            nrhs = int32(size(B, 2))
+            ldb  = int32(stride(B, 2))
+            ccall(dlsym(_jl_liblapack, $string(gtsv)),
+                  Void,
+                  (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{Int32}),
+                  &n, &nrhs, M.dl, M.d, M.du, B, &ldb, &info)
+            if info != 0 throw(LapackException(info)) end
+            M, B
+        end
+        function _jl_lapack_ptsv(M::Tridiagonal{$elty}, B::StridedVecOrMat{$elty})
+            if stride(B,1) != 1
+                error("_jl_lapack_ptsv: matrix columns must have contiguous elements");
+            end
+            info = zero(Int32)
+            n    = int32(length(M.d))
+            nrhs = int32(size(B, 2))
+            ldb  = int32(stride(B, 2))
+            ccall(dlsym(_jl_liblapack, $string(ptsv)),
+                  Void,
+                  (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{Int32}),
+                  &n, &nrhs, M.d, M.dl, B, &ldb, &info)
+            if info != 0 throw(LapackException(info)) end
+            M, B
+        end
+    end
+end
+# Solvers using decompositions
+for (gttrs, pttrs, elty) in
+    ((:dgttrs_,:dpttrs_,:Float64),
+     (:sgttrs_,:spttrs,:Float32),
+     (:zgttrs_,:zpttrs,:Complex128),
+     (:cgttrs_,:cpttrs,:Complex64))
+    @eval begin
+        function _jl_lapack_gttrs(trans::LapackChar, M::Tridiagonal{$elty}, ipiv::Vector{Int32}, B::StridedVecOrMat{$elty})
+            if stride(B,1) != 1
+                error("_jl_lapack_gttrs: matrix columns must have contiguous elements");
+            end
+            info = zero(Int32)
+            n    = int32(length(M.d))
+            nrhs = int32(size(B, 2))
+            ldb  = int32(stride(B, 2))
+            ccall(dlsym(_jl_liblapack, $string(gttrs)),
+                  Void,
+                  (Ptr{Uint8}, Ptr{Int32}, Ptr{Int32},
+                   Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, Ptr{Int32}),
+                  &trans, &n, &nrhs, M.dl, M.d, M.du, M.dutmp, ipiv, B, &ldb, &info)
+            if info != 0 throw(LapackException(info)) end
+            B
+        end
+        function _jl_lapack_pttrs(D::Vector{$elty}, E::Vector{$elty}, B::StridedVecOrMat{$elty})
+            if stride(B,1) != 1
+                error("_jl_lapack_pttrs: matrix columns must have contiguous elements");
+            end
+            info = zero(Int32)
+            n    = int32(length(D))
+            if length(E) != n-1
+                error("subdiagonal must be one element shorter than diagonal")
+            end
+            nrhs = int32(size(B, 2))
+            ldb  = int32(stride(B, 2))
+            ccall(dlsym(_jl_liblapack, $string(pttrs)),
+                  Void,
+                  (Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{Int32}),
+                  &n, &nrhs, D, E, B, &ldb, &info)
+            if info != 0 throw(LapackException(info)) end
+            B
+        end
+    end
+end
+# Eigenvalue-eigenvector (symmetric only)
+for (stev, elty) in
+    ((:dstev_,:Float64),
+     (:sstev_,:Float32),
+     (:zstev_,:Complex128),
+     (:cstev_,:Complex64))
+    @eval begin
+        function _jl_lapack_stev(Z::Array, M::Tridiagonal{$elty})
+            n    = int32(length(M.d))
+            if isempty(Z)
+                job = 'N'
+                ldz = 1
+                work = Array($elty, 0)
+                Ztmp = work
+            else
+                if stride(Z,1) != 1
+                    error("_jl_lapack_stev: eigenvector matrix columns must have contiguous elements");
+                end
+                if size(Z, 1) != n
+                    error("_jl_lapack_stev: eigenvector matrix columns are not of the correct size")
+                end
+                Ztmp = Z
+                job = 'V'
+                ldz  = int32(stride(Z, 2))
+                work = Array($elty, max(1, 2*n-2))
+            end
+            info = zero(Int32)
+            ccall(dlsym(_jl_liblapack, $string(stev)),
+                  Void,
+                  (Ptr{Uint8}, Ptr{Int32},
+                   Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
+                   Ptr{Int32}, Ptr{$elty}, Ptr{Int32}),
+                  &job, &n, M.d, M.dl, Ztmp, &ldz, work, &info)
+            if info != 0 throw(LapackException(info)) end
+            M.d
+        end
+    end
+end
+function _jl_lapack_stev(job::LapackChar, M::Tridiagonal)
+    if job == 'N' || job == 'n'
+        Z = []
+    elseif job == 'V' || job == 'v'
+        n = length(M.d)
+        Z = Array(eltype(M), n, n)
+    else
+        error("Job type not recognized")
+    end
+    D = _jl_lapack_stev(Z, M)
+    return D, Z
+end
diff --git a/base/linalg_specialized.jl b/base/linalg_specialized.jl
new file mode 100644
index 0000000000000..349db1c8cd2c8
--- /dev/null
+++ b/base/linalg_specialized.jl
@@ -0,0 +1,299 @@
+#### Specialized matrix types ####
+
+# Some of these also have important routines defined in factorizations.jl,
+# linalg_lapack.jl, etc.
+
+#### Tridiagonal matrices ####
+type Tridiagonal{T} <: AbstractMatrix{T}
+    dl::Vector{T}   # sub-diagonal
+    d::Vector{T}   # diagonal
+    du::Vector{T}   # sup-diagonal
+    dutmp::Vector{T}  # scratch space for vector RHS solver, sup-diagonal
+    rhstmp::Vector{T}  # scratch space, rhs
+
+    function Tridiagonal(N::Int)
+        dutmp = Array(T, N-1)
+        rhstmp = Array(T, N)
+        new(dutmp, rhstmp, dutmp, dutmp, rhstmp)  # first three will be overwritten
+    end
+end
+function Tridiagonal{T<:Number}(dl::Vector{T}, d::Vector{T}, du::Vector{T})
+    N = length(d)
+    if length(dl) != N-1 || length(du) != N-1
+        error("The sub- and super-diagonals must have length N-1")
+    end
+    M = Tridiagonal{T}(N)
+    M.dl = copy(dl)
+    M.d = copy(d)
+    M.du = copy(du)
+    return M
+end
+function Tridiagonal{Tl<:Number, Td<:Number, Tu<:Number}(dl::Vector{Tl}, d::Vector{Td}, du::Vector{Tu})
+    R = promote(Tl, Td, Tu)
+    Tridiagonal(convert(Vector{R}, dl), convert(Vector{R}, d), convert(Vector{R}, du))
+end
+size(M::Tridiagonal) = (length(M.d), length(M.d))
+function show(io, M::Tridiagonal)
+    println(io, summary(M), ":")
+    print(io, " sub: ")
+    print_matrix(io, (M.dl)')
+    print(io, "\ndiag: ")
+    print_matrix(io, (M.d)')
+    print(io, "\n sup: ")
+    print_matrix(io, (M.du)')
+end
+full{T}(M::Tridiagonal{T}) = convert(Matrix{T}, M)
+function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
+    A = zeros(T, size(M))
+    for i = 1:length(M.d)
+        A[i,i] = M.d[i]
+    end
+    for i = 1:length(M.d)-1
+        A[i+1,i] = M.dl[i]
+        A[i,i+1] = M.du[i]
+    end
+    return A
+end
+function similar(M::Tridiagonal, T, dims::Dims)
+    if length(dims) != 2 || dims[1] != dims[2]
+        error("Tridiagonal matrices must be square")
+    end
+    return Tridiagonal{T}(dims[1])
+end
+copy(M::Tridiagonal) = Tridiagonal(M.dl, M.d, M.du)
+
+# Operations on Tridiagonal matrices
+round(M::Tridiagonal) = Tridiagonal(round(M.dl), round(M.d), round(M.du))
+iround(M::Tridiagonal) = Tridiagonal(iround(M.dl), iround(M.d), iround(M.du))
+
+## Solvers
+# Allocation-free variants
+# Note that solve is non-aliasing, so you can use the same array for
+# input and output
+function solve(x::AbstractArray, xrng::Ranges{Int}, M::Tridiagonal, rhs::AbstractArray, rhsrng::Ranges{Int})
+    d = M.d
+    N = length(d)
+    if length(xrng) != N || length(rhsrng) != N
+        error("dimension mismatch")
+    end
+    dl = M.dl
+    du = M.du
+    dutmp = M.dutmp
+    rhstmp = M.rhstmp
+    xstart = first(xrng)
+    xstride = step(xrng)
+    rhsstart = first(rhsrng)
+    rhsstride = step(rhsrng)
+    # Forward sweep
+    denom = d[1]
+    dulast = du[1] / denom
+    dutmp[1] = dulast
+    rhslast = rhs[rhsstart] / denom
+    rhstmp[1] = rhslast
+    irhs = rhsstart+rhsstride
+    for i in 2:N-1
+        dltmp = dl[i-1]
+        denom = d[i] - dltmp*dulast
+        dulast = du[i] / denom
+        dutmp[i] = dulast
+        rhslast = (rhs[irhs] - dltmp*rhslast)/denom
+        rhstmp[i] = rhslast
+        irhs += rhsstride
+    end
+    dltmp = dl[N-1]
+    denom = d[N] - dltmp*dulast
+    xlast = (rhs[irhs] - dltmp*rhslast)/denom
+    # Backward sweep
+    ix = xstart + (N-2)*xstride
+    x[ix+xstride] = xlast
+    for i in N-1:-1:1
+        xlast = rhstmp[i] - dutmp[i]*xlast
+        x[ix] = xlast
+        ix -= xstride
+    end
+    return x
+end
+solve(x::StridedVector, M::Tridiagonal, rhs::StridedVector) = solve(x, 1:length(x), M, rhs, 1:length(rhs))
+function solve(M::Tridiagonal, rhs::StridedVector)
+    x = similar(rhs)
+    solve(x, M, rhs)
+end
+function solve(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
+    if size(B, 1) != size(M, 1)
+        error("dimension mismatch")
+    end
+    if size(X) != size(B)
+        error("dimension mismatch in output")
+    end
+    m, n = size(B)
+    r = 1:m
+    for j = 1:n
+        r.start = (j-1)*m+1
+        solve(X, r, M, B, r)
+    end
+    return X
+end
+function solve(M::Tridiagonal, B::StridedMatrix)
+    X = similar(B)
+    solve(X, M, B)
+end
+
+# User-friendly solver
+\(M::Tridiagonal, rhs::Union(StridedVector,StridedMatrix)) = solve(M, rhs)
+
+# Tridiagonal multiplication
+function mult(x::AbstractArray, xrng::Ranges{Int}, M::Tridiagonal, v::AbstractArray, vrng::Ranges{Int})
+    dl = M.dl
+    d = M.d
+    du = M.du
+    N = length(d)
+    xi = first(xrng)
+    xstride = step(xrng)
+    vi = first(vrng)
+    vstride = step(vrng)
+    x[xi] = d[1]*v[vi] + du[1]*v[vi+vstride]
+    xi += xstride
+    for i = 2:N-1
+        x[xi] = dl[i-1]*v[vi] + d[i]*v[vi+vstride] + du[i]*v[vi+2*vstride]
+        xi += xstride
+        vi += vstride
+    end
+    x[xi] = dl[N-1]*v[vi] + d[N]*v[vi+vstride]
+    return x
+end
+mult(x::StridedVector, M::Tridiagonal, v::StridedVector) = mult(x, 1:length(x), M, v, 1:length(v))
+function mult(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
+    if size(B, 1) != size(M, 1)
+        error("dimension mismatch")
+    end
+    if size(X) != size(B)
+        error("dimension mismatch in output")
+    end
+    m, n = size(B)
+    r = 1:m
+    for j = 1:n
+        r.start = (j-1)*m+1
+        solve(X, r, M, B, r)
+    end
+    return X
+end
+mult(X::StridedMatrix, M1::Tridiagonal, M2::Tridiagonal) = mult(X, M1, full(M2))
+function *(M::Tridiagonal, B::Union(StridedVector,StridedMatrix))
+    X = similar(B)
+    mult(X, M, B)
+end
+
+
+
+#### Woodbury matrices ####
+# This type provides support for the Woodbury matrix identity
+type Woodbury{T} <: AbstractMatrix{T}
+    A
+    U::Matrix{T}
+    C
+    Cp
+    V::Matrix{T}
+    tmpN1::Vector{T}
+    tmpN2::Vector{T}
+    tmpk1::Vector{T}
+    tmpk2::Vector{T}
+
+    function Woodbury(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T})
+        N = size(A, 1)
+        k = size(U, 2)
+        if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
+            error("Sizes do not match")
+        end
+        if k > 1
+            if size(C, 1) != k || size(C, 2) != k
+                error("Size of C is incorrect")
+            end
+        end
+        Cp = inv(inv(C) + V*(A\U))
+        # temporary space for allocation-free solver
+        tmpN1 = Array(T, N)
+        tmpN2 = Array(T, N)
+        tmpk1 = Array(T, k)
+        tmpk2 = Array(T, k)
+        # don't copy A, it could be huge
+        new(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
+    end
+end
+Woodbury{T}(A::AbstractMatrix{T}, U::Matrix{T}, C, V::Matrix{T}) = Woodbury{T}(A, U, C, V)
+Woodbury{T}(A::AbstractMatrix{T}, U::Vector{T}, C, V::Matrix{T}) = Woodbury{T}(A, reshape(U, length(U), 1), C, V)
+
+size(W::Woodbury) = size(W.A)
+function show(io, W::Woodbury)
+    println(io, summary(W), ":")
+    print(io, "A: ", W.A)
+    print(io, "\nU:\n")
+    print_matrix(io, W.U)
+    if isa(W.C, Matrix)
+        print(io, "\nC:\n")
+        print_matrix(io, W.C)
+    else
+        print(io, "\nC: ", W.C)
+    end
+    print(io, "\nV:\n")
+    print_matrix(io, W.V)
+end
+full{T}(W::Woodbury{T}) = convert(Matrix{T}, W)
+convert{T}(::Type{Matrix{T}}, W::Woodbury{T}) = full(W.A) + W.U*W.C*W.V
+function similar(W::Woodbury, T, dims::Dims)
+    if length(dims) != 2 || dims[1] != dims[2]
+        error("Woodbury matrices must be square")
+    end
+    n = size(W, 1)
+    k = size(W.U, 2)
+    return Woodbury{T}(similar(W.A), Array(T, n, k), Array(T, k, k), Array(T, k, n))
+end
+copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)
+
+## Woodbury matrix routines ##
+
+function *(W::Woodbury, B::StridedVecOrMat)
+    return W.A*B + W.U*(W.C*(W.V*B))
+end
+function \(W::Woodbury, R::StridedVecOrMat)
+    AinvR = W.A\R
+    return AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
+end
+# Allocation-free solver for arbitrary strides (requires that W.A has a
+# non-aliasing "solve" routine, e.g., is Tridiagonal)
+function solve(x::AbstractArray, xrng::Ranges{Int}, W::Woodbury, rhs::AbstractArray, rhsrng::Ranges{Int})
+    solve(W.tmpN1, 1:length(W.tmpN1), W.A, rhs, rhsrng)
+    A_mul_B(W.tmpk1, W.V, W.tmpN1)
+    A_mul_B(W.tmpk2, W.Cp, W.tmpk1)
+    A_mul_B(W.tmpN2, W.U, W.tmpk2)
+    solve(W.tmpN2, W.A, W.tmpN2)
+    indx = first(xrng)
+    xinc = step(xrng)
+    for i = 1:length(W.tmpN2)
+        x[indx] = W.tmpN1[i] - W.tmpN2[i]
+        indx += xinc
+    end
+end
+solve(x::AbstractVector, W::Woodbury, rhs::AbstractVector) = solve(x, 1:length(x), W, rhs, 1:length(rhs))
+function solve(W::Woodbury, rhs::AbstractVector)
+    x = similar(rhs)
+    solve(x, W, rhs)
+end
+function solve(X::StridedMatrix, W::Woodbury, B::StridedMatrix)
+    if size(B, 1) != size(W, 1)
+        error("dimension mismatch")
+    end
+    if size(X) != size(B)
+        error("dimension mismatch in output")
+    end
+    m, n = size(B)
+    r = 1:m
+    for j = 1:n
+        r.start = (j-1)*m+1
+        solve(X, r, W, B, r)
+    end
+    return X
+end
+function solve(W::Woodbury, B::StridedMatrix)
+    X = similar(B)
+    solve(X, W, B)
+end
diff --git a/base/sysimg.jl b/base/sysimg.jl
index f6746abb661ff..8b95ea9ba7198 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -130,6 +130,7 @@ include("build_h.jl")
 # linear algebra
 include("linalg.jl")
 include("linalg_dense.jl")
+include("linalg_specialized.jl")
 include("linalg_blas.jl")
 include("linalg_lapack.jl")
 include("factorizations.jl")
diff --git a/extras/linalg_sparse.jl b/extras/linalg_sparse.jl
index 62796500c6219..940b0fc2059ce 100644
--- a/extras/linalg_sparse.jl
+++ b/extras/linalg_sparse.jl
@@ -572,64 +572,3 @@ diagmm{Tv,Ti,T}(A::SparseMatrixCSC{Tv,Ti}, b::Vector{T}) =
 
 diagmm{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) =
     diagmm!(SparseMatrixCSC(size(A,1),size(A,2),Ti[],Ti[],promote_type(Tv,T)[]), b, A)
-
-
-# Tridiagonal solver
-# Allocation-free variants
-function solve(x::Array, xstart::Int, xstride::Int, M::Tridiagonal, d::Array, dstart::Int, dstride::Int)
-    # Grab refs to members (for efficiency)
-    a = M.a
-    b = M.b
-    c = M.c
-    cp = M.cp
-    dp = M.dp
-    N = length(b)
-    # Forward sweep
-    cp[1] = c[1] / b[1]
-    dp[1] = d[dstart] / b[1]
-    id = dstart+dstride
-    for i = 2:N
-        atmp = a[i]
-        temp = b[i] - atmp*cp[i-1]
-        cp[i] = c[i] / temp
-        dp[i] = (d[id] - atmp*dp[i-1])/temp
-        id += dstride
-    end
-    # Backward sweep
-    ix = xstart + (N-2)*xstride
-    x[ix+xstride] = dp[N]
-    for i = N-1:-1:1
-        x[ix] = dp[i] - cp[i]*x[ix+xstride]
-        ix -= xstride
-    end
-end
-function solve(x::Vector, M::Tridiagonal, d::Vector)
-    if length(d) != length(M.b)
-        error("Size mismatch between matrix and rhs")
-    end
-    solve(x, 1, 1, M, d, 1, 1)
-end
-# User-friendly solver
-function \(M::Tridiagonal, d::Vector)
-    x = similar(d)
-    solve(x, M, d)
-    return x
-end
-
-# Tridiagonal multiplication
-function mult(x::Vector, M::Tridiagonal, v::Vector)
-    a = M.a
-    b = M.b
-    c = M.c
-    N = length(b)
-    x[1] = b[1]*v[1] + c[1]*v[2]
-    for i = 2:N-1
-        x[i] = a[i]*v[i-1] + b[i]*v[i] + c[i]*v[i+1]
-    end
-    x[N] = a[N]*v[N-1] + b[N]*v[n]
-end
-function *(M::Tridiagonal, v::Vector)
-    x = similar(v)
-    mult(x, M, v)
-    return x
-end
diff --git a/extras/sparse.jl b/extras/sparse.jl
index ad18d6db847e9..f1c97d7f59f2e 100644
--- a/extras/sparse.jl
+++ b/extras/sparse.jl
@@ -1023,48 +1023,3 @@ function assign(S::SparseAccumulator, v, i::Integer)
     end
     return S
 end
-
-type Tridiagonal{T<:FloatingPoint} <: AbstractMatrix{T}
-    a::Vector{T}   # sub-diagonal
-    b::Vector{T}   # diagonal
-    c::Vector{T}   # sup-diagonal
-    cp::Vector{T}  # scratch space, sup-diagonal
-    dp::Vector{T}  # scratch space, rhs
-
-    function Tridiagonal(N::Int)
-        cp = Array(T, N)
-        dp = Array(T, N)
-        new(cp, cp, cp, cp, dp)  # first three will be overwritten
-    end
-end
-function Tridiagonal{T}(a::Vector{T}, b::Vector{T}, c::Vector{T})
-    N = length(b)
-    if length(a) != N || length(c) != N
-        error("All three vectors must have the same length")
-    end
-    M = Tridiagonal{T}(N)
-    M.a = copy(a)
-    M.b = copy(b)
-    M.c = copy(c)
-    return M
-end
-size(M::Tridiagonal) = (length(M.b), length(M.b))
-function show(io, M::Tridiagonal)
-    println(io, summary(M), ":")
-    print_matrix(io, vcat((M.a)', (M.b)', (M.c)'))
-#    println(io, " sub: ", (M.a)')
-#    println(io, "diag: ", (M.b)')
-#    println(io, " sup: ", (M.c)')
-end
-full{T}(M::Tridiagonal{T}) = convert(Matrix{T}, M)
-function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
-    A = zeros(T, size(M))
-    for i = 1:length(M.b)
-        A[i,i] = M.b[i]
-    end
-    for i = 1:length(M.b)-1
-        A[i+1,i] = M.a[i+1]
-        A[i,i+1] = M.c[i]
-    end
-    return A
-end
diff --git a/test/lapack.jl b/test/lapack.jl
index 916742749e66c..37c2d7a1e3b7e 100644
--- a/test/lapack.jl
+++ b/test/lapack.jl
@@ -140,4 +140,60 @@ eA3 = [-1.50964415879218 -5.6325707998812  -4.934938326092;
         0.367879439109187 1.47151775849686  1.10363831732856;
         0.135335281175235 0.406005843524598 0.541341126763207]'
 @assert norm((expm(A3) - eA3) ./ eA3) < 50000*eps()
+
+
+# basic tridiagonal operations
+n = 5
+d = 1 + rand(n)
+dl = -rand(n-1)
+du = -rand(n-1)
+T = Tridiagonal(dl, d, du)
+@assert size(T, 1) == n
+@assert size(T) == (n, n)
+F = diagm(d)
+for i = 1:n-1
+    F[i,i+1] = du[i]
+    F[i+1,i] = dl[i]
+end
+@assert full(T) == F
+
+# tridiagonal linear algebra
+Eps = sqrt(eps())
+v = randn(n)
+@assert T*v == F*v
+invFv = F\v
+@assert norm(T\v - invFv) < Eps
+@assert norm(solve(T,v) - invFv) < Eps
+B = randn(n,2)
+@assert norm(solve(T, B) - F\B) < Eps
+Tlu = LU(T)
+x = Tlu\v
+@assert norm(x - invFv) < Eps
+
+# symmetric tridiagonal
+Ts = Tridiagonal(dl, d, dl)
+Fs = full(Ts)
+invFsv = Fs\v
+Tldlt = LDLT(Ts)
+x = Tldlt\v
+@assert norm(x - invFsv) < Eps
+
+# eigenvalues/eigenvectors of symmetric tridiagonal
+DT, VT = eig(Ts)
+D, V = eig(Fs)
+@assert norm(DT - D) < Eps
+@assert norm(VT - V) < Eps
+
+
+# Woodbury
+U = randn(n,2)
+V = randn(2,n)
+C = randn(2,2)
+W = Woodbury(T, U, C, V)
+F = full(W)
+
+@assert norm(W*v - F*v) < Eps
+@assert norm(W\v - F\v) < Eps
+
+
 end
diff --git a/test/sparse.jl b/test/sparse.jl
index afa168a1202bd..86c97e372d702 100644
--- a/test/sparse.jl
+++ b/test/sparse.jl
@@ -41,5 +41,4 @@ for i = 1 : 10
     a = sprand(5, 4, 0.5)
     @assert all([a[1:2,1:2] a[1:2,3:4]; a[3:5,1] [a[3:4,2:4]; a[5,2:4]]] == a)
 end
-
 end # cd