@@ -13,12 +13,54 @@ struct SparseArray{T,N} <: AbstractArray{T,N}
13
13
new {T,N} (copy (A. data), A. dims)
14
14
end
15
15
end
16
+ SparseArray {T} (:: UndefInitializer , dims... ) where {T} = SparseArray {T} (undef, dims)
17
+
18
+ nonzeros (A:: SparseArray{<:Any,1} ) = (first (k)=> v for (k,v) in A. data)
19
+ nonzeros (A:: SparseArray ) = (CartesianIndex (k)=> v for (k,v) in A. data)
20
+
21
+ function SparseArray (A:: Adjoint{T,<:SparseArray{T,2}} ) where T
22
+ B = SparseArray {T} (undef, size (A))
23
+ for (I, v) in nonzeros (parent (A))
24
+ B[I[2 ], I[1 ]] = v
25
+ end
26
+ return B
27
+ end
16
28
17
29
memsize (A:: SparseArray ) = memsize (A. data)
18
30
@inline function Base. getindex (A:: SparseArray{T,N} , I:: Vararg{Int,N} ) where {T,N}
19
31
@boundscheck checkbounds (A, I... )
20
32
return get (A. data, I, zero (T))
21
33
end
34
+
35
+ _newindex (i:: Int , range:: Int ) = i == range ? () : nothing
36
+ function _newindex (i:: Int , range:: AbstractVector{Int} )
37
+ k = findfirst (== (i), range)
38
+ k === nothing ? nothing : (k,)
39
+ end
40
+ _newindices (I:: Tuple{} , indices:: Tuple{} ) = ()
41
+ function _newindices (I:: Tuple , indices:: Tuple )
42
+ i = _newindex (I[1 ], indices[1 ])
43
+ Itail = _newindices (Base. tail (I), Base. tail (indices))
44
+ (i === nothing || Itail === nothing ) && return nothing
45
+ return (i... , Itail... )
46
+ end
47
+
48
+ _findfirstvalue (v, r) = findfirst (== (v), r)
49
+ # slicing should produce SparseArray
50
+ function Base. _unsafe_getindex (:: IndexCartesian , A:: SparseArray{T,N} ,
51
+ I:: Vararg{<:Union{Int,AbstractVector{Int}},N} ) where {T,N}
52
+ @boundscheck checkbounds (A, I... )
53
+ indices = Base. to_indices (A, I)
54
+ B = SparseArray {T} (undef, length .(Base. index_shape (indices... )))
55
+ for (k, v) in A. data
56
+ newI = _newindices (k, indices)
57
+ if newI != = nothing
58
+ B[newI... ] = v
59
+ end
60
+ end
61
+ return B
62
+ end
63
+
22
64
@inline function Base. setindex! (A:: SparseArray{T,N} , v, I:: Vararg{Int,N} ) where {T,N}
23
65
@boundscheck checkbounds (A, I... )
24
66
if v != zero (v)
@@ -59,32 +101,32 @@ Base.similar(A::SparseArray, ::Type{S}, dims::Dims{N}) where {S,N} =
59
101
60
102
# Vector space functions
61
103
# ------------------------
62
- function LinearAlgebra. lmul! (a:: Number , d :: SparseArray )
63
- lmul! (a, d . data. vals)
104
+ function LinearAlgebra. lmul! (a:: Number , x :: SparseArray )
105
+ lmul! (a, x . data. vals)
64
106
# typical occupation in a dict is about 30% from experimental testing
65
107
# the benefits of scaling all values (e.g. SIMD) largely outweight the extra work
66
- return d
108
+ return x
67
109
end
68
- function LinearAlgebra. rmul! (d :: SparseArray , a:: Number )
69
- rmul! (d . data. vals, a)
70
- return d
110
+ function LinearAlgebra. rmul! (x :: SparseArray , a:: Number )
111
+ rmul! (x . data. vals, a)
112
+ return x
71
113
end
72
114
function LinearAlgebra. axpby! (α, x:: SparseArray , β, y:: SparseArray )
73
115
lmul! (y, β)
74
- for (k, v) in x
116
+ for (k, v) in nonzeros (x)
75
117
y[k] += α* v
76
118
end
77
119
return y
78
120
end
79
121
function LinearAlgebra. axpy! (α, x:: SparseArray , y:: SparseArray )
80
- for (k, v) in x
122
+ for (k, v) in nonzeros (x)
81
123
y[k] += α* v
82
124
end
83
125
return y
84
126
end
85
127
86
128
function LinearAlgebra. norm (x:: SparseArray , p:: Real = 2 )
87
- norm (Base. Generator (last, A . data), p)
129
+ norm (Base. Generator (last, x . data), p)
88
130
end
89
131
90
132
function LinearAlgebra. dot (x:: SparseArray , y:: SparseArray )
0 commit comments