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