From 98a1979edc6e716d77b0650a695b3783f816fb28 Mon Sep 17 00:00:00 2001 From: hlw Date: Fri, 26 Jul 2019 21:33:17 +0800 Subject: [PATCH 01/18] add ad support for (tri/bi-)diagonal matrices --- Project.toml | 1 + src/SparseDiffTools.jl | 1 + src/differentiation/compute_jacobian_ad.jl | 6 +++--- test/test_ad.jl | 5 +++++ 4 files changed, 10 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index b9c7add4..f3f4ca73 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" VertexSafeGraphs = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" +ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" [compat] julia = "1" diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index 7e3a8b42..e4442850 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -10,6 +10,7 @@ using VertexSafeGraphs using LinearAlgebra using SparseArrays +using ArrayInterface using BlockBandedMatrices: blocksize, nblocks using ForwardDiff: Dual, jacobian, partials, DEFAULT_CHUNK_THRESHOLD diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index e00d9e53..55fc3a3c 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -21,7 +21,7 @@ getsize(N::Integer) = N function ForwardColorJacCache(f,x,_chunksize = nothing; dx = nothing, color=1:length(x), - sparsity::Union{SparseMatrixCSC,Nothing}=nothing) + sparsity::Union{AbstractArray,Nothing}=nothing) if _chunksize === nothing chunksize = default_chunk_size(maximum(color)) @@ -103,8 +103,8 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, partial_i = p[i] t .= Dual{typeof(f)}.(x, partial_i) f(fx,t) - if sparsity isa SparseMatrixCSC - rows_index, cols_index, val = findnz(sparsity) + if has_sparsestruct(sparsity) + rows_index, cols_index = findstructralnz(sparsity) for j in 1:chunksize dx .= partials.(fx, j) for k in 1:length(cols_index) diff --git a/test/test_ad.jl b/test/test_ad.jl index 3055ce4f..70978535 100644 --- a/test/test_ad.jl +++ b/test/test_ad.jl @@ -1,6 +1,7 @@ using SparseDiffTools using ForwardDiff: Dual, jacobian using SparseArrays, Test +using LinearAlgebra fcalls = 0 function f(dx,x) @@ -56,3 +57,7 @@ jac_cache = ForwardColorJacCache(f,x,color = repeat(1:3,10), sparsity = _J1) forwarddiff_color_jacobian!(_denseJ1, f, x, jac_cache) @test _denseJ1 ≈ J @test fcalls == 1 + +_Jt = similar(Tridiagonal(J)) +forwarddiff_color_jacobian!(_Jt, f, x, color = repeat(1:3,10), sparsity = _Jt) +@test _Jt ≈ J \ No newline at end of file From 804c5593252cb8a8385a7957962c84def693141b Mon Sep 17 00:00:00 2001 From: hlw Date: Fri, 26 Jul 2019 22:01:38 +0800 Subject: [PATCH 02/18] allow diagonal matrix act as sparsity by default --- src/differentiation/compute_jacobian_ad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 55fc3a3c..16c62726 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -81,7 +81,7 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, x::AbstractArray{<:Number}; dx = nothing, color = eachindex(x), - sparsity = J isa SparseMatrixCSC ? J : nothing) + sparsity = has_sparsestruct(J) ? J : nothing) forwarddiff_color_jacobian!(J,f,x,ForwardColorJacCache(f,x,dx=dx,color=color,sparsity=sparsity)) end From 0ab6c06974156423b3156bc58fe0b87e1450de7a Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Mon, 29 Jul 2019 07:34:21 -0400 Subject: [PATCH 03/18] update to ArrayInterface 1.0 --- src/differentiation/compute_jacobian_ad.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 16c62726..60c05295 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -81,7 +81,7 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, x::AbstractArray{<:Number}; dx = nothing, color = eachindex(x), - sparsity = has_sparsestruct(J) ? J : nothing) + sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing) forwarddiff_color_jacobian!(J,f,x,ForwardColorJacCache(f,x,dx=dx,color=color,sparsity=sparsity)) end @@ -103,8 +103,8 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, partial_i = p[i] t .= Dual{typeof(f)}.(x, partial_i) f(fx,t) - if has_sparsestruct(sparsity) - rows_index, cols_index = findstructralnz(sparsity) + if ArrayInterface.has_sparsestruct(sparsity) + rows_index, cols_index = ArrayInterface.findstructralnz(sparsity) for j in 1:chunksize dx .= partials.(fx, j) for k in 1:length(cols_index) From b51b942d2f85e6c6fbf037231171cff2cc59ab48 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Mon, 29 Jul 2019 07:34:45 -0400 Subject: [PATCH 04/18] set lower bound --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index f3f4ca73..88507ed9 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" [compat] julia = "1" +ArrayInterface = "1" [extras] IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" From 8242436da4d026499aea7045d9acb392c42c51a0 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Mon, 29 Jul 2019 07:57:38 -0400 Subject: [PATCH 05/18] update for ArrayInterface 1.0 matrix_colors --- src/coloring/high_level.jl | 76 ++++---------------------------------- 1 file changed, 7 insertions(+), 69 deletions(-) diff --git a/src/coloring/high_level.jl b/src/coloring/high_level.jl index 08978e03..83a16c44 100644 --- a/src/coloring/high_level.jl +++ b/src/coloring/high_level.jl @@ -1,9 +1,9 @@ -abstract type ColoringAlgorithm end -struct GreedyD1Color <: ColoringAlgorithm end -struct BacktrackingColor <: ColoringAlgorithm end -struct ContractionColor <: ColoringAlgorithm end -struct GreedyStar1Color <: ColoringAlgorithm end -struct GreedyStar2Color <: ColoringAlgorithm end +abstract type SparseDiffToolsColoringAlgorithm <: ArrayInterface.ColoringAlgorithm end +struct GreedyD1Color <: SparseDiffToolsColoringAlgorithm end +struct BacktrackingColor <: SparseDiffToolsColoringAlgorithm end +struct ContractionColor <: SparseDiffToolsColoringAlgorithm end +struct GreedyStar1Color <: SparseDiffToolsColoringAlgorithm end +struct GreedyStar2Color <: SparseDiffToolsColoringAlgorithm end """ matrix_colors(A,alg::ColoringAlgorithm = GreedyD1Color()) @@ -13,70 +13,8 @@ struct GreedyStar2Color <: ColoringAlgorithm end The coloring defaults to a greedy distance-1 coloring. """ -function matrix_colors(A::AbstractMatrix,alg::ColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false) +function ArrayInterface.matrix_colors(A::AbstractMatrix,alg::SparseDiffToolsColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false) _A = A isa SparseMatrixCSC ? A : sparse(A) # Avoid the copy A_graph = matrix2graph(_A, partition_by_rows) color_graph(A_graph,alg) end - -""" - matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular}) - - The color vector for dense matrix and triangular matrix is simply - `[1,2,3,...,size(A,2)]` -""" -function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular}) - eachindex(1:size(A,2)) # Vector size matches number of rows -end - -function _cycle(repetend,len) - repeat(repetend,div(len,length(repetend))+1)[1:len] -end - -function matrix_colors(A::Diagonal) - fill(1,size(A,2)) -end - -function matrix_colors(A::Bidiagonal) - _cycle(1:2,size(A,2)) -end - -function matrix_colors(A::Union{Tridiagonal,SymTridiagonal}) - _cycle(1:3,size(A,2)) -end - -function matrix_colors(A::BandedMatrix) - u,l=bandwidths(A) - width=u+l+1 - _cycle(1:width,size(A,2)) -end - -function matrix_colors(A::BlockBandedMatrix) - l,u=blockbandwidths(A) - blockwidth=l+u+1 - nblock=nblocks(A,2) - cols=[blocksize(A,(1,i))[2] for i in 1:nblock] - blockcolors=_cycle(1:blockwidth,nblock) - #the reserved number of colors of a block is the maximum length of columns of blocks with the same block color - ncolors=[maximum(cols[i:blockwidth:nblock]) for i in 1:blockwidth] - endinds=cumsum(ncolors) - startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth] - colors=[(startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]] for i in 1:nblock] - vcat(colors...) -end - -function matrix_colors(A::BandedBlockBandedMatrix) - l,u=blockbandwidths(A) - lambda,mu=subblockbandwidths(A) - blockwidth=l+u+1 - subblockwidth=lambda+mu+1 - nblock=nblocks(A,2) - cols=[blocksize(A,(1,i))[2] for i in 1:nblock] - blockcolors=_cycle(1:blockwidth,nblock) - #the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color - ncolors=[min(subblockwidth,maximum(cols[i:blockwidth:nblock])) for i in 1:min(blockwidth,nblock)] - endinds=cumsum(ncolors) - startinds=[endinds[i]-ncolors[i]+1 for i in 1:min(blockwidth,nblock)] - colors=[_cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]],cols[i]) for i in 1:nblock] - vcat(colors...) -end From fa5de4329b1728bb88604056da2d3ca0afe48b2d Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Mon, 29 Jul 2019 08:00:03 -0400 Subject: [PATCH 06/18] update matrix_coloring docs --- README.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 19c344c3..b4c1a15d 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,9 @@ gmres!(res,J,v) ### Matrix Coloring +This library extends the common `ArrayInterface.matrix_colors` function to allow +for coloring sparse matrices using graphical techniques. + Matrix coloring allows you to reduce the number of times finite differencing requires an `f` call to `maximum(colors)+1`, or reduces automatic differentiation to using `maximum(colors)` partials. Since normally these values are `length(x)`, @@ -101,7 +104,8 @@ this can be significant savings. The API for computing the color vector is: ```julia -matrix_colors(A::AbstractMatrix,alg::ColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false) +matrix_colors(A::AbstractMatrix,alg::ColoringAlgorithm = GreedyD1Color(); + partition_by_rows::Bool = false) ``` The first argument is the abstract matrix which represents the sparsity pattern From 88ccb815e982ed07f3966fce04fe600dd089156d Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 26 Jun 2019 22:32:05 -0400 Subject: [PATCH 07/18] vectorized partial decompression add a comment to the very non-obvious vectorized code --- src/differentiation/compute_jacobian_ad.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 60c05295..569be173 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -98,6 +98,7 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, sparsity = jac_cache.sparsity color_i = 1 chunksize = length(first(first(jac_cache.p))) + fill!(J,zero(eltype(J))) for i in 1:length(p) partial_i = p[i] @@ -107,11 +108,20 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, rows_index, cols_index = ArrayInterface.findstructralnz(sparsity) for j in 1:chunksize dx .= partials.(fx, j) + #= for k in 1:length(cols_index) if color[cols_index[k]] == color_i J[rows_index[k], cols_index[k]] = dx[rows_index[k]] end end + + or + + J[rows_index, cols_index] .+= (color[cols_index] .== color_i) .* dx[rows_index] + + += means requires a zero'd out start + =# + @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index, cols_index) color_i += 1 (color_i > maximum(color)) && continue end From 46812c93d552ddc45f5e0123589bfb3910e83187 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 09:03:19 -0400 Subject: [PATCH 08/18] faster sparse matrix decompression and trait based choice --- Project.toml | 1 + src/SparseDiffTools.jl | 3 +- src/differentiation/compute_jacobian_ad.jl | 33 ++++++++++++++-------- 3 files changed, 23 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 88507ed9..39f100a9 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ VertexSafeGraphs = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" [compat] +ArrayInterface = "1.1" julia = "1" ArrayInterface = "1" diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index e4442850..f62ad4f9 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -9,8 +9,7 @@ using Requires using VertexSafeGraphs using LinearAlgebra -using SparseArrays -using ArrayInterface +using SparseArrays, ArrayInterface using BlockBandedMatrices: blocksize, nblocks using ForwardDiff: Dual, jacobian, partials, DEFAULT_CHUNK_THRESHOLD diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 569be173..f39c22a7 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -108,20 +108,29 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, rows_index, cols_index = ArrayInterface.findstructralnz(sparsity) for j in 1:chunksize dx .= partials.(fx, j) - #= - for k in 1:length(cols_index) - if color[cols_index[k]] == color_i - J[rows_index[k], cols_index[k]] = dx[rows_index[k]] + if ArrayInterface.fast_scalar_indexing(x1) + for i in 1:length(cols_index) + if color[cols_index[i]] == color_i + if J isa SparseMatrixCSC + J.nzval[i] = dx[rows_index[i]] + else + J[rows_index[i],cols_index[i]] = dx[rows_index[i]] + end + end + end + else + #= + J.nzval[rows_index] .+= (color[cols_index] .== color_i) .* dx[rows_index] + or + J[rows_index, cols_index] .+= (color[cols_index] .== color_i) .* dx[rows_index] + += means requires a zero'd out start + =# + if J isa SparseMatrixCSC + @.. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index) + else + @.. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index, cols_index) end end - - or - - J[rows_index, cols_index] .+= (color[cols_index] .== color_i) .* dx[rows_index] - - += means requires a zero'd out start - =# - @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index, cols_index) color_i += 1 (color_i > maximum(color)) && continue end From 4708ce291214c00998221af5108f76540c251976 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 30 Jul 2019 09:07:57 -0400 Subject: [PATCH 09/18] Update src/differentiation/compute_jacobian_ad.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Mathieu Besançon --- src/differentiation/compute_jacobian_ad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index f39c22a7..15ebf380 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -98,7 +98,7 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, sparsity = jac_cache.sparsity color_i = 1 chunksize = length(first(first(jac_cache.p))) - fill!(J,zero(eltype(J))) + fill!(J, zero(eltype(J))) for i in 1:length(p) partial_i = p[i] From b40f408f34cdabdf110f72fff40ad49df282a797 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 30 Jul 2019 09:08:55 -0400 Subject: [PATCH 10/18] Update src/differentiation/compute_jacobian_ad.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Mathieu Besançon --- src/differentiation/compute_jacobian_ad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 15ebf380..5be8c5b8 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -100,7 +100,7 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, chunksize = length(first(first(jac_cache.p))) fill!(J, zero(eltype(J))) - for i in 1:length(p) + for i in eachindex(p) partial_i = p[i] t .= Dual{typeof(f)}.(x, partial_i) f(fx,t) From c98efb3b2f3d90db5e0322a5533a3fd54f8d9cd0 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 09:10:03 -0400 Subject: [PATCH 11/18] don't reuse i --- src/differentiation/compute_jacobian_ad.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 5be8c5b8..c8399eb3 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -109,12 +109,12 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, for j in 1:chunksize dx .= partials.(fx, j) if ArrayInterface.fast_scalar_indexing(x1) - for i in 1:length(cols_index) - if color[cols_index[i]] == color_i + for j in 1:length(cols_index) + if color[cols_index[j]] == color_i if J isa SparseMatrixCSC - J.nzval[i] = dx[rows_index[i]] + J.nzval[j] = dx[rows_index[j]] else - J[rows_index[i],cols_index[i]] = dx[rows_index[i]] + J[rows_index[j],cols_index[j]] = dx[rows_index[j]] end end end From d7ce4ef38043890c5d7e300969bcbf9786860840 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 09:11:39 -0400 Subject: [PATCH 12/18] oops k --- src/differentiation/compute_jacobian_ad.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index c8399eb3..6d0c0782 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -109,12 +109,12 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, for j in 1:chunksize dx .= partials.(fx, j) if ArrayInterface.fast_scalar_indexing(x1) - for j in 1:length(cols_index) - if color[cols_index[j]] == color_i + for k in eachindex(cols_index) + if color[cols_index[k]] == color_i if J isa SparseMatrixCSC - J.nzval[j] = dx[rows_index[j]] + J.nzval[k] = dx[rows_index[k]] else - J[rows_index[j],cols_index[j]] = dx[rows_index[j]] + J[rows_index[k],cols_index[k]] = dx[rows_index[k]] end end end From 3aefa2ebedb2fb683a8603b5720b9fb554c5a3f9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 30 Jul 2019 09:54:31 -0400 Subject: [PATCH 13/18] Update test/test_ad.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Mathieu Besançon --- test/test_ad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ad.jl b/test/test_ad.jl index 70978535..29207476 100644 --- a/test/test_ad.jl +++ b/test/test_ad.jl @@ -60,4 +60,4 @@ forwarddiff_color_jacobian!(_denseJ1, f, x, jac_cache) _Jt = similar(Tridiagonal(J)) forwarddiff_color_jacobian!(_Jt, f, x, color = repeat(1:3,10), sparsity = _Jt) -@test _Jt ≈ J \ No newline at end of file +@test _Jt ≈ J From 8a20d3a0536db2dcdbc6ffe04e2f617292687c8e Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 18:52:13 -0400 Subject: [PATCH 14/18] remove duplicated ArrayInterface compat --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 39f100a9..bebbbe47 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,6 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" [compat] ArrayInterface = "1.1" julia = "1" -ArrayInterface = "1" [extras] IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" From 8eef6a569fa7dd2293533aa9f68c17778baff40b Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 20:20:18 -0400 Subject: [PATCH 15/18] continue to export matrix_colors --- src/SparseDiffTools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index f62ad4f9..fa01f83f 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -20,7 +20,7 @@ export contract_color, greedy_star1_coloring, greedy_star2_coloring, matrix2graph, - matrix_colors, + ArrayInterface.matrix_colors, forwarddiff_color_jacobian!, ForwardColorJacCache, auto_jacvec,auto_jacvec!, From ca2a2385daef5621ad4d66c2617984b7d67e6eb6 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 30 Jul 2019 21:17:53 -0400 Subject: [PATCH 16/18] import matrix_colors --- src/SparseDiffTools.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index fa01f83f..fc1e8a0e 100644 --- a/src/SparseDiffTools.jl +++ b/src/SparseDiffTools.jl @@ -14,13 +14,14 @@ using SparseArrays, ArrayInterface using BlockBandedMatrices: blocksize, nblocks using ForwardDiff: Dual, jacobian, partials, DEFAULT_CHUNK_THRESHOLD +using ArrayInterface: matrix_colors export contract_color, greedy_d1, greedy_star1_coloring, greedy_star2_coloring, matrix2graph, - ArrayInterface.matrix_colors, + matrix_colors, forwarddiff_color_jacobian!, ForwardColorJacCache, auto_jacvec,auto_jacvec!, From 8e3a06e2dfeeb3c189a2158bab92043c09390a49 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 31 Jul 2019 05:46:09 -0400 Subject: [PATCH 17/18] No @.. --- src/differentiation/compute_jacobian_ad.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 6d0c0782..076ba0e3 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -126,9 +126,9 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, += means requires a zero'd out start =# if J isa SparseMatrixCSC - @.. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index) + @. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index) else - @.. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index, cols_index) + @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((color,),cols_index) == color_i) * getindex((dx,),rows_index),rows_index, cols_index) end end color_i += 1 From 70cb8b7a1ae3d85b6f7fddd04104c17d0c000920 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 31 Jul 2019 21:41:08 -0400 Subject: [PATCH 18/18] fix decompression indexing --- src/differentiation/compute_jacobian_ad.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/differentiation/compute_jacobian_ad.jl b/src/differentiation/compute_jacobian_ad.jl index 076ba0e3..468fe7d7 100644 --- a/src/differentiation/compute_jacobian_ad.jl +++ b/src/differentiation/compute_jacobian_ad.jl @@ -108,8 +108,8 @@ function forwarddiff_color_jacobian!(J::AbstractMatrix{<:Number}, rows_index, cols_index = ArrayInterface.findstructralnz(sparsity) for j in 1:chunksize dx .= partials.(fx, j) - if ArrayInterface.fast_scalar_indexing(x1) - for k in eachindex(cols_index) + if ArrayInterface.fast_scalar_indexing(dx) + for k in 1:length(cols_index) if color[cols_index[k]] == color_i if J isa SparseMatrixCSC J.nzval[k] = dx[rows_index[k]]