spdiagm.jl

function spdiagm{T<:Integer}(B, d::Array{T, 1}, m::Integer, n::Integer)
    ndiags = length(d)
    ncoeffs = 0
    for vec in B
        ncoeffs += length(vec)
    end
    I = Array(T, ncoeffs)
    J = Array(T, ncoeffs)
    V = Array(Float64, ncoeffs)
    id = 0
    i = 0
    for vec in B
        id += 1
        diag = d[id]
        numel = length(vec)
        if diag < 0
            row = -diag
            col = 0
        elseif diag > 0
            row = 0
            col = diag
        else
            row = 0
            col = 0
        end
        range = 1+i:numel+i
        I[range] = row+1:row+numel
        J[range] = col+1:col+numel
        V[range] = vec[1:numel]
        i += numel
    end
    return sparse(I, J, V, m, n)
end

test_spdiags.jl

using Base.Test

x = rand(100)

B1 = (-1*ones(99), 2*ones(100), -1*ones(99))
B2 = (x[1:99], x, x[2:100])
B3 = (-3*ones(328), -5*ones(331), 7*ones(333), 5*ones(331), 3*ones(328))

d1 = [-1, 0, 1]
d2 = [-1, 0, 1]
d3 = [-5, -2, 0, 2, 5]

m1, n1 = 100, 100
m2, n2 = 100, 100
m3, n3 = 333, 333

@test full(spdiagm(B1, d1, m1, n1)) == full(Tridiagonal(-1*ones(99), 2*ones(100), -1*ones(99)))
@test full(spdiagm(B2, d2, m2, n2)) == full(Tridiagonal(x[1:99], x, x[2:100]))
@test squeeze(full(spdiagm(B3, d3, m3, n3)[6, 1:11])', 2) == [-3.0, 0.0, 0.0, -5.0, 0.0, 7.0, 0.0, 5.0, 0.0, 0.0, 3.0]