Line data Source code
1 : # TODO: We could work directly on the components of SparseMatrixCSC,
2 : # see spmap!
3 :
4 26 : function spmap(f, A::SparseMatrixCSC; dropzeros::Bool=false)
5 13 : m, n = size(A)
6 13 : i, j, v = findnz(A)
7 13 : w = map(f, v)
8 13 : if dropzeros
9 4 : k = findall(isnonzero, w)
10 4 : sparse(i[k], j[k], w[k], m, n)
11 : else
12 9 : sparse(i, j, w, m, n)
13 : end
14 : end
15 :
16 8 : spmap(f, A::SparseT; kwargs...) = spmap(f, A'; kwargs...)'
17 :
18 : """
19 : B = spmap(f, A[; dropzeros=false])
20 :
21 : Get sparse matrix with values of `A` mapped by `f`.
22 :
23 : If `A` is a sparse matrix, `spmap(A)` gives a sparse matrix and
24 : `spmap(A')` gives an *adjoint* of a sparse matrix.
25 :
26 : Unless `dropzeros=true`, the result will include structural zeros
27 : (i.e., all values that mapped to zero).
28 :
29 : """ spmap
30 :
31 :
32 : """
33 :
34 : B = spfilter(predicate, A[; dropzeros=true])
35 :
36 : Get sparse matrix with entries `predicate(A[i,j])`.
37 :
38 : Calls [`spmap`](@ref) but with default setting `dropzeros=true`.
39 :
40 : See also [`spmap`](@ref)
41 :
42 : """
43 12 : function spfilter(predicate, A::SparseOp; dropzeros::Bool=true)
44 : # TODO: alternative: SparseArrays.fkeep! if dropzeros
45 :
46 6 : ZERO = zero(eltype(A))
47 60 : spmap(v->predicate(v) ? v : ZERO, A, dropzeros=dropzeros)
48 : end
49 :
50 : """
51 : S = spones(A[; typ=eltype(A), dropzeros=false])
52 :
53 : Get sparse matrix with nonzero entries replaced by `1`. Calls
54 : [`spmap`](@ref). Note that any structural zeros are *not* replaced by
55 : `1`.
56 :
57 : See also [`spmap`](@ref)
58 : """
59 10 : function spones(A::SparseOp; typ::Type=eltype(A), dropzeros::Bool=false)
60 5 : ONE = one(typ)
61 : # default: keep structural zeros
62 1404359 : spmap(v->iszero(v) ? v : ONE, A, dropzeros=dropzeros)
63 : end
64 :
65 : """
66 : P = sppositive(A[; dropzeros=false])
67 :
68 : Get sparse matrix with positive entries of `A`. Calls
69 : [`spfilter`](@ref).
70 :
71 : See also [`spmap`](@ref)
72 : """
73 12 : sppositive(A::SparseOp; dropzeros::Bool=false) =
74 : spfilter(ispositive, A, dropzeros=dropzeros)
75 :
76 : #-----------------------------------------------------------------------------
77 :
78 4923 : function _mapv!(f, v)
79 4923 : Threads.@threads for k in eachindex(v)
80 345112907 : @inbounds v[k] = f(v[k])
81 0 : end
82 :
83 : # @turbo inline = true unroll = 2 thread = true @. v = f(v)
84 :
85 4923 : v
86 : end
87 :
88 5042 : function _mapiv!(f, i, v)
89 5042 : @assert length(v) == length(i)
90 :
91 5042 : Threads.@threads for k in eachindex(v)
92 353608186 : @inbounds v[k] = f(i[k], v[k])
93 0 : end
94 :
95 : # @turbo inline = true unroll = 2 thread = true @. v = f(i, v)
96 :
97 5042 : v
98 : end
99 :
100 : """
101 : spmap!(f, A::SparseOp)
102 :
103 : Apply `f(v)` to all nonzero values `v` in-place and return `A`.
104 :
105 : See also [`spmapi!`](@ref), [`spmapj!`](@ref), [`spmapij!`](@ref), [`scale!`](@ref)
106 : """
107 4923 : function spmap!(f, A::SparseOp)
108 4923 : _mapv!(f, parent(A).nzval)
109 4923 : A
110 : end
111 :
112 5039 : function spmapi!(f, A::SparseMatrixCSC)
113 5039 : _mapiv!(f, A.rowval, A.nzval)
114 5039 : A
115 : end
116 :
117 3 : function spmapi!(f, A::SparseT)
118 12129 : spmapij!((i, _, v) -> f(i, v), A)
119 3 : A
120 : end
121 :
122 : """
123 : spmapi!(f, A::SparseOp) -> A
124 :
125 : Apply `f(i, v)` to all nonzero values `v` providing the row index `i`
126 : and return `A`.
127 :
128 : !!! note
129 : This operation is more expensive on transposed matrices.
130 :
131 : See also [`spmap!`](@ref), [`spmapj!`](@ref), [`spmapij!`](@ref), [`scale!`](@ref)
132 : """ spmapi!
133 :
134 :
135 3 : function spmapj!(f, A::SparseT)
136 3 : _mapiv!(f, parent(A).rowval, parent(A).nzval)
137 3 : A
138 : end
139 :
140 5493 : function spmapj!(f, A::SparseMatrixCSC)
141 385485405 : spmapij!((_, j, v) -> f(j, v), A)
142 5493 : A
143 : end
144 :
145 : """
146 : spmapj!(f, A::SparseOp) -> A
147 :
148 : Apply `f(j, v)` to all nonzero values `v` providing the column index `j`
149 : and return `A`.
150 :
151 : !!! note
152 : This operation is less expensive on transposed matrices.
153 :
154 : See also [`spmap!`](@ref), [`spmapi!`](@ref), [`spmapij!`](@ref), [`scale!`](@ref)
155 : """ spmapj!
156 :
157 : """
158 : spmapij!(f, A::SparseOp) -> A
159 :
160 : Apply `f(i, j, v)` to all nonzero values `v` providing the row and
161 : column indices `i, j` and return `A`.
162 :
163 : See also [`spmapi!`](@ref), [`spmapj!`](@ref), [`spmapij!`](@ref), [`scale!`](@ref)
164 : """
165 5498 : function spmapij!(f, A::SparseOp)
166 5498 : M = parent(A)
167 5498 : _, n = size(M)
168 5498 : @assert length(M.colptr) == n+1
169 :
170 385485810 : g(::SparseMatrixCSC, i, j, v) = f(i, j, v)
171 16168 : g(::SparseT, i, j, v) = f(j, i, v)
172 :
173 5498 : for j in 1:n, k in M.colptr[j]:(M.colptr[j+1]-1)
174 385486877 : i = M.rowval[k]
175 385486877 : v = M.nzval[k]
176 :
177 385486877 : M.nzval[k] = g(A, i, j, v)
178 390976487 : end
179 :
180 5498 : A
181 : end
182 :
183 4921 : function scale!(A::SparseOp, α)
184 345112591 : spmap!(v -> (α * v), A)
185 4921 : A
186 : end
187 :
188 : """
189 : scalerows!(A::SpaceOp, w[, α=1])
190 :
191 : Scale rows of `A` in-place as for `A = A * (α * Diagonal(w))` and
192 : return `A`.
193 :
194 : See also [`scale!`](@ref), [`scalecols`](@ref)
195 : """
196 5040 : function scalerows!(A::SparseOp, w::AbstractVector, α=eltype(A)(1))
197 5040 : @assert length(w) == size(A, 1)
198 :
199 353607882 : spmapi!((i, v) -> (α * v * w[i]), A)
200 : end
201 :
202 : # #
203 : # function scalerows!(A::SparseMatrixCSC, w::AbstractVector, α=eltype(A)(1))
204 : # @assert length(w) == size(A, 1)
205 :
206 : # W = Vector{eltype(A)}(undef, nnz(A))
207 : # for i in eachindex(A.nzval)
208 : # W[i] = w[A.rowval[i]] * α
209 : # end
210 :
211 : # A.nzval .*= W
212 :
213 : # # for i in eachindex(A.nzval)
214 : # # A.nzval[i] *= α * w[A.rowval[i]]
215 : # # end
216 :
217 : # A
218 : # end
219 : # #
220 :
221 : """
222 : scalecols!(A::SpaceOp, w[, α=1])
223 :
224 : Scale columns of `A` in-place as for `A = α * Diagonal(w) * A` and
225 : return `A`.
226 :
227 : See also [`scale!`](@ref), [`scalerows`](@ref)
228 : """
229 5494 : function scalecols!(A::SparseOp, w::AbstractVector{T}, α=eltype(A)(1)) where {T}
230 5494 : @assert length(w) == size(A, 2)
231 :
232 385485810 : spmapj!((j, v) -> α * v * w[j], A)
233 : end
234 :
235 2516 : scale!(D::Diagonal, A::SparseOp, α=eltype(A)(1)) = scalerows!(A, parent(D), α)
236 :
237 8237 : scale!(A::SparseOp, D::Diagonal, α=eltype(A)(1)) = scalecols!(A, parent(D), α)
238 :
239 2 : scale!(::UniformScaling, A::SparseOp, α) = scale!(A, α)
240 :
241 2 : scale!(A::SparseOp, ::UniformScaling, α) = scale!(A, α)
242 :
243 : """
244 : In-place methods for `A = A * α`
245 :
246 : scale!(A, α)
247 : scale!(I, A, α)
248 : scale!(A, I, α)
249 :
250 : In-place methods for `A = A * (α * Diagonal(w))`, i.e., scale rows
251 :
252 : scalerows!(A, w[, α=1])
253 : scale!(A, Diagonal(w)[, α=1])
254 :
255 : In-place methods for `A = α * Diagonal(w) * A`, i.e., scale rows
256 :
257 : scalecols!(A, w[, α=1])
258 : scale!(Diagonal(w), A, [, α=1])
259 :
260 : Scale rows or columns of sparse operator `A` in-place and
261 : return `A`.
262 :
263 : !!! note
264 : `A::SparseOp` can be a sparse matrix or a transposed matrix.
265 : Mind different cost: scaling rows (or columns of `A'`) is
266 : less expensive than scaling columns (or rows of `A'`).
267 :
268 : !!! note
269 : All methods return `A`, i.e., the object that was modified.
270 :
271 : See also [`scalerows`](@ref), [`scalecols`](@ref), [`spmap!`](@ref),
272 : [`spmapi!`](@ref), [`spmapj!`](@ref)
273 : """ scale!
|
