From 474938128deb1075e6d92be201aa97f1e7e6bf83 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 10 Oct 2022 14:38:53 -0500 Subject: [PATCH 01/14] Add `undefmatrix` (similar to `zeromatrix`) This is potentially much faster since it doesn't have to initialize values. For constructing a 1000x1000 `Matrix{Float64}` I benchmark this to be 10x faster which is very important because in cases like an autoswitching solver, you often want to preallocate a jacobian matrix, but might never use it if the equation isn't stiff. --- .../src/ArrayInterfaceCore.jl | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl b/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl index 905dac535..7903aad06 100644 --- a/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl +++ b/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl @@ -512,6 +512,27 @@ function zeromatrix(u::Array{T}) where {T} fill!(out, false) end +""" + undefmatrix(u::AbstractVector) + +Creates the matrix version of `u` with possibly undefined values. Note that this is unique because +`similar(u,length(u),length(u))` returns a mutable type, so it is not type-matching, +while `fill(zero(eltype(u)),length(u),length(u))` doesn't match the array type, +i.e., you'll get a CPU array from a GPU array. The generic fallback is +`u .* u'`, which works on a surprising number of types, but can be broken +with weird (recursive) broadcast overloads. For higher-order tensors, this +returns the matrix linear operator type which acts on the `vec` of the array. +""" +function undefmatrix(u) + x = safevec(u) + x .* x' +end + +# Reduces compile time burdens +function undefematrix(u::Array{T}) where {T} + out = Matrix{T}(undef, length(u), length(u)) +end + """ restructure(x,y) From 937b222fd336bbc91446b32a9edca1fd63dcd865 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 10 Oct 2022 14:51:53 -0500 Subject: [PATCH 02/14] Update lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl Co-authored-by: Chris Elrod --- lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl b/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl index 7903aad06..5e2887edc 100644 --- a/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl +++ b/lib/ArrayInterfaceCore/src/ArrayInterfaceCore.jl @@ -524,13 +524,7 @@ with weird (recursive) broadcast overloads. For higher-order tensors, this returns the matrix linear operator type which acts on the `vec` of the array. """ function undefmatrix(u) - x = safevec(u) - x .* x' -end - -# Reduces compile time burdens -function undefematrix(u::Array{T}) where {T} - out = Matrix{T}(undef, length(u), length(u)) + similar(u, length(u), length(u)) end """ From 324253a232342ecbc6e6285aec0b9d0717521406 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 10 Oct 2022 20:20:06 -0500 Subject: [PATCH 03/14] add tests --- lib/ArrayInterfaceCore/test/runtests.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index 58dbbfecc..82b23a5a8 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -11,7 +11,22 @@ using Test using Aqua Aqua.test_all(ArrayInterfaceCore) -@test zeromatrix(rand(4,4,4)) == zeros(4*4*4,4*4*4) +@testset "zeromatrix and unsafematrix" begin + for T in (Int, Float32, Float64) + for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrix{T})) + v = vecttype(rand(4)) + um = undefmatrix(v) + @test size(um) == (length(v),length(v)) + @test typeof(um) == mattype + @test zeromatrix(v) == zeros(T,length(v),length(v)) + end + v = rand(T,4,4,4) + um = undefmatrix(v) + @test size(um) == (length(v),length(v)) + @test typeof(um) == Matrix{T} + @test zeromatrix(v) == zeros(T,4*4*4,4*4*4) + end +end @testset "matrix colors" begin @test ArrayInterfaceCore.fast_matrix_colors(1) == false From b8475f4340102041e6473b9379bae14f1b54f444 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 10 Oct 2022 20:48:23 -0500 Subject: [PATCH 04/14] fix test --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index 82b23a5a8..6ab84d6ec 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -13,7 +13,7 @@ Aqua.test_all(ArrayInterfaceCore) @testset "zeromatrix and unsafematrix" begin for T in (Int, Float32, Float64) - for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrix{T})) + for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrixCSC{T, Int})) v = vecttype(rand(4)) um = undefmatrix(v) @test size(um) == (length(v),length(v)) From 5788c1e8d1923a21fac23b0b33bfe4a0fee6956d Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 11 Oct 2022 19:36:41 -0500 Subject: [PATCH 05/14] fix typo --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index 6ab84d6ec..69bb4d780 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -14,7 +14,7 @@ Aqua.test_all(ArrayInterfaceCore) @testset "zeromatrix and unsafematrix" begin for T in (Int, Float32, Float64) for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrixCSC{T, Int})) - v = vecttype(rand(4)) + v = vectype(rand(4)) um = undefmatrix(v) @test size(um) == (length(v),length(v)) @test typeof(um) == mattype From bb01cbdd310723f5b77bdc4bd7c0a1a81182036e Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 11 Oct 2022 19:59:53 -0500 Subject: [PATCH 06/14] add `SArray` and `MArray` implimentations --- .../src/ArrayInterfaceStaticArrays.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl index f5c78b668..b1bd6015c 100644 --- a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl +++ b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl @@ -9,6 +9,14 @@ import ArrayInterfaceStaticArraysCore const CanonicalInt = Union{Int,StaticInt} +function ArrayInterface.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L} + return MMatrix{L, L, T, L*L}(undef) +end +# SArray doesn't have an undef constructor and is going to be small enough that this is fine. +function ArrayInterface.undefmatrix(s::SArray) + v = vec(s) + return v.*v' +end ArrayInterface.known_first(::Type{<:StaticArrays.SOneTo}) = 1 ArrayInterface.known_last(::Type{StaticArrays.SOneTo{N}}) where {N} = N ArrayInterface.known_length(::Type{StaticArrays.SOneTo{N}}) where {N} = N From 0114919b0f0001bfa021902de57fc12d8e7014fb Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 14:14:26 -0400 Subject: [PATCH 07/14] typo. --- .../src/ArrayInterfaceStaticArrays.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl index b1bd6015c..2defb4eaf 100644 --- a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl +++ b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl @@ -9,11 +9,11 @@ import ArrayInterfaceStaticArraysCore const CanonicalInt = Union{Int,StaticInt} -function ArrayInterface.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L} +function ArrayInterfaceCore.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L} return MMatrix{L, L, T, L*L}(undef) end # SArray doesn't have an undef constructor and is going to be small enough that this is fine. -function ArrayInterface.undefmatrix(s::SArray) +function ArrayInterfaceCore.undefmatrix(s::SArray) v = vec(s) return v.*v' end From 18484e6254b124d9646ab78846e09c44844a4c9d Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 14:41:12 -0400 Subject: [PATCH 08/14] export undefmatrix --- src/ArrayInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ArrayInterface.jl b/src/ArrayInterface.jl index 26e173974..043b310bf 100644 --- a/src/ArrayInterface.jl +++ b/src/ArrayInterface.jl @@ -4,7 +4,7 @@ using ArrayInterfaceCore import ArrayInterfaceCore: allowed_getindex, allowed_setindex!, aos_to_soa, buffer, parent_type, fast_matrix_colors, findstructralnz, has_sparsestruct, issingular, isstructured, matrix_colors, restructure, lu_instance, - safevec, zeromatrix, ColoringAlgorithm, fast_scalar_indexing, parameterless_type, + safevec, zeromatrix, undefmatrix, ColoringAlgorithm, fast_scalar_indexing, parameterless_type, ndims_index, ndims_shape, is_splat_index, is_forwarding_wrapper, IndicesInfo, childdims, parentdims, map_tuple_type, flatten_tuples, GetIndex, SetIndex!, defines_strides, stride_preserving_index From 8fc2008da8af28dad9c0f00586fc90befd07d728 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 14:41:40 -0400 Subject: [PATCH 09/14] Update ArrayInterfaceStaticArrays.jl --- .../src/ArrayInterfaceStaticArrays.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl index 2defb4eaf..b1bd6015c 100644 --- a/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl +++ b/lib/ArrayInterfaceStaticArrays/src/ArrayInterfaceStaticArrays.jl @@ -9,11 +9,11 @@ import ArrayInterfaceStaticArraysCore const CanonicalInt = Union{Int,StaticInt} -function ArrayInterfaceCore.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L} +function ArrayInterface.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L} return MMatrix{L, L, T, L*L}(undef) end # SArray doesn't have an undef constructor and is going to be small enough that this is fine. -function ArrayInterfaceCore.undefmatrix(s::SArray) +function ArrayInterface.undefmatrix(s::SArray) v = vec(s) return v.*v' end From b427dfcb15a5e0613a37c8df4bbb893f19b575c7 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 15:32:51 -0400 Subject: [PATCH 10/14] fix test --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index 69bb4d780..c1b41832f 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -14,7 +14,7 @@ Aqua.test_all(ArrayInterfaceCore) @testset "zeromatrix and unsafematrix" begin for T in (Int, Float32, Float64) for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrixCSC{T, Int})) - v = vectype(rand(4)) + v = vectype(rand(3:10, 4)) um = undefmatrix(v) @test size(um) == (length(v),length(v)) @test typeof(um) == mattype From 9cd38b1a76c210b01830fe474a63b8a195eda000 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 15:33:36 -0400 Subject: [PATCH 11/14] add docs --- docs/src/api.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/api.md b/docs/src/api.md index ad43f97bc..2f6347d34 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -39,6 +39,7 @@ ArrayInterfaceCore.promote_eltype ArrayInterfaceCore.restructure ArrayInterfaceCore.safevec ArrayInterfaceCore.zeromatrix +ArrayInterfaceCore.undefmatrix ``` ### Types From 8157475c0830b128c47224157f3ce357a6d8f52f Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 17:06:19 -0400 Subject: [PATCH 12/14] fix test --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index c1b41832f..5ad455406 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -1,5 +1,5 @@ using ArrayInterfaceCore -using ArrayInterfaceCore: zeromatrix +using ArrayInterfaceCore: zeromatrix, zeromatrix import ArrayInterfaceCore: has_sparsestruct, findstructralnz, fast_scalar_indexing, lu_instance, parent_type, zeromatrix, IndicesInfo using Base: setindex From 55d514a91a78f70fba03fef14079bc657c858ff5 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 16:33:22 -0500 Subject: [PATCH 13/14] Update lib/ArrayInterfaceCore/test/runtests.jl Co-authored-by: Christopher Rackauckas --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index 5ad455406..dbd4c20e0 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -1,5 +1,5 @@ using ArrayInterfaceCore -using ArrayInterfaceCore: zeromatrix, zeromatrix +using ArrayInterfaceCore: zeromatrix, undefmatrix import ArrayInterfaceCore: has_sparsestruct, findstructralnz, fast_scalar_indexing, lu_instance, parent_type, zeromatrix, IndicesInfo using Base: setindex From 2b2e40c7504e92ba86c85729631164a52101e7fc Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 12 Oct 2022 18:38:31 -0400 Subject: [PATCH 14/14] actually fix tests --- lib/ArrayInterfaceCore/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ArrayInterfaceCore/test/runtests.jl b/lib/ArrayInterfaceCore/test/runtests.jl index dbd4c20e0..3e172d404 100644 --- a/lib/ArrayInterfaceCore/test/runtests.jl +++ b/lib/ArrayInterfaceCore/test/runtests.jl @@ -14,7 +14,7 @@ Aqua.test_all(ArrayInterfaceCore) @testset "zeromatrix and unsafematrix" begin for T in (Int, Float32, Float64) for (vectype, mattype) in ((Vector{T}, Matrix{T}), (SparseVector{T}, SparseMatrixCSC{T, Int})) - v = vectype(rand(3:10, 4)) + v = vectype(rand(T, 4)) um = undefmatrix(v) @test size(um) == (length(v),length(v)) @test typeof(um) == mattype