Skip to content
This repository was archived by the owner on Aug 22, 2025. It is now read-only.

Commit e4cb322

Browse files
A few more array specializations
1 parent 0d6f365 commit e4cb322

File tree

1 file changed

+44
-8
lines changed

1 file changed

+44
-8
lines changed

src/differentiation/compute_jacobian_ad.jl

Lines changed: 44 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -27,11 +27,17 @@ function ForwardColorJacCache(f::F,x,_chunksize = nothing;
2727

2828
if x isa Array
2929
p = generate_chunked_partials(x,colorvec,chunksize)
30+
t = similar(x,Dual{typeof(ForwardDiff.Tag(f,eltype(vec(x))))})
31+
for i in eachindex(t)
32+
t[i] = Dual{typeof(ForwardDiff.Tag(f,eltype(vec(x))))}(x[i],first(p)[1])
33+
end
3034
else
3135
p = adapt.(parameterless_type(x),generate_chunked_partials(x,colorvec,chunksize))
36+
_t = Dual{typeof(ForwardDiff.Tag(f,eltype(vec(x))))}.(vec(x),first(p))
37+
t = ArrayInterface.restructure(x,_t)
3238
end
33-
_t = Dual{typeof(ForwardDiff.Tag(f,eltype(vec(x))))}.(vec(x),first(p))
34-
t = ArrayInterface.restructure(x,_t)
39+
40+
3541
if dx isa Nothing
3642
fx = similar(t)
3743
_dx = similar(x)
@@ -50,16 +56,25 @@ function generate_chunked_partials(x,colorvec,::Val{chunksize}) where chunksize
5056
maxcolor = maximum(colorvec)
5157
num_of_chunks = Int(ceil(maxcolor / chunksize))
5258
padding_size = (chunksize - (maxcolor % chunksize)) % chunksize
53-
partials = colorvec .== (1:maxcolor)'
59+
60+
# partials = colorvec .== (1:maxcolor)'
61+
partials = BitMatrix(undef, length(colorvec), maxcolor)
62+
for j in 1:length(colorvec), i in 1:maxcolor
63+
partials[i,j] = colorvec[j] == i
64+
end
65+
5466
padding_matrix = BitMatrix(undef, length(x), padding_size)
5567
partials = hcat(partials, padding_matrix)
5668

5769
chunked_partials = Vector{Vector{NTuple{chunksize,eltype(x)}}}(undef, num_of_chunks)
5870
for i in 1:num_of_chunks
59-
chunked_partials[i] = Tuple.(eachrow(@view(partials[:,(i-1)*chunksize+1:i*chunksize])))
71+
tmp = Vector{NTuple{chunksize,eltype(x)}}(undef, chunksize)
72+
for j in 1:size(partials,1)
73+
tmp[j] = Tuple(@view partials[j,(i-1)*chunksize+1:i*chunksize])
74+
end
75+
chunked_partials[i] = tmp
6076
end
6177
chunked_partials
62-
6378
end
6479

6580
@inline function forwarddiff_color_jacobian(f,
@@ -287,11 +302,26 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number},
287302

288303
for i in eachindex(p)
289304
partial_i = p[i]
290-
vect .= Dual{typeof(ForwardDiff.Tag(f,eltype(vecx)))}.(vecx, partial_i)
305+
306+
if vect isa Array
307+
@inbounds @simd ivdep for j in eachindex(vect)
308+
vect = Dual{typeof(ForwardDiff.Tag(f,eltype(vecx)))}(vecx[j], partial_i[j])
309+
end
310+
else
311+
vect .= Dual{typeof(ForwardDiff.Tag(f,eltype(vecx)))}.(vecx, partial_i)
312+
end
313+
291314
f(fx,t)
292315
if !(sparsity isa Nothing)
293316
for j in 1:chunksize
294-
dx .= partials.(fx, j)
317+
318+
if dx isa Array
319+
@inbounds @simd ivdep for k in eachindex(dx)
320+
dx[k] = partials(fx[k], j)
321+
end
322+
else
323+
dx .= partials.(fx, j)
324+
end
295325

296326
if ArrayInterface.fast_scalar_indexing(dx)
297327
#dx is implicitly used in vecdx
@@ -320,7 +350,13 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number},
320350
for j in 1:chunksize
321351
col_index = (i-1)*chunksize + j
322352
(col_index > ncols) && return J
323-
J[:, col_index] .= partials.(vecfx, j)
353+
if J isa Array
354+
@inbounds @simd for k in 1:size(J,1)
355+
J[k, col_index] = partials(vecfx[k], j)
356+
end
357+
else
358+
J[:, col_index] .= partials.(vecfx, j)
359+
end
324360
end
325361
end
326362
end

0 commit comments

Comments
 (0)