Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Address all issues discussed in JuliaLang#4547
Add test for spdiagm and update docs
  • Loading branch information
Viral B. Shah committed Oct 17, 2013
1 parent 7c59a49 commit 82d96bd
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 24 deletions.
46 changes: 24 additions & 22 deletions base/sparse/sparsematrix.jl
Expand Up @@ -1288,32 +1288,19 @@ function istril(A::SparseMatrixCSC)
return true
end

function spdiagm{T}(v::Union(AbstractVector{T},AbstractMatrix{T}), k::Integer=0)
if isa(v, AbstractMatrix)
if (size(v,1) != 1 && size(v,2) != 1)
error("Input should be nx1 or 1xn")
end
end

nv = length(v)
nr = nc = nv + abs(k)

x = diagind(nr, nc, k)
I, J = ind2sub((nr,nc), x)
# Create a sparse diagonal matrix by specifying multiple diagonals
# packed into a tuple, alongside their diagonal offsets and matrix shape

sparse(I,J,v,nr,nc)
end

## extend spdiagm by specifying multiple diagonals packed into a tuple alongside their diagonal offsets and matrix shape
function spdiagm{T<:Integer}(B, d::Array{T, 1}, m::Integer, n::Integer)
function spdiagm_int(B::Tuple, d::Tuple)
ndiags = length(d)
if length(B) != ndiags; throw(ArgumentError); end
ncoeffs = 0
for vec in B
ncoeffs += length(vec)
end
I = Array(T, ncoeffs)
J = Array(T, ncoeffs)
V = Array(Float64, ncoeffs)
I = Array(Int, ncoeffs)
J = Array(Int, ncoeffs)
V = Array(eltype(B[1]), ncoeffs)
id = 0
i = 0
for vec in B
Expand All @@ -1333,12 +1320,27 @@ function spdiagm{T<:Integer}(B, d::Array{T, 1}, m::Integer, n::Integer)
range = 1+i:numel+i
I[range] = row+1:row+numel
J[range] = col+1:col+numel
V[range] = vec[1:numel]
copy!(sub(V, range), vec)
i += numel
end
return sparse(I, J, V, m, n)

return (I,J,V)
end

function spdiagm(B::Tuple, d::Tuple, m::Integer, n::Integer)
(I,J,V) = spdiagm_int(B, d)
return sparse(I,J,V,m,n)
end

function spdiagm(B::Tuple, d::Tuple)
(I,J,V) = spdiagm_int(B, d)
return sparse(I,J,V)
end

spdiagm(B::AbstractVector, d::Number, m, n) = spdiagm((B,), (d,), m, n)

spdiagm(B::AbstractVector, d::Number) = spdiagm((B,), (d,))

## expand a colptr or rowptr into a dense index vector
function expandptr{T<:Integer}(V::Vector{T})
if V[1] != 1 error("expandptr: first index must be one") end
Expand Down
4 changes: 2 additions & 2 deletions doc/stdlib/sparse.rst
Expand Up @@ -51,9 +51,9 @@ Sparse matrices support much of the same set of operations as dense matrices. Th

Create a sparse identity matrix of specified type of size ``m x m``. In case ``n`` is supplied, create a sparse identity matrix of size ``m x n``.

.. function:: spdiagm(v)
.. function:: spdiagm(B, d[, m, n])

Construct a sparse diagonal matrix and place "v" on the diagonal.
Construct a sparse diagonal matrix. ``B`` is a tuple of vectors containing the diagonals and ``d`` is a tuple containing the positions of the diagonals. In the case the input contains only one diagonaly, ``B`` can be a vector (instead of a tuple) and ``d`` can be the diagonal position (instead of a tuple). Optionally, ``m`` and ``n`` specify the size of the resulting sparse matrix.

.. function:: sprand(m,n,density[,rng])

Expand Down
4 changes: 4 additions & 0 deletions test/sparse.jl
Expand Up @@ -76,6 +76,10 @@ end
@test prod(se33, 1) == [0.0 0.0 0.0]
@test prod(se33, 2) == [0.0 0.0 0.0]'

# spdiagm
@test full(spdiagm((ones(2), ones(2)), (0, -1), 3, 3)) ==
[1.0 0.0 0.0; 1.0 1.0 0.0; 0.0 1.0 0.0]

# elimination tree
## upper triangle of the pattern test matrix from Figure 4.2 of
## "Direct Methods for Sparse Linear Systems" by Tim Davis, SIAM, 2006
Expand Down

0 comments on commit 82d96bd

Please sign in to comment.