In [1]:
## Intro to Sparse Matrices

In [2]:
## Coordinate format.
I = [1,2,1,2,5,1,4,6,3]
J = [1,1,2,2,1,4,2,2,3]
V = [2,-1,1,1,-1,-1,1,2,2]

Out[2]:
9-element Array{Int64,1}:
2
-1
1
1
-1
-1
1
2
2
In [3]:
##
A = zeros(Int64,6,4)
for k = 1:length(I)
A[I[k],J[k]] = V[k]
end

In [4]:
##
A = sparse(I,J,V,6,4)

Out[4]:
6×4 SparseMatrixCSC{Int64,Int64} with 9 stored entries:
[1, 1]  =  2
[2, 1]  =  -1
[5, 1]  =  -1
[1, 2]  =  1
[2, 2]  =  1
[4, 2]  =  1
[6, 2]  =  2
[3, 3]  =  2
[1, 4]  =  -1
In [5]:
##
"""
spmatvec_coord Compute a sparse matrix vector product with A in coord form
y = spmatvec(Acoord, x) compute y=A*x where A is represented as a matrix in
coordinate form:
Acoord is a nnz-by-3 array where
Acoord[k,1] is the row index of the kth non-zero
Acoord[k,2] is the column index of the kth non-zero
Acoord[k,3] is the value of the kth non-zero

Example:
n = 10
A = sprand(n,n,0.2)
i,j,v = findnz(A)
x = randn(n)
z = A*x
C = SparseTripletFormat(i,j,v)
y = spmatvec_coord(C,x)
@show norm(z-y)


"""
type SparseTripletFormat{T}
i::Array{Int64}
j::Array{Int64}
v::Array{T}
end
function spmatvec_coord(C, x)
m = maximum(C.i)
y = zeros(m)
for k = 1:length(C.i)
y[C.i[k]] = y[C.i[k]] + C.v[k]*x[C.j[k]]
end
return y
end

Out[5]:
spmatvec_coord (generic function with 1 method)
In [6]:
##
"""
spmatvec_csc Compute a sparse matrix vector product with A in CSC form
y = spmatvec_csc(A,x) multiply a sparse matrix stored in compressed


Example:
n = 10
A = sprand(n,n,0.2)
x = randn(n)
z = A*x
y = spmatvec_csc(A,x)
@show norm(z-y)

"""
function spmatvec_csc{T}(A::SparseMatrixCSC{T,Int64},x)
assert(A.n == length(x))
y = zeros(A.m)
for j = 1:A.n
for k = A.colptr[j]:(A.colptr[j+1]-1)
y[A.rowval[k]] = y[A.rowval[k]] + A.nzval[k]*x[j]
end
end
return y
end

Out[6]:
spmatvec_csc
In [7]:
##
x = randn(4,1)
C = SparseTripletFormat(I,J,V)
y = spmatvec_coord(C,x)

Out[7]:
6-element Array{Float64,1}:
-1.2844
-1.52415
1.35238
-1.5007
-0.0234462
-3.0014   
In [8]:
##
yjulia = A*x

Out[8]:
6×1 Array{Float64,2}:
-1.2844
-1.52415
1.35238
-1.5007
-0.0234462
-3.0014   
In [9]:
##
y - yjulia

Out[9]:
6×1 Array{Float64,2}:
0.0
0.0
0.0
0.0
0.0
0.0
In [10]:
## Convert to a sparse matrix
A = sparse(I,J,V,6,4)

Out[10]:
6×4 SparseMatrixCSC{Int64,Int64} with 9 stored entries:
[1, 1]  =  2
[2, 1]  =  -1
[5, 1]  =  -1
[1, 2]  =  1
[2, 2]  =  1
[4, 2]  =  1
[6, 2]  =  2
[3, 3]  =  2
[1, 4]  =  -1
In [11]:
## Convert to a dense matrix.
full(A)

Out[11]:
6×4 Array{Int64,2}:
2  1  0  -1
-1  1  0   0
0  0  2   0
0  1  0   0
-1  0  0   0
0  2  0   0
In [12]:
## Convert sparse A back to coordinate
i,j,v = findnz(A)

Out[12]:
([1, 2, 5, 1, 2, 4, 6, 3, 1], [1, 1, 1, 2, 2, 2, 2, 3, 4], [2, -1, -1, 1, 1, 1, 2, 2, -1])
In [13]:
## Compressed sparse column
I = [1,2,5,1,2,4,6,3,1]
V = [2,-1,-1,1,1,1,2,2,-1]
Jp = [1,4,8,9,10]

Out[13]:
5-element Array{Int64,1}:
1
4
8
9
10
In [14]:
##
A = SparseMatrixCSC(6,4,Jp,I,V)

Out[14]:
6×4 SparseMatrixCSC{Int64,Int64} with 9 stored entries:
[1, 1]  =  2
[2, 1]  =  -1
[5, 1]  =  -1
[1, 2]  =  1
[2, 2]  =  1
[4, 2]  =  1
[6, 2]  =  2
[3, 3]  =  2
[1, 4]  =  -1
In [15]:
##
y = spmatvec_csc(A,x)
yjulia = A*x
y - yjulia

Out[15]:
6×1 Array{Float64,2}:
0.0
0.0
0.0
0.0
0.0
0.0
In [16]:
## julia stores CSC
@show typeof(A)
whos()

typeof(A) = SparseMatrixCSC{Int64,Int64}
A    224 bytes  6×4 SparseMatrixCSC{Int64,Int64}
Base               Module
C    240 bytes  SparseTripletFormat{Int64}
Compat  19368 KB     Module
In [17]:
##
# So this takes 184 + 40 bytes
# it stores: rowval = 64-bit integers, nnz of them
#            nzval = 64-bit doubles, nnz of them
#            colptr = 64-bit integers, n + 1 of them
# 9 non-zeros, n = 4
# 9*8*2 (bytes for I, V) + 5*8 (bytes for Jp)
# sizeof(A) = 40
9*8*2 + 5*8 + 40

Out[17]:
224
In [18]:
## Random sparse matrices
A = sprand(1_000_000,100_000,10/1_000_000)

Out[18]:
1000000×100000 SparseMatrixCSC{Float64,Int64} with 1001420 stored entries:
[67022  ,       1]  =  0.57332
[145433 ,       1]  =  0.0720053
[281961 ,       1]  =  0.921276
[341823 ,       1]  =  0.043949
[363376 ,       1]  =  0.692532
[379345 ,       1]  =  0.632874
[533276 ,       1]  =  0.733147
[820502 ,       1]  =  0.0624762
[832311 ,       1]  =  0.523175
[114761 ,       2]  =  0.826103
⋮
[698317 ,   99999]  =  0.871107
[224514 ,  100000]  =  0.0332364
[259887 ,  100000]  =  0.304153
[466180 ,  100000]  =  0.75642
[690933 ,  100000]  =  0.883819
[719820 ,  100000]  =  0.186653
[767738 ,  100000]  =  0.328502
[792087 ,  100000]  =  0.877643
[821042 ,  100000]  =  0.0590822
[836920 ,  100000]  =  0.134507
[971410 ,  100000]  =  0.480095