From c5a5b87cf77920ff3d6e69a6a28889650deec9cb Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 27 Dec 2019 21:23:06 -0500 Subject: [PATCH 1/8] update BlockBandedMatrix iteration for new release --- src/iteration_utils.jl | 119 +++++++++++++++++++++-------------------- 1 file changed, 60 insertions(+), 59 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index cc06c05..729575e 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -22,84 +22,85 @@ _use_findstructralnz(::SparseMatrixCSC) = false function __init__() @require BlockBandedMatrices="ffab5731-97b5-5995-9138-79e8c1846df0" begin - _use_findstructralnz(::BlockBandedMatrices.BandedBlockBandedMatrix) = false - _use_findstructralnz(::BlockBandedMatrices.BlockBandedMatrix) = false + @require BlockArrays="8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin + _use_findstructralnz(::BlockBandedMatrices.BandedBlockBandedMatrix) = false + _use_findstructralnz(::BlockBandedMatrices.BlockBandedMatrix) = false - @inline function _colorediteration!(Jac::BlockBandedMatrices.BandedBlockBandedMatrix, - sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, - rows_index,cols_index,vfx,colorvec,color_i,ncols) - λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) - rs = BlockBandedMatrices.BlockSizes((BlockBandedMatrices.cumulsizes(Jac,1),)) # column block sizes - cs = BlockBandedMatrices.BlockSizes((BlockBandedMatrices.cumulsizes(Jac,2),)) - b = BlockBandedMatrices.BlockArray(vfx,rs) - c = BlockBandedMatrices.BlockArray(colorvec,cs) - @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) - c_v = c.blocks[J.n[1]] - @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) - V = view(Jac,K,J) - b_v = b.blocks[K.n[1]] - data = BlockBandedMatrices.bandeddata(V) - p = pointer(data) - st = stride(data,2) - m,n = size(V) - @inbounds for j=1:n - if c_v[j] == color_i - @inbounds for k=max(1,j-μ):min(m,j+λ) - unsafe_store!(p, b_v[k], (j-1)*st + μ + k - j + 1) + @inline function _colorediteration!(Jac::BlockBandedMatrices.BandedBlockBandedMatrix, + sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, + rows_index,cols_index,vfx,colorvec,color_i,ncols) + λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) + rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes + cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + b = BlockBandedMatrices.BlockArray(vfx,rs) + c = BlockBandedMatrices.BlockArray(colorvec,cs) + @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) + c_v = c.blocks[J.n[1]] + @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) + V = view(Jac,K,J) + b_v = b.blocks[K.n[1]] + data = BlockBandedMatrices.bandeddata(V) + p = pointer(data) + st = stride(data,2) + m,n = size(V) + @inbounds for j=1:n + if c_v[j] == color_i + @inbounds for k=max(1,j-μ):min(m,j+λ) + unsafe_store!(p, b_v[k], (j-1)*st + μ + k - j + 1) + end end end end end end - end - @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, - sparsity::BlockBandedMatrices.BlockBandedMatrix, - rows_index,cols_index,vfx,colorvec,color_i,ncols) - rs = BlockBandedMatrices.BlockSizes((BlockBandedMatrices.cumulsizes(Jac,1),)) # column block sizes - cs = BlockBandedMatrices.BlockSizes((BlockBandedMatrices.cumulsizes(Jac,2),)) - b = BlockBandedMatrices.BlockArray(vfx,rs) - c = BlockBandedMatrices.BlockArray(colorvec,cs) - @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) - c_v = c.blocks[J.n[1]] - blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) - _,n = BlockBandedMatrices.blocksize(Jac,(blockcolrange[1].n[1],J.n[1])) - @inbounds for j = 1:n - if c_v[j] == color_i - @inbounds for K = blockcolrange - V = view(Jac,K,J) - b_v = b.blocks[K.n[1]] - m = size(V,1) - @inbounds for k = 1:m - V[k,j] = b_v[k] + @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, + sparsity::BlockBandedMatrices.BlockBandedMatrix, + rows_index,cols_index,vfx,colorvec,color_i,ncols) + rs = BlockArrays.axes((BlockBandedMatrices.cumulsizes(Jac,1),)) # column block sizes + cs = BlockArrays.axes((BlockBandedMatrices.cumulsizes(Jac,2),)) + b = BlockBandedMatrices.BlockArray(vfx,rs) + c = BlockBandedMatrices.BlockArray(colorvec,cs) + @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) + c_v = c.blocks[J.n[1]] + blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) + _,n = BlockBandedMatrices.blocksize(Jac,(blockcolrange[1].n[1],J.n[1])) + @inbounds for j = 1:n + if c_v[j] == color_i + @inbounds for K = blockcolrange + V = view(Jac,K,J) + b_v = b.blocks[K.n[1]] + m = size(V,1) + @inbounds for k = 1:m + V[k,j] = b_v[k] + end end end end end end - end - end + end - @require BandedMatrices = "aae01518-5342-5314-be14-df237901396f" begin + @require BandedMatrices = "aae01518-5342-5314-be14-df237901396f" begin - _use_findstructralnz(::BandedMatrices.BandedMatrix) = false + _use_findstructralnz(::BandedMatrices.BandedMatrix) = false - @inline function _colorediteration!(Jac::BandedMatrices.BandedMatrix, - sparsity::BandedMatrices.BandedMatrix, - rows_index,cols_index,vfx,colorvec,color_i,ncols) - nrows = size(Jac,1) - l,u = BandedMatrices.bandwidths(Jac) - #data = BandedMatrices.bandeddata(Jac) - @inbounds for col_index in max(1,1-l):min(ncols,ncols+u) - if colorvec[col_index] == color_i - @inbounds for row_index in max(1,col_index-u):min(nrows,col_index+l) - #data[u+row_index-col_index+1,col_index] = vfx[row_index] - Jac[row_index,col_index]=vfx[row_index] + @inline function _colorediteration!(Jac::BandedMatrices.BandedMatrix, + sparsity::BandedMatrices.BandedMatrix, + rows_index,cols_index,vfx,colorvec,color_i,ncols) + nrows = size(Jac,1) + l,u = BandedMatrices.bandwidths(Jac) + #data = BandedMatrices.bandeddata(Jac) + @inbounds for col_index in max(1,1-l):min(ncols,ncols+u) + if colorvec[col_index] == color_i + @inbounds for row_index in max(1,col_index-u):min(nrows,col_index+l) + #data[u+row_index-col_index+1,col_index] = vfx[row_index] + Jac[row_index,col_index]=vfx[row_index] + end end end end end end - end From b8500ac87a203a58dd7a17eac05a46e3fd899d7a Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 27 Dec 2019 22:00:29 -0500 Subject: [PATCH 2/8] try some translations --- src/iteration_utils.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index 729575e..18363ec 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -30,11 +30,11 @@ function __init__() sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) - rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes - cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + rs = BlockArrays.blockaxes(Jac,1) # column block sizes + cs = BlockArrays.blockaxes(Jac,2) b = BlockBandedMatrices.BlockArray(vfx,rs) c = BlockBandedMatrices.BlockArray(colorvec,cs) - @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) + @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) V = view(Jac,K,J) @@ -57,11 +57,11 @@ function __init__() @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, sparsity::BlockBandedMatrices.BlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) - rs = BlockArrays.axes((BlockBandedMatrices.cumulsizes(Jac,1),)) # column block sizes - cs = BlockArrays.axes((BlockBandedMatrices.cumulsizes(Jac,2),)) + rs = BlockArrays.blockaxes(Jac,1) # column block sizes + cs = BlockArrays.blockaxes(Jac,2) b = BlockBandedMatrices.BlockArray(vfx,rs) c = BlockBandedMatrices.BlockArray(colorvec,cs) - @inbounds for J=BlockBandedMatrices.Block.(1:BlockBandedMatrices.nblocks(Jac,2)) + @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) _,n = BlockBandedMatrices.blocksize(Jac,(blockcolrange[1].n[1],J.n[1])) From aae3222c4e343dbe43af7d839df64c07187f78cc Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 27 Dec 2019 22:05:51 -0500 Subject: [PATCH 3/8] fix a few things --- src/iteration_utils.jl | 4 ++-- test/coloring_tests.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index 18363ec..ca35a3c 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -32,8 +32,8 @@ function __init__() λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) rs = BlockArrays.blockaxes(Jac,1) # column block sizes cs = BlockArrays.blockaxes(Jac,2) - b = BlockBandedMatrices.BlockArray(vfx,rs) - c = BlockBandedMatrices.BlockArray(colorvec,cs) + b = BlockBandedMatrices.BlockArray(vfx,collect(rs)) + c = BlockBandedMatrices.BlockArray(colorvec,collect(cs)) @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) diff --git a/test/coloring_tests.jl b/test/coloring_tests.jl index 0d3a5b2..230ff93 100644 --- a/test/coloring_tests.jl +++ b/test/coloring_tests.jl @@ -106,13 +106,13 @@ function f(out, x) return vec(out) end x = rand(10000) -Jbbb = BandedBlockBandedMatrix(Ones(10000, 10000), (fill(100, 100), fill(100, 100)), (1, 1), (1, 1)) +Jbbb = BandedBlockBandedMatrix(Ones(10000, 10000), fill(100, 100), fill(100, 100), (1, 1), (1, 1)) Jsparse = sparse(Jbbb) colorsbbb = ArrayInterface.matrix_colors(Jbbb) DiffEqDiffTools.finite_difference_jacobian!(Jbbb, f, x, colorvec=colorsbbb) DiffEqDiffTools.finite_difference_jacobian!(Jsparse, f, x, colorvec=colorsbbb) @test Jbbb ≈ Jsparse -Jbb = BlockBandedMatrix(similar(Jsparse),(fill(100, 100), fill(100, 100)),(1,1)); +Jbb = BlockBandedMatrix(similar(Jsparse),fill(100, 100), fill(100, 100),(1,1)); colorsbb = ArrayInterface.matrix_colors(Jbb) DiffEqDiffTools.finite_difference_jacobian!(Jbb, f, x, colorvec=colorsbb) @test Jbb ≈ Jsparse From f8e189d2a6a24da22c332a944fc29d44bd6dc74a Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 27 Dec 2019 22:49:33 -0500 Subject: [PATCH 4/8] try conversion of ranges of Block to ranges of Int --- src/iteration_utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index ca35a3c..4803d50 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -32,8 +32,8 @@ function __init__() λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) rs = BlockArrays.blockaxes(Jac,1) # column block sizes cs = BlockArrays.blockaxes(Jac,2) - b = BlockBandedMatrices.BlockArray(vfx,collect(rs)) - c = BlockBandedMatrices.BlockArray(colorvec,collect(cs)) + b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) + c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) @@ -59,8 +59,8 @@ function __init__() rows_index,cols_index,vfx,colorvec,color_i,ncols) rs = BlockArrays.blockaxes(Jac,1) # column block sizes cs = BlockArrays.blockaxes(Jac,2) - b = BlockBandedMatrices.BlockArray(vfx,rs) - c = BlockBandedMatrices.BlockArray(colorvec,cs) + b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) + c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) From a0bad18b72d1b0ebf219e31924a03d3157d477c6 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Fri, 27 Dec 2019 23:52:30 -0500 Subject: [PATCH 5/8] blocklasts? --- src/iteration_utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index 4803d50..e5058ae 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -30,8 +30,8 @@ function __init__() sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) - rs = BlockArrays.blockaxes(Jac,1) # column block sizes - cs = BlockArrays.blockaxes(Jac,2) + rs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,1)) # column block sizes + cs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,2)) b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) @inbounds for J=BlockArrays.blockaxes(Jac,2) @@ -57,8 +57,8 @@ function __init__() @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, sparsity::BlockBandedMatrices.BlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) - rs = BlockArrays.blockaxes(Jac,1) # column block sizes - cs = BlockArrays.blockaxes(Jac,2) + rs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,1)) # column block sizes + cs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,2)) b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) @inbounds for J=BlockArrays.blockaxes(Jac,2) From f3cdfb953aa9420bb45a80c767b9e8dccb00606b Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Sat, 28 Dec 2019 00:09:26 -0500 Subject: [PATCH 6/8] try some more combinations of blocklasts and axes --- src/iteration_utils.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index e5058ae..e57b250 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -30,10 +30,10 @@ function __init__() sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) - rs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,1)) # column block sizes - cs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,2)) - b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) - c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) + rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes + cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + b = BlockBandedMatrices.BlockArray(vfx,rs) + c = BlockBandedMatrices.BlockArray(colorvec,cs) @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] @inbounds for K=BlockBandedMatrices.blockcolrange(Jac,J) @@ -57,10 +57,10 @@ function __init__() @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, sparsity::BlockBandedMatrices.BlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) - rs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,1)) # column block sizes - cs = BlockArrays.blocklasts.(BlockArrays.blockaxes(Jac,2)) - b = BlockBandedMatrices.BlockArray(vfx,Int.(rs)) - c = BlockBandedMatrices.BlockArray(colorvec,Int.(cs)) + rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes + cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + b = BlockBandedMatrices.BlockArray(vfx,rs) + c = BlockBandedMatrices.BlockArray(colorvec,cs) @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) From 221898fa64821def02fa3f7420e4bcee541aca30 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Sat, 28 Dec 2019 02:38:18 -0500 Subject: [PATCH 7/8] one step closer --- src/iteration_utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index e57b250..6423e88 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -30,8 +30,8 @@ function __init__() sparsity::BlockBandedMatrices.BandedBlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) λ,μ = BlockBandedMatrices.subblockbandwidths(Jac) - rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes - cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + rs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) # column block sizes + cs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) b = BlockBandedMatrices.BlockArray(vfx,rs) c = BlockBandedMatrices.BlockArray(colorvec,cs) @inbounds for J=BlockArrays.blockaxes(Jac,2) @@ -57,8 +57,8 @@ function __init__() @inline function _colorediteration!(Jac::BlockBandedMatrices.BlockBandedMatrix, sparsity::BlockBandedMatrices.BlockBandedMatrix, rows_index,cols_index,vfx,colorvec,color_i,ncols) - rs = BlockArrays.blocklasts(BlockArrays.axes(Jac,1)) # column block sizes - cs = BlockArrays.blocklasts(BlockArrays.axes(Jac,2)) + rs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) # column block sizes + cs = BlockArrays.blocklengths(BlockArrays.axes(Jac,1)) b = BlockBandedMatrices.BlockArray(vfx,rs) c = BlockBandedMatrices.BlockArray(colorvec,cs) @inbounds for J=BlockArrays.blockaxes(Jac,2) From de2efdd0588334186d7f30d7ff53075c3933e375 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Sat, 28 Dec 2019 10:27:02 -0500 Subject: [PATCH 8/8] iteration complete --- src/iteration_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/iteration_utils.jl b/src/iteration_utils.jl index 6423e88..cdcbda7 100644 --- a/src/iteration_utils.jl +++ b/src/iteration_utils.jl @@ -64,7 +64,7 @@ function __init__() @inbounds for J=BlockArrays.blockaxes(Jac,2) c_v = c.blocks[J.n[1]] blockcolrange = BlockBandedMatrices.blockcolrange(Jac,J) - _,n = BlockBandedMatrices.blocksize(Jac,(blockcolrange[1].n[1],J.n[1])) + _,n = length.(getindex.(axes(Jac), (blockcolrange[1], J))) @inbounds for j = 1:n if c_v[j] == color_i @inbounds for K = blockcolrange