diff --git a/Project.toml b/Project.toml
index a71078c..8136d00 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,39 +1,32 @@
 name = "SparseArraysBase"
 uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
-Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BroadcastMapConversion = "4a4adec5-520f-4750-bb37-d5e66b4ddeb2"
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+Derive = "a07dfc7f-7d04-4eb5-84cc-a97f051f655a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-NestedPermutedDimsArrays = "2c2a8ec4-3cfc-4276-aa3e-1307b4294e58"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
-VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [compat]
-Accessors = "0.1.38"
-Adapt = "4.1.1"
 Aqua = "0.8.9"
-ArrayLayouts = "1.10.4"
-Compat = "4.16.0"
+ArrayLayouts = "1.11.0"
+BroadcastMapConversion = "0.1.0"
+Derive = "0.3.0"
 Dictionaries = "0.4.3"
 LinearAlgebra = "1.10"
-MacroTools = "0.5.13"
-SparseArrays = "1.10"
+SafeTestsets = "0.1"
+Suppressor = "0.2"
 Test = "1.10"
-VectorInterface = "0.5.0"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["Aqua", "Test", "Suppressor", "SafeTestsets"]
diff --git a/README.md b/README.md
index 80a0543..bdf69d5 100644
--- a/README.md
+++ b/README.md
@@ -32,10 +32,127 @@ julia> Pkg.add("SparseArraysBase")
 ## Examples
 
 ````julia
-using SparseArraysBase: SparseArraysBase
+using SparseArraysBase:
+  SparseArrayDOK,
+  SparseMatrixDOK,
+  SparseVectorDOK,
+  eachstoredindex,
+  getstoredindex,
+  getunstoredindex,
+  isstored,
+  setstoredindex!,
+  setunstoredindex!,
+  storedlength,
+  storedpairs,
+  storedvalues,
+  zero!
+using Test: @test, @test_throws
+
+a = SparseArrayDOK{Float64}(2, 2)
+````
+
+AbstractArray interface:
+
+````julia
+@test iszero(a)
+@test iszero(sum(a))
+@test iszero(storedlength(a))
+
+a[1, 2] = 12
+@test a == [0 12; 0 0]
+@test a[1, 1] == 0
+@test a[2, 1] == 0
+@test a[1, 2] == 12
+@test a[2, 2] == 0
+````
+
+SparseArraysBase interface:
+
+````julia
+using Dictionaries: IndexError
+@test issetequal(eachstoredindex(a), [CartesianIndex(1, 2)])
+@test getstoredindex(a, 1, 2) == 12
+@test_throws IndexError getstoredindex(a, 1, 1)
+@test getunstoredindex(a, 1, 1) == 0
+@test getunstoredindex(a, 1, 2) == 0
+@test !isstored(a, 1, 1)
+@test isstored(a, 1, 2)
+@test setstoredindex!(copy(a), 21, 1, 2) == [0 21; 0 0]
+@test_throws IndexError setstoredindex!(copy(a), 21, 2, 1)
+@test setunstoredindex!(copy(a), 21, 1, 2) == [0 21; 0 0]
+@test storedlength(a) == 1
+@test issetequal(storedpairs(a), [CartesianIndex(1, 2) => 12])
+@test issetequal(storedvalues(a), [12])
+@test sum(a) == 12
+@test isreal(a)
+@test !iszero(a)
+@test mapreduce(x -> 2x, +, a) == 24
+````
+
+AbstractArray functionality:
+
+````julia
+b = a .+ 2 .* a'
+@test b isa SparseMatrixDOK{Float64}
+@test b == [0 12; 24 0]
+@test storedlength(b) == 2
+@test sum(b) == 36
+@test isreal(b)
+@test !iszero(b)
+@test mapreduce(x -> 2x, +, b) == 72
+
+b = permutedims(a, (2, 1))
+@test b isa SparseMatrixDOK{Float64}
+@test b[1, 1] == a[1, 1]
+@test b[2, 1] == a[1, 2]
+@test b[1, 2] == a[2, 1]
+@test b[2, 2] == a[2, 2]
+
+b = a * a'
+@test b isa SparseMatrixDOK{Float64}
+@test b == [144 0; 0 0]
+@test storedlength(b) == 1
 ````
 
-Examples go here.
+Second column.
+
+````julia
+b = a[1:2, 2]
+@test b isa SparseVectorDOK{Float64}
+@test b == [12, 0]
+@test storedlength(b) == 1
+
+a = SparseArrayDOK{Float64}(2, 2)
+a .= 2
+for I in eachindex(a)
+  @test a[I] == 2
+end
+@test storedlength(a) == length(a)
+
+a = SparseArrayDOK{Float64}(2, 2)
+fill!(a, 2)
+for I in eachindex(a)
+  @test a[I] == 2
+end
+@test storedlength(a) == length(a)
+
+a = SparseArrayDOK{Float64}(2, 2)
+fill!(a, 0)
+@test iszero(a)
+@test iszero(storedlength(a))
+
+a = SparseArrayDOK{Float64}(2, 2)
+a[1, 2] = 12
+zero!(a)
+@test iszero(a)
+@test iszero(storedlength(a))
+
+a = SparseArrayDOK{Float64}(2, 2)
+a[1, 2] = 12
+b = zero(a)
+@test iszero(b)
+@test iszero(storedlength(b))
+````
 
 ---
 
diff --git a/TODO.md b/TODO.md
index 300cf62..97b0318 100644
--- a/TODO.md
+++ b/TODO.md
@@ -1,124 +1 @@
-# To-do
-
-- Remove Compat.jl from dependencies and test dependencies.
-- Unregistered dependencies: BroadcastMapConversion, NestedPermutedDimsArrays, TypeParameterAccessors.
-- Create SparseArraysBaseLinearAlgebraExt, SparseArraysBaseNestedPermutedDimsArraysExt, SparseArraysBaseVectorInterfaceExt.
-- Change [sources] entries from paths to urls.
-
-# SparseArraysBase
-
-SparseArraysBase is a package that aims to expand on the sparse array functionality that is currently in Julia Base.
-While SparseArrays.jl is centered mostly around `SparseMatrixCSC` and the SuiteSparse library, here we wish to broaden the scope a bit, and consider generic sparse arrays.
-Abstractly, the mental model can be considered as a storage object that holds the stored values, and a bijection between the array indices and the indices of the storage.
-For now, we focus on providing efficient implementations of Dictionary of Key (DOK) type sparse storage formats, but may expand upon this in the future.
-As a result, for typical linear algebra routines, we still expect `SparseMatrixCSC` to be the object of choice.
-
-The design consists of roughly three components:
-- `AbstractSparseArray` interface functions
-- Overloaded Julia base methods
-- `SparseArrayDOK` struct that implements this
-
-## AbstractSparseArray
-
-The first part consists of typical functions that are useful in the context of sparse arrays.
-The minimal interface, which enables the usage of the rest of this package, consists of the following functions:
-
-| Signature | Description | Default |
-|-----------|-------------|---------|
-| `sparse_storage(a::AbstractArray)` | Returns the storage object of the sparse array | `a` |
-| `storage_index_to_index(a::AbstractArray, I)` | Converts a storage index to an array index | `I` |
-| `index_to_storage_index(a::AbstractArray, I)` | Converts an array index to a storage index | `I` |
-
-Using these primitives, several convenience functions are defined to facilitate the writing of sparse array algorithms.
-
-| Signature | Description | Default |
-|-----------|-------------|---------|
-| `storage_indices(a)` | Returns the indices of the storage | `eachindex(sparse_storage(a))` |
-| `stored_indices(a)` | Returns the indices of the stored values | `Iterators.map(Base.Fix1(storage_index_to_index, a), storage_indices(a))` |
-| `stored_length(a)` | Returns the number of stored values | `length(storage_indices(a))` |
-
-<!-- TODO: `getindex!`, `increaseindex!`, `sparse_map`, expose "zero" functionality?  -->
-
-Interesting to note here is that the design is such that we can define sparse arrays without having to subtype `AbstractSparseArray`.
-To achieve this, each function `f` is defined in terms of `sparse_f`, rather than directly overloading `f`.
-<!--
-TODO:
-In order to opt-in to the sparse array functionality, one needs to dispatch the functions through `sparse_f` instead of `f`.
-For convenience, you can automatically dispatch all functions through `sparse_f` by using the following macro:
-
-```julia
-@abstractsparsearray MySparseArrayType
-```
--->
-
-## Overloaded Julia base methods
-
-The second part consists of overloading Julia base methods to work with sparse arrays.
-In particular, specialised implementations exist for the following functions:
-
-- `sparse_similar`
-- `sparse_reduce`
-- `sparse_map`
-- `sparse_map!`
-- `sparse_all`
-- `sparse_any`
-- `sparse_isequal`
-- `sparse_fill!`
-- `sparse_zero`, `sparse_zero!`, `sparse_iszero`
-- `sparse_one`, `sparse_one!`, `sparse_isone`
-- `sparse_reshape`, `sparse_reshape!`
-- `sparse_cat`, `sparse_cat!`
-- `sparse_copy!`, `sparse_copyto!`
-- `sparse_permutedims`, `sparse_permutedims!`
-- `sparse_mul!`, `sparse_dot`
-
-## SparseArrayDOK
-
-Finally, the `SparseArrayDOK` struct is provided as a concrete implementation of the `AbstractSparseArray` interface.
-It is a dictionary of keys (DOK) type sparse array, which stores the values in a `Dictionaries.jl` dictionary, and maps the indices to the keys of the dictionary.
-This model is particularly useful for sparse arrays with a small number of non-zero elements, or for arrays that are constructed incrementally, as it boasts fast random accesses and insertions.
-The drawback is that sequential iteration is slower than for other sparse array types, leading to slower linear algebra operations.
-For the purposes of `SparseArraysBase`, this struct will serve as the canonical example of a sparse array, and will be returned by default when new sparse arrays are created.
-
-One particular feature of `SparseArrayDOK` is that it can be used in cases where the non-stored entries have to be constructed in a non-trivial way.
-Typically, sparse arrays use `zero(eltype(a))` to construct the non-stored entries, but this is not always sufficient.
-A concrete example is found in `BlockSparseArrays.jl`, where initialization of the non-stored entries requires the construction of a block of zeros of appropriate size.
-
-<!-- TODO: update TODOs -->
-
-## TODO
-Still need to implement `Base` functions:
-```julia
-[x] sparse_zero(a::AbstractArray) = similar(a)
-[x] sparse_iszero(a::AbstractArray) = iszero(nonzero_length(a)) # Uses `all`, make `sparse_all`?
-[x] sparse_one(a::AbstractArray) = ...
-[x] sparse_isreal(a::AbstractArray) = ... # Uses `all`, make `sparse_all`?
-[x] sparse_isequal(a1::AbstractArray, a2::AbstractArray) = ...
-[x] sparse_conj!(a::AbstractArray) = conj!(nonzeros(a))
-[x] sparse_reshape(a::AbstractArray, dims) = ...
-[ ] sparse_all(f, a::AbstractArray) = ...
-[ ] sparse_getindex(a::AbstractArray, 1:2, 2:3) = ... # Slicing
-```
-`LinearAlgebra` functions:
-```julia
-[ ] sparse_mul!
-[ ] sparse_lmul!
-[ ] sparse_ldiv!
-[ ] sparse_rdiv!
-[ ] sparse_axpby!
-[ ] sparse_axpy!
-[ ] sparse_norm
-[ ] sparse_dot/sparse_inner
-[ ] sparse_adoint!
-[ ] sparse_transpose!
-
-# Using conversion to `SparseMatrixCSC`:
-[ ] sparse_qr
-[ ] sparse_eigen
-[ ] sparse_svd
-```
-`TensorAlgebra` functions:
-```julia
-[ ] add!
-[ ] contract!
-```
+- Updates for latest Derive.
diff --git a/docs/Project.toml b/docs/Project.toml
index 8e69378..4a57f54 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,5 @@
 [deps]
-BroadcastMapConversion = "4a4adec5-520f-4750-bb37-d5e66b4ddeb2"
+Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-NestedPermutedDimsArrays = "2c2a8ec4-3cfc-4276-aa3e-1307b4294e58"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
-TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
diff --git a/examples/Project.toml b/examples/Project.toml
index 2b28dad..5727855 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-BroadcastMapConversion = "4a4adec5-520f-4750-bb37-d5e66b4ddeb2"
-NestedPermutedDimsArrays = "2c2a8ec4-3cfc-4276-aa3e-1307b4294e58"
+Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
-TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/examples/README.jl b/examples/README.jl
index 9ebae3b..da64f99 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -37,5 +37,114 @@ julia> Pkg.add("SparseArraysBase")
 
 # ## Examples
 
-using SparseArraysBase: SparseArraysBase
-# Examples go here.
+using SparseArraysBase:
+  SparseArrayDOK,
+  SparseMatrixDOK,
+  SparseVectorDOK,
+  eachstoredindex,
+  getstoredindex,
+  getunstoredindex,
+  isstored,
+  setstoredindex!,
+  setunstoredindex!,
+  storedlength,
+  storedpairs,
+  storedvalues,
+  zero!
+using Test: @test, @test_throws
+
+a = SparseArrayDOK{Float64}(2, 2)
+
+# AbstractArray interface:
+
+@test iszero(a)
+@test iszero(sum(a))
+@test iszero(storedlength(a))
+
+a[1, 2] = 12
+@test a == [0 12; 0 0]
+@test a[1, 1] == 0
+@test a[2, 1] == 0
+@test a[1, 2] == 12
+@test a[2, 2] == 0
+
+# SparseArraysBase interface:
+
+using Dictionaries: IndexError
+@test issetequal(eachstoredindex(a), [CartesianIndex(1, 2)])
+@test getstoredindex(a, 1, 2) == 12
+@test_throws IndexError getstoredindex(a, 1, 1)
+@test getunstoredindex(a, 1, 1) == 0
+@test getunstoredindex(a, 1, 2) == 0
+@test !isstored(a, 1, 1)
+@test isstored(a, 1, 2)
+@test setstoredindex!(copy(a), 21, 1, 2) == [0 21; 0 0]
+@test_throws IndexError setstoredindex!(copy(a), 21, 2, 1)
+@test setunstoredindex!(copy(a), 21, 1, 2) == [0 21; 0 0]
+@test storedlength(a) == 1
+@test issetequal(storedpairs(a), [CartesianIndex(1, 2) => 12])
+@test issetequal(storedvalues(a), [12])
+@test sum(a) == 12
+@test isreal(a)
+@test !iszero(a)
+@test mapreduce(x -> 2x, +, a) == 24
+
+# AbstractArray functionality:
+
+b = a .+ 2 .* a'
+@test b isa SparseMatrixDOK{Float64}
+@test b == [0 12; 24 0]
+@test storedlength(b) == 2
+@test sum(b) == 36
+@test isreal(b)
+@test !iszero(b)
+@test mapreduce(x -> 2x, +, b) == 72
+
+b = permutedims(a, (2, 1))
+@test b isa SparseMatrixDOK{Float64}
+@test b[1, 1] == a[1, 1]
+@test b[2, 1] == a[1, 2]
+@test b[1, 2] == a[2, 1]
+@test b[2, 2] == a[2, 2]
+
+b = a * a'
+@test b isa SparseMatrixDOK{Float64}
+@test b == [144 0; 0 0]
+@test storedlength(b) == 1
+
+# Second column.
+b = a[1:2, 2]
+@test b isa SparseVectorDOK{Float64}
+@test b == [12, 0]
+@test storedlength(b) == 1
+
+a = SparseArrayDOK{Float64}(2, 2)
+a .= 2
+for I in eachindex(a)
+  @test a[I] == 2
+end
+@test storedlength(a) == length(a)
+
+a = SparseArrayDOK{Float64}(2, 2)
+fill!(a, 2)
+for I in eachindex(a)
+  @test a[I] == 2
+end
+@test storedlength(a) == length(a)
+
+a = SparseArrayDOK{Float64}(2, 2)
+fill!(a, 0)
+@test iszero(a)
+@test iszero(storedlength(a))
+
+a = SparseArrayDOK{Float64}(2, 2)
+a[1, 2] = 12
+zero!(a)
+@test iszero(a)
+@test iszero(storedlength(a))
+
+a = SparseArrayDOK{Float64}(2, 2)
+a[1, 2] = 12
+b = zero(a)
+@test iszero(b)
+@test iszero(storedlength(b))
diff --git a/src/SparseArraysBase.jl b/src/SparseArraysBase.jl
index 6d264c6..a7231ca 100644
--- a/src/SparseArraysBase.jl
+++ b/src/SparseArraysBase.jl
@@ -1,38 +1,9 @@
 module SparseArraysBase
 
-include("sparsearrayinterface/arraylayouts.jl")
-include("sparsearrayinterface/densearray.jl")
-include("sparsearrayinterface/vectorinterface.jl")
-include("sparsearrayinterface/interface.jl")
-include("sparsearrayinterface/interface_optional.jl")
-include("sparsearrayinterface/indexing.jl")
-include("sparsearrayinterface/base.jl")
-include("sparsearrayinterface/map.jl")
-include("sparsearrayinterface/copyto.jl")
-include("sparsearrayinterface/broadcast.jl")
-include("sparsearrayinterface/conversion.jl")
-include("sparsearrayinterface/wrappers.jl")
-include("sparsearrayinterface/zero.jl")
-include("sparsearrayinterface/cat.jl")
-include("sparsearrayinterface/SparseArraysBaseLinearAlgebraExt.jl")
-include("abstractsparsearray/abstractsparsearray.jl")
-include("abstractsparsearray/abstractsparsematrix.jl")
-include("abstractsparsearray/abstractsparsevector.jl")
-include("abstractsparsearray/wrappedabstractsparsearray.jl")
-include("abstractsparsearray/arraylayouts.jl")
-include("abstractsparsearray/sparsearrayinterface.jl")
-include("abstractsparsearray/base.jl")
-include("abstractsparsearray/broadcast.jl")
-include("abstractsparsearray/map.jl")
-include("abstractsparsearray/baseinterface.jl")
-include("abstractsparsearray/convert.jl")
-include("abstractsparsearray/cat.jl")
-include("abstractsparsearray/SparseArraysBaseSparseArraysExt.jl")
-include("abstractsparsearray/SparseArraysBaseLinearAlgebraExt.jl")
-include("sparsearraydok/defaults.jl")
-include("sparsearraydok/sparsearraydok.jl")
-include("sparsearraydok/sparsematrixdok.jl")
-include("sparsearraydok/sparsevectordok.jl")
-include("sparsearraydok/arraylayouts.jl")
+include("abstractsparsearrayinterface.jl")
+include("sparsearrayinterface.jl")
+include("wrappers.jl")
+include("abstractsparsearray.jl")
+include("sparsearraydok.jl")
 
 end
diff --git a/src/abstractsparsearray.jl b/src/abstractsparsearray.jl
new file mode 100644
index 0000000..b803e82
--- /dev/null
+++ b/src/abstractsparsearray.jl
@@ -0,0 +1,26 @@
+abstract type AbstractSparseArray{T,N} <: AbstractArray{T,N} end
+
+using Derive: @array_aliases
+# Define AbstractSparseVector, AnyAbstractSparseArray, etc.
+@array_aliases AbstractSparseArray
+
+using Derive: Derive
+function Derive.interface(::Type{<:AbstractSparseArray})
+  return SparseArrayInterface()
+end
+
+using Derive: @derive
+
+# TODO: These need to be loaded since `AbstractArrayOps`
+# includes overloads of functions from these modules.
+# Ideally that wouldn't be needed and can be circumvented
+# with `GlobalRef`.
+using ArrayLayouts: ArrayLayouts
+using LinearAlgebra: LinearAlgebra
+
+# Derive `Base.getindex`, `Base.setindex!`, etc.
+# TODO: Define `AbstractMatrixOps` and overload for
+# `AnyAbstractSparseMatrix` and `AnyAbstractSparseVector`,
+# which is where matrix multiplication and factorizations
+# should go.
+@derive AnyAbstractSparseArray AbstractArrayOps
diff --git a/src/abstractsparsearray/SparseArraysBaseLinearAlgebraExt.jl b/src/abstractsparsearray/SparseArraysBaseLinearAlgebraExt.jl
deleted file mode 100644
index b90dfb5..0000000
--- a/src/abstractsparsearray/SparseArraysBaseLinearAlgebraExt.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-using LinearAlgebra: LinearAlgebra
-
-LinearAlgebra.norm(a::AbstractSparseArray, p::Real=2) = sparse_norm(a, p)
-
-# a1 * a2 * α + a_dest * β
-function LinearAlgebra.mul!(
-  a_dest::AbstractMatrix,
-  a1::AbstractSparseMatrix,
-  a2::AbstractSparseMatrix,
-  α::Number=true,
-  β::Number=false,
-)
-  sparse_mul!(a_dest, a1, a2, α, β)
-  return a_dest
-end
diff --git a/src/abstractsparsearray/SparseArraysBaseSparseArraysExt.jl b/src/abstractsparsearray/SparseArraysBaseSparseArraysExt.jl
deleted file mode 100644
index 7a07003..0000000
--- a/src/abstractsparsearray/SparseArraysBaseSparseArraysExt.jl
+++ /dev/null
@@ -1,20 +0,0 @@
-using Base: Forward
-using SparseArrays: SparseArrays, SparseMatrixCSC, findnz, getcolptr, nonzeros, rowvals
-using SparseArraysBase: stored_length
-
-# Julia Base `AbstractSparseArray` interface
-SparseArrays.nnz(a::AbstractSparseArray) = stored_length(a)
-
-sparse_storage(a::SparseMatrixCSC) = nonzeros(a)
-function storage_index_to_index(a::SparseMatrixCSC, I)
-  I1s, I2s = findnz(a)
-  return CartesianIndex(I1s[I], I2s[I])
-end
-function index_to_storage_index(a::SparseMatrixCSC, I::CartesianIndex{2})
-  i0, i1 = Tuple(I)
-  r1 = getcolptr(a)[i1]
-  r2 = getcolptr(a)[i1 + 1] - 1
-  (r1 > r2) && return nothing
-  r1 = searchsortedfirst(rowvals(a), i0, r1, r2, Forward)
-  return ((r1 > r2) || (rowvals(a)[r1] != i0)) ? nothing : r1
-end
diff --git a/src/abstractsparsearray/abstractsparsearray.jl b/src/abstractsparsearray/abstractsparsearray.jl
deleted file mode 100644
index 8c73f7e..0000000
--- a/src/abstractsparsearray/abstractsparsearray.jl
+++ /dev/null
@@ -1,3 +0,0 @@
-using ArrayLayouts: LayoutArray
-
-abstract type AbstractSparseArray{T,N} <: LayoutArray{T,N} end
diff --git a/src/abstractsparsearray/abstractsparsematrix.jl b/src/abstractsparsearray/abstractsparsematrix.jl
deleted file mode 100644
index a545c56..0000000
--- a/src/abstractsparsearray/abstractsparsematrix.jl
+++ /dev/null
@@ -1 +0,0 @@
-const AbstractSparseMatrix{T} = AbstractSparseArray{T,2}
diff --git a/src/abstractsparsearray/abstractsparsevector.jl b/src/abstractsparsearray/abstractsparsevector.jl
deleted file mode 100644
index 033df09..0000000
--- a/src/abstractsparsearray/abstractsparsevector.jl
+++ /dev/null
@@ -1 +0,0 @@
-const AbstractSparseVector{T} = AbstractSparseArray{T,1}
diff --git a/src/abstractsparsearray/arraylayouts.jl b/src/abstractsparsearray/arraylayouts.jl
deleted file mode 100644
index 5b9baf8..0000000
--- a/src/abstractsparsearray/arraylayouts.jl
+++ /dev/null
@@ -1,41 +0,0 @@
-using ArrayLayouts: ArrayLayouts, Dot, DualLayout, MatMulMatAdd, MatMulVecAdd, MulAdd
-using LinearAlgebra: Adjoint, Transpose
-using TypeParameterAccessors: parenttype
-
-function ArrayLayouts.MemoryLayout(arraytype::Type{<:AnyAbstractSparseArray})
-  return SparseLayout()
-end
-
-# TODO: Generalize to `SparseVectorLike`/`AnySparseVector`.
-function ArrayLayouts.MemoryLayout(arraytype::Type{<:Adjoint{<:Any,<:AbstractSparseVector}})
-  return DualLayout{typeof(MemoryLayout(parenttype(arraytype)))}()
-end
-# TODO: Generalize to `SparseVectorLike`/`AnySparseVector`.
-function ArrayLayouts.MemoryLayout(
-  arraytype::Type{<:Transpose{<:Any,<:AbstractSparseVector}}
-)
-  return DualLayout{typeof(MemoryLayout(parenttype(arraytype)))}()
-end
-
-function sparse_matmul!(m::MulAdd)
-  α, a1, a2, β, a_dest = m.α, m.A, m.B, m.β, m.C
-  sparse_mul!(a_dest, a1, a2, α, β)
-  return a_dest
-end
-
-function ArrayLayouts.materialize!(
-  m::MatMulMatAdd{<:AbstractSparseLayout,<:AbstractSparseLayout,<:AbstractSparseLayout}
-)
-  sparse_matmul!(m)
-  return m.C
-end
-function ArrayLayouts.materialize!(
-  m::MatMulVecAdd{<:AbstractSparseLayout,<:AbstractSparseLayout,<:AbstractSparseLayout}
-)
-  sparse_matmul!(m)
-  return m.C
-end
-
-function Base.copy(d::Dot{<:SparseLayout,<:SparseLayout})
-  return sparse_dot(d.A, d.B)
-end
diff --git a/src/abstractsparsearray/base.jl b/src/abstractsparsearray/base.jl
deleted file mode 100644
index 0ea8265..0000000
--- a/src/abstractsparsearray/base.jl
+++ /dev/null
@@ -1,16 +0,0 @@
-# Base
-function Base.:(==)(a1::AnyAbstractSparseArray, a2::AnyAbstractSparseArray)
-  return sparse_isequal(a1, a2)
-end
-
-function Base.reshape(a::AnyAbstractSparseArray, dims::Tuple{Vararg{Int}})
-  return sparse_reshape(a, dims)
-end
-
-function Base.zero(a::AnyAbstractSparseArray)
-  return sparse_zero(a)
-end
-
-function Base.one(a::AnyAbstractSparseArray)
-  return sparse_one(a)
-end
diff --git a/src/abstractsparsearray/baseinterface.jl b/src/abstractsparsearray/baseinterface.jl
deleted file mode 100644
index ba16da7..0000000
--- a/src/abstractsparsearray/baseinterface.jl
+++ /dev/null
@@ -1,33 +0,0 @@
-Base.size(a::AbstractSparseArray) = error("Not implemented")
-
-function Base.similar(a::AbstractSparseArray, elt::Type, dims::Tuple{Vararg{Int}})
-  return error("Not implemented")
-end
-
-function Base.getindex(a::AbstractSparseArray, I...)
-  return sparse_getindex(a, I...)
-end
-
-# Fixes ambiguity error with `ArrayLayouts`.
-function Base.getindex(a::AbstractSparseMatrix, I1::AbstractVector, I2::AbstractVector)
-  return sparse_getindex(a, I1, I2)
-end
-
-# Fixes ambiguity error with `ArrayLayouts`.
-function Base.getindex(
-  a::AbstractSparseMatrix, I1::AbstractUnitRange, I2::AbstractUnitRange
-)
-  return sparse_getindex(a, I1, I2)
-end
-
-function Base.isassigned(a::AbstractSparseArray, I::Integer...)
-  return sparse_isassigned(a, I...)
-end
-
-function Base.setindex!(a::AbstractSparseArray, I...)
-  return sparse_setindex!(a, I...)
-end
-
-function Base.fill!(a::AbstractSparseArray, value)
-  return sparse_fill!(a, value)
-end
diff --git a/src/abstractsparsearray/broadcast.jl b/src/abstractsparsearray/broadcast.jl
deleted file mode 100644
index 565fccb..0000000
--- a/src/abstractsparsearray/broadcast.jl
+++ /dev/null
@@ -1,4 +0,0 @@
-# Broadcasting
-function Broadcast.BroadcastStyle(arraytype::Type{<:AnyAbstractSparseArray})
-  return SparseArraysBase.SparseArrayStyle{ndims(arraytype)}()
-end
diff --git a/src/abstractsparsearray/cat.jl b/src/abstractsparsearray/cat.jl
deleted file mode 100644
index 3d04751..0000000
--- a/src/abstractsparsearray/cat.jl
+++ /dev/null
@@ -1,4 +0,0 @@
-# TODO: Change to `AnyAbstractSparseArray`.
-function Base.cat(as::AnyAbstractSparseArray...; dims)
-  return sparse_cat(as...; dims)
-end
diff --git a/src/abstractsparsearray/convert.jl b/src/abstractsparsearray/convert.jl
deleted file mode 100644
index 8b532ed..0000000
--- a/src/abstractsparsearray/convert.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-Base.convert(type::Type{<:AbstractSparseArray}, a::AbstractArray) = type(a)
-
-Base.convert(::Type{T}, a::T) where {T<:AbstractSparseArray} = a
-
-function (::Type{T})(a::T) where {T<:AbstractSparseArray}
-  return copy(a)
-end
diff --git a/src/abstractsparsearray/map.jl b/src/abstractsparsearray/map.jl
deleted file mode 100644
index 106cfff..0000000
--- a/src/abstractsparsearray/map.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-using ArrayLayouts: LayoutArray
-
-# Map
-function Base.map!(f, a_dest::AbstractArray, a_srcs::Vararg{AnyAbstractSparseArray})
-  SparseArraysBase.sparse_map!(f, a_dest, a_srcs...)
-  return a_dest
-end
-
-function Base.copy!(a_dest::AbstractArray, a_src::AnyAbstractSparseArray)
-  SparseArraysBase.sparse_copy!(a_dest, a_src)
-  return a_dest
-end
-
-function Base.copyto!(a_dest::AbstractArray, a_src::AnyAbstractSparseArray)
-  SparseArraysBase.sparse_copyto!(a_dest, a_src)
-  return a_dest
-end
-
-# Fix ambiguity error
-function Base.copyto!(a_dest::LayoutArray, a_src::AnyAbstractSparseArray)
-  SparseArraysBase.sparse_copyto!(a_dest, a_src)
-  return a_dest
-end
-
-function Base.permutedims!(a_dest::AbstractArray, a_src::AnyAbstractSparseArray, perm)
-  SparseArraysBase.sparse_permutedims!(a_dest, a_src, perm)
-  return a_dest
-end
-
-function Base.mapreduce(f, op, as::Vararg{AnyAbstractSparseArray}; kwargs...)
-  return SparseArraysBase.sparse_mapreduce(f, op, as...; kwargs...)
-end
-
-# TODO: Why isn't this calling `mapreduce` already?
-function Base.iszero(a::AnyAbstractSparseArray)
-  return SparseArraysBase.sparse_iszero(a)
-end
-
-# TODO: Why isn't this calling `mapreduce` already?
-function Base.isreal(a::AnyAbstractSparseArray)
-  return SparseArraysBase.sparse_isreal(a)
-end
diff --git a/src/abstractsparsearray/sparsearrayinterface.jl b/src/abstractsparsearray/sparsearrayinterface.jl
deleted file mode 100644
index 8b86d1b..0000000
--- a/src/abstractsparsearray/sparsearrayinterface.jl
+++ /dev/null
@@ -1,46 +0,0 @@
-using Dictionaries: set!
-
-sparse_storage(::AbstractSparseArray) = error("Not implemented")
-
-function index_to_storage_index(
-  a::AbstractSparseArray{<:Any,N}, I::CartesianIndex{N}
-) where {N}
-  !isassigned(sparse_storage(a), I) && return nothing
-  return I
-end
-
-function setindex_notstored!(
-  a::AbstractSparseArray{<:Any,N}, value, I::CartesianIndex{N}
-) where {N}
-  iszero(value) && return a
-  return error("Setting the specified unstored index is not supported.")
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-# TODO: Check if this is efficient, or determine if this mapping should
-# be performed in `storage_index_to_index` and/or `index_to_storage_index`.
-function sparse_storage(a::SubArray{<:Any,<:Any,<:AbstractSparseArray})
-  parent_storage = sparse_storage(parent(a))
-  all_sliced_storage_indices = map(keys(parent_storage)) do I
-    return map_index(a.indices, I)
-  end
-  sliced_storage_indices = filter(!isnothing, all_sliced_storage_indices)
-  sliced_parent_storage = map(I -> parent_storage[I], keys(sliced_storage_indices))
-  return typeof(parent_storage)(sliced_storage_indices, sliced_parent_storage)
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-function stored_indices(
-  a::AnyPermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:AbstractSparseArray}
-)
-  return Iterators.map(
-    I -> CartesianIndex(map(i -> I[i], perm(a))), stored_indices(parent(a))
-  )
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-function sparse_storage(
-  a::AnyPermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:AbstractSparseArray}
-)
-  return sparse_storage(parent(a))
-end
diff --git a/src/abstractsparsearray/wrappedabstractsparsearray.jl b/src/abstractsparsearray/wrappedabstractsparsearray.jl
deleted file mode 100644
index a4ee4be..0000000
--- a/src/abstractsparsearray/wrappedabstractsparsearray.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-using Adapt: WrappedArray
-using LinearAlgebra: Adjoint, Transpose
-
-const WrappedAbstractSparseArray{T,N,A} = WrappedArray{
-  T,N,<:AbstractSparseArray,<:AbstractSparseArray{T,N}
-}
-
-const AnyAbstractSparseArray{T,N} = Union{
-  <:AbstractSparseArray{T,N},<:WrappedAbstractSparseArray{T,N}
-}
-
-function stored_indices(a::Adjoint)
-  return Iterators.map(I -> CartesianIndex(reverse(Tuple(I))), stored_indices(parent(a)))
-end
-stored_length(a::Adjoint) = stored_length(parent(a))
-sparse_storage(a::Adjoint) = Iterators.map(adjoint, sparse_storage(parent(a)))
-
-function stored_indices(a::Transpose)
-  return Iterators.map(I -> CartesianIndex(reverse(Tuple(I))), stored_indices(parent(a)))
-end
-stored_length(a::Transpose) = stored_length(parent(a))
-sparse_storage(a::Transpose) = Iterators.map(transpose, sparse_storage(parent(a)))
diff --git a/src/abstractsparsearrayinterface.jl b/src/abstractsparsearrayinterface.jl
new file mode 100644
index 0000000..d3d10d0
--- /dev/null
+++ b/src/abstractsparsearrayinterface.jl
@@ -0,0 +1,311 @@
+using Derive: Derive, @derive, @interface, AbstractArrayInterface
+
+# This is to bring `ArrayLayouts.zero!` into the namespace
+# since it is considered part of the sparse array interface.
+using ArrayLayouts: zero!
+
+function eachstoredindex end
+function getstoredindex end
+function getunstoredindex end
+function isstored end
+function setstoredindex! end
+function setunstoredindex! end
+function storedlength end
+function storedpairs end
+function storedvalues end
+
+# Generic functionality for converting to a
+# dense array, trying to preserve information
+# about the array (such as which device it is on).
+# TODO: Maybe call `densecopy`?
+# TODO: Make sure this actually preserves the device,
+# maybe use `TypeParameterAccessors.unwrap_array_type`.
+# TODO: Turn into an `@interface` function.
+function densearray(a::AbstractArray)
+  # TODO: `set_ndims(unwrap_array_type(a), ndims(a))(a)`
+  # Maybe define `densetype(a) = set_ndims(unwrap_array_type(a), ndims(a))`.
+  # Or could use `unspecify_parameters(unwrap_array_type(a))(a)`.
+  return Array(a)
+end
+
+# Minimal interface for `SparseArrayInterface`.
+# Fallbacks for dense/non-sparse arrays.
+@interface ::AbstractArrayInterface isstored(a::AbstractArray, I::Int...) = true
+@interface ::AbstractArrayInterface eachstoredindex(a::AbstractArray) = eachindex(a)
+@interface ::AbstractArrayInterface getstoredindex(a::AbstractArray, I::Int...) =
+  getindex(a, I...)
+@interface ::AbstractArrayInterface function setstoredindex!(
+  a::AbstractArray, value, I::Int...
+)
+  setindex!(a, value, I...)
+  return a
+end
+# TODO: Should this error by default if the value at the index
+# is stored? It could be disabled with something analogous
+# to `checkbounds`, like `checkstored`/`checkunstored`.
+@interface ::AbstractArrayInterface function setunstoredindex!(
+  a::AbstractArray, value, I::Int...
+)
+  # TODO: Make this a `MethodError`?
+  return error("Not implemented.")
+end
+
+# TODO: Use `Base.to_indices`?
+isstored(a::AbstractArray, I::CartesianIndex) = isstored(a, Tuple(I)...)
+# TODO: Use `Base.to_indices`?
+getstoredindex(a::AbstractArray, I::CartesianIndex) = getstoredindex(a, Tuple(I)...)
+# TODO: Use `Base.to_indices`?
+getunstoredindex(a::AbstractArray, I::CartesianIndex) = getunstoredindex(a, Tuple(I)...)
+# TODO: Use `Base.to_indices`?
+function setstoredindex!(a::AbstractArray, value, I::CartesianIndex)
+  return setstoredindex!(a, value, Tuple(I)...)
+end
+# TODO: Use `Base.to_indices`?
+function setunstoredindex!(a::AbstractArray, value, I::CartesianIndex)
+  return setunstoredindex!(a, value, Tuple(I)...)
+end
+
+# Interface defaults.
+# TODO: Have a fallback that handles element types
+# that don't define `zero(::Type)`.
+@interface ::AbstractArrayInterface getunstoredindex(a::AbstractArray, I::Int...) =
+  zero(eltype(a))
+
+# Derived interface.
+@interface ::AbstractArrayInterface storedlength(a::AbstractArray) = length(storedvalues(a))
+@interface ::AbstractArrayInterface storedpairs(a::AbstractArray) =
+  map(I -> I => getstoredindex(a, I), eachstoredindex(a))
+
+@interface ::AbstractArrayInterface function eachstoredindex(as::AbstractArray...)
+  return eachindex(as...)
+end
+
+@interface ::AbstractArrayInterface storedvalues(a::AbstractArray) = a
+
+# Automatically derive the interface for all `AbstractArray` subtypes.
+# TODO: Define `SparseArrayInterfaceOps` derivable trait and rewrite this
+# as `@derive AbstractArray SparseArrayInterfaceOps`.
+@derive (T=AbstractArray,) begin
+  SparseArraysBase.eachstoredindex(::T)
+  SparseArraysBase.eachstoredindex(::T...)
+  SparseArraysBase.getstoredindex(::T, ::Int...)
+  SparseArraysBase.getunstoredindex(::T, ::Int...)
+  SparseArraysBase.isstored(::T, ::Int...)
+  SparseArraysBase.setstoredindex!(::T, ::Any, ::Int...)
+  SparseArraysBase.setunstoredindex!(::T, ::Any, ::Int...)
+  SparseArraysBase.storedlength(::T)
+  SparseArraysBase.storedpairs(::T)
+  SparseArraysBase.storedvalues(::T)
+end
+
+# TODO: Add `ndims` type parameter, like `Base.Broadcast.AbstractArrayStyle`.
+# TODO: This isn't used to define interface functions right now.
+# Currently, `@interface` expects an instance, probably it should take a
+# type instead so fallback functions can use abstract types.
+abstract type AbstractSparseArrayInterface <: AbstractArrayInterface end
+
+to_vec(x) = vec(collect(x))
+to_vec(x::AbstractArray) = vec(x)
+
+# A view of the stored values of an array.
+# Similar to: `@view a[collect(eachstoredindex(a))]`, but the issue
+# with that is it returns a `SubArray` wrapping a sparse array, which
+# is then interpreted as a sparse array so it can lead to recursion.
+# Also, that involves extra logic for determining if the indices are
+# stored or not, but we know the indices are stored so we can use
+# `getstoredindex` and `setstoredindex!`.
+# Most sparse arrays should overload `storedvalues` directly
+# and avoid this wrapper since it adds extra indirection to
+# access stored values.
+struct StoredValues{T,A<:AbstractArray{T},I} <: AbstractVector{T}
+  array::A
+  storedindices::I
+end
+StoredValues(a::AbstractArray) = StoredValues(a, to_vec(eachstoredindex(a)))
+Base.size(a::StoredValues) = size(a.storedindices)
+Base.getindex(a::StoredValues, I::Int) = getstoredindex(a.array, a.storedindices[I])
+function Base.setindex!(a::StoredValues, value, I::Int)
+  return setstoredindex!(a.array, value, a.storedindices[I])
+end
+
+@interface ::AbstractSparseArrayInterface storedvalues(a::AbstractArray) = StoredValues(a)
+
+@interface ::AbstractSparseArrayInterface function eachstoredindex(
+  a1::AbstractArray, a2::AbstractArray, a_rest::AbstractArray...
+)
+  # TODO: Make this more customizable, say with a function
+  # `combine/promote_storedindices(a1, a2)`.
+  return union(eachstoredindex.((a1, a2, a_rest...))...)
+end
+
+@interface ::AbstractSparseArrayInterface function eachstoredindex(a::AbstractArray)
+  # TODO: Use `MethodError`?
+  return error("Not implemented.")
+end
+
+# We restrict to `I::Vararg{Int,N}` to allow more general functions to handle trailing
+# indices and linear indices.
+@interface ::AbstractSparseArrayInterface function Base.getindex(
+  a::AbstractArray{<:Any,N}, I::Vararg{Int,N}
+) where {N}
+  !isstored(a, I...) && return getunstoredindex(a, I...)
+  return getstoredindex(a, I...)
+end
+
+# We restrict to `I::Vararg{Int,N}` to allow more general functions to handle trailing
+# indices and linear indices.
+@interface ::AbstractSparseArrayInterface function Base.setindex!(
+  a::AbstractArray{<:Any,N}, value, I::Vararg{Int,N}
+) where {N}
+  if !isstored(a, I...)
+    # Don't set the value if it is zero, but only check
+    # if it is zero if the elements are numbers since otherwise
+    # it may be nontrivial to check.
+    eltype(a) <: Number && iszero(value) && return a
+    setunstoredindex!(a, value, I...)
+    return a
+  end
+  setstoredindex!(a, value, I...)
+  return a
+end
+
+# TODO: This may need to be defined in `sparsearraydok.jl`, after `SparseArrayDOK`
+# is defined. And/or define `default_type(::SparseArrayStyle, T::Type) = SparseArrayDOK{T}`.
+@interface ::AbstractSparseArrayInterface function Base.similar(
+  a::AbstractArray, T::Type, size::Tuple{Vararg{Int}}
+)
+  # TODO: Define `default_similartype` or something like that?
+  return SparseArrayDOK{T}(size...)
+end
+
+@interface ::AbstractSparseArrayInterface function Base.map!(
+  f, a_dest::AbstractArray, as::AbstractArray...
+)
+  # TODO: Define a function `preserves_unstored(a_dest, f, as...)`
+  # to determine if a function preserves the stored values
+  # of the destination sparse array.
+  # The current code may be inefficient since it actually
+  # accesses an unstored element, which in the case of a
+  # sparse array of arrays can allocate an array.
+  # Sparse arrays could be expected to define a cheap
+  # unstored element allocator, for example
+  # `get_prototypical_unstored(a::AbstractArray)`.
+  I = first(eachindex(as...))
+  preserves_unstored = iszero(f(map(a -> getunstoredindex(a, I), as)...))
+  if !preserves_unstored
+    # Doesn't preserve unstored values, loop over all elements.
+    for I in eachindex(as...)
+      a_dest[I] = map(f, map(a -> a[I], as)...)
+    end
+    return a_dest
+  end
+  # Define `eachstoredindex` promotion.
+  for I in eachstoredindex(as...)
+    a_dest[I] = f(map(a -> a[I], as)...)
+  end
+  return a_dest
+end
+
+# `f::typeof(norm)`, `op::typeof(max)` used by `norm`.
+function reduce_init(f, op, as...)
+  # TODO: Generalize this.
+  @assert isone(length(as))
+  a = only(as)
+  ## TODO: Make this more efficient for block sparse
+  ## arrays, in that case it allocates a block. Maybe
+  ## it can use `FillArrays.Zeros`.
+  return f(getunstoredindex(a, first(eachindex(a))))
+end
+
+@interface ::AbstractSparseArrayInterface function Base.mapreduce(
+  f, op, as::AbstractArray...; init=reduce_init(f, op, as...), kwargs...
+)
+  # TODO: Generalize this.
+  @assert isone(length(as))
+  a = only(as)
+  output = mapreduce(f, op, storedvalues(a); init, kwargs...)
+  ## TODO: Bring this check back, or make the function more general.
+  ## f_notstored = apply_notstored(f, a)
+  ## @assert isequal(op(output, eltype(output)(f_notstored)), output)
+  return output
+end
+
+abstract type AbstractSparseArrayStyle{N} <: Broadcast.AbstractArrayStyle{N} end
+
+@derive AbstractSparseArrayStyle AbstractArrayStyleOps
+
+struct SparseArrayStyle{N} <: AbstractSparseArrayStyle{N} end
+
+SparseArrayStyle{M}(::Val{N}) where {M,N} = SparseArrayStyle{N}()
+
+@interface ::AbstractSparseArrayInterface function Broadcast.BroadcastStyle(type::Type)
+  return SparseArrayStyle{ndims(type)}()
+end
+
+using ArrayLayouts: ArrayLayouts, MatMulMatAdd
+
+abstract type AbstractSparseLayout <: ArrayLayouts.MemoryLayout end
+
+function ArrayLayouts.sub_materialize(::AbstractSparseLayout, a::AbstractArray, axes::Tuple)
+  a_dest = similar(a)
+  a_dest .= a
+  return a_dest
+end
+
+function mul_indices(I1::CartesianIndex{2}, I2::CartesianIndex{2})
+  if I1[2] ≠ I2[1]
+    return nothing
+  end
+  return CartesianIndex(I1[1], I2[2])
+end
+
+using LinearAlgebra: mul!
+function default_mul!!(
+  a_dest::AbstractMatrix,
+  a1::AbstractMatrix,
+  a2::AbstractMatrix,
+  α::Number=true,
+  β::Number=false,
+)
+  mul!(a_dest, a1, a2, α, β)
+  return a_dest
+end
+
+function default_mul!!(
+  a_dest::Number, a1::Number, a2::Number, α::Number=true, β::Number=false
+)
+  return a1 * a2 * α + a_dest * β
+end
+
+# a1 * a2 * α + a_dest * β
+function sparse_mul!(
+  a_dest::AbstractArray,
+  a1::AbstractArray,
+  a2::AbstractArray,
+  α::Number=true,
+  β::Number=false;
+  (mul!!)=(default_mul!!),
+)
+  for I1 in eachstoredindex(a1)
+    for I2 in eachstoredindex(a2)
+      I_dest = mul_indices(I1, I2)
+      if !isnothing(I_dest)
+        a_dest[I_dest] = mul!!(a_dest[I_dest], a1[I1], a2[I2], α, β)
+      end
+    end
+  end
+  return a_dest
+end
+
+function ArrayLayouts.materialize!(
+  m::MatMulMatAdd{<:AbstractSparseLayout,<:AbstractSparseLayout,<:AbstractSparseLayout}
+)
+  sparse_mul!(m.C, m.A, m.B, m.α, m.β)
+  return m.C
+end
+
+struct SparseLayout <: AbstractSparseLayout end
+
+@interface ::AbstractSparseArrayInterface function ArrayLayouts.MemoryLayout(type::Type)
+  return SparseLayout()
+end
diff --git a/src/sparsearraydok.jl b/src/sparsearraydok.jl
new file mode 100644
index 0000000..9249008
--- /dev/null
+++ b/src/sparsearraydok.jl
@@ -0,0 +1,70 @@
+using Dictionaries: Dictionary, IndexError, set!
+
+function default_getunstoredindex(a::AbstractArray, I::Int...)
+  return zero(eltype(a))
+end
+
+struct SparseArrayDOK{T,N,F} <: AbstractSparseArray{T,N}
+  storage::Dictionary{CartesianIndex{N},T}
+  size::NTuple{N,Int}
+  getunstoredindex::F
+end
+
+using Derive: Derive
+# This defines the destination type of various operations in Derive.jl.
+Derive.arraytype(::AbstractSparseArrayInterface, T::Type) = SparseArrayDOK{T}
+
+function SparseArrayDOK{T,N}(size::Vararg{Int,N}) where {T,N}
+  getunstoredindex = default_getunstoredindex
+  F = typeof(getunstoredindex)
+  return SparseArrayDOK{T,N,F}(Dictionary{CartesianIndex{N},T}(), size, getunstoredindex)
+end
+
+function SparseArrayDOK{T}(::UndefInitializer, size::Tuple{Vararg{Int}}) where {T}
+  return SparseArrayDOK{T,length(size)}(size...)
+end
+
+function SparseArrayDOK{T}(size::Int...) where {T}
+  return SparseArrayDOK{T,length(size)}(size...)
+end
+
+using Derive: @array_aliases
+# Define `SparseMatrixDOK`, `AnySparseArrayDOK`, etc.
+@array_aliases SparseArrayDOK
+
+storage(a::SparseArrayDOK) = a.storage
+Base.size(a::SparseArrayDOK) = a.size
+
+storedvalues(a::SparseArrayDOK) = values(storage(a))
+function isstored(a::SparseArrayDOK, I::Int...)
+  return CartesianIndex(I) in keys(storage(a))
+end
+function eachstoredindex(a::SparseArrayDOK)
+  return keys(storage(a))
+end
+function getstoredindex(a::SparseArrayDOK, I::Int...)
+  return storage(a)[CartesianIndex(I)]
+end
+function getunstoredindex(a::SparseArrayDOK, I::Int...)
+  return a.getunstoredindex(a, I...)
+end
+function setstoredindex!(a::SparseArrayDOK, value, I::Int...)
+  # TODO: Have a way to disable this check, analogous to `checkbounds`,
+  # since this is already checked in `setindex!`.
+  isstored(a, I...) || throw(IndexError("key $(CartesianIndex(I)) not found"))
+  # TODO: If `iszero(value)`, unstore the index.
+  storage(a)[CartesianIndex(I)] = value
+  return a
+end
+function setunstoredindex!(a::SparseArrayDOK, value, I::Int...)
+  set!(storage(a), CartesianIndex(I), value)
+  return a
+end
+
+# Optional, but faster than the default.
+storedpairs(a::SparseArrayDOK) = pairs(storage(a))
+
+function ArrayLayouts.zero!(a::SparseArrayDOK)
+  empty!(storage(a))
+  return a
+end
diff --git a/src/sparsearraydok/arraylayouts.jl b/src/sparsearraydok/arraylayouts.jl
deleted file mode 100644
index fc2740d..0000000
--- a/src/sparsearraydok/arraylayouts.jl
+++ /dev/null
@@ -1,13 +0,0 @@
-using ArrayLayouts: ArrayLayouts, MemoryLayout, MulAdd
-
-ArrayLayouts.MemoryLayout(::Type{<:SparseArrayDOK}) = SparseLayout()
-
-# Default sparse array type for `AbstractSparseLayout`.
-default_sparsearraytype(elt::Type) = SparseArrayDOK{elt}
-
-# TODO: Preserve GPU memory! Implement `CuSparseArrayLayout`, `MtlSparseLayout`?
-function Base.similar(
-  ::MulAdd{<:AbstractSparseLayout,<:AbstractSparseLayout}, elt::Type, axes
-)
-  return similar(default_sparsearraytype(elt), axes)
-end
diff --git a/src/sparsearraydok/defaults.jl b/src/sparsearraydok/defaults.jl
deleted file mode 100644
index 1dc674b..0000000
--- a/src/sparsearraydok/defaults.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-using Dictionaries: Dictionary
-
-default_zero() = Zero()
-default_data(type::Type, ndims::Int) = Dictionary{default_keytype(ndims),type}()
-default_keytype(ndims::Int) = CartesianIndex{ndims}
diff --git a/src/sparsearraydok/sparsearraydok.jl b/src/sparsearraydok/sparsearraydok.jl
deleted file mode 100644
index 95a5b14..0000000
--- a/src/sparsearraydok/sparsearraydok.jl
+++ /dev/null
@@ -1,146 +0,0 @@
-using Accessors: @set
-using Dictionaries: Dictionary, set!
-using MacroTools: @capture
-
-# TODO: Parametrize by `data`?
-struct SparseArrayDOK{T,N,Zero} <: AbstractSparseArray{T,N}
-  data::Dictionary{CartesianIndex{N},T}
-  dims::Ref{NTuple{N,Int}}
-  zero::Zero
-  function SparseArrayDOK{T,N,Zero}(data, dims::NTuple{N,Int}, zero) where {T,N,Zero}
-    return new{T,N,Zero}(data, Ref(dims), zero)
-  end
-end
-
-# Constructors
-function SparseArrayDOK(data, dims::Tuple{Vararg{Int}}, zero)
-  return SparseArrayDOK{eltype(data),length(dims),typeof(zero)}(data, dims, zero)
-end
-
-function SparseArrayDOK{T,N,Zero}(dims::Tuple{Vararg{Int}}, zero) where {T,N,Zero}
-  return SparseArrayDOK{T,N,Zero}(default_data(T, N), dims, zero)
-end
-
-function SparseArrayDOK{T,N}(dims::Tuple{Vararg{Int}}, zero) where {T,N}
-  return SparseArrayDOK{T,N,typeof(zero)}(dims, zero)
-end
-
-function SparseArrayDOK{T,N}(dims::Tuple{Vararg{Int}}) where {T,N}
-  return SparseArrayDOK{T,N}(dims, default_zero())
-end
-
-function SparseArrayDOK{T}(dims::Tuple{Vararg{Int}}) where {T}
-  return SparseArrayDOK{T,length(dims)}(dims)
-end
-
-function SparseArrayDOK{T}(dims::Int...) where {T}
-  return SparseArrayDOK{T}(dims)
-end
-
-# Specify zero function
-function SparseArrayDOK{T}(dims::Tuple{Vararg{Int}}, zero) where {T}
-  return SparseArrayDOK{T,length(dims)}(dims, zero)
-end
-
-# undef
-function SparseArrayDOK{T,N,Zero}(
-  ::UndefInitializer, dims::Tuple{Vararg{Int}}, zero
-) where {T,N,Zero}
-  return SparseArrayDOK{T,N,Zero}(dims, zero)
-end
-
-function SparseArrayDOK{T,N}(::UndefInitializer, dims::Tuple{Vararg{Int}}, zero) where {T,N}
-  return SparseArrayDOK{T,N}(dims, zero)
-end
-
-function SparseArrayDOK{T,N}(::UndefInitializer, dims::Tuple{Vararg{Int}}) where {T,N}
-  return SparseArrayDOK{T,N}(dims)
-end
-
-function SparseArrayDOK{T}(::UndefInitializer, dims::Tuple{Vararg{Int}}) where {T}
-  return SparseArrayDOK{T}(dims)
-end
-
-# Axes version
-function SparseArrayDOK{T}(
-  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange}}
-) where {T}
-  @assert all(isone, first.(axes))
-  return SparseArrayDOK{T}(length.(axes))
-end
-
-function SparseArrayDOK{T}(::UndefInitializer, dims::Int...) where {T}
-  return SparseArrayDOK{T}(dims...)
-end
-
-function SparseArrayDOK{T}(::UndefInitializer, dims::Tuple{Vararg{Int}}, zero) where {T}
-  return SparseArrayDOK{T}(dims, zero)
-end
-
-# Base `AbstractArray` interface
-Base.size(a::SparseArrayDOK) = a.dims[]
-
-getindex_zero_function(a::SparseArrayDOK) = a.zero
-function set_getindex_zero_function(a::SparseArrayDOK, f)
-  return @set a.zero = f
-end
-
-function setindex_notstored!(
-  a::SparseArrayDOK{<:Any,N}, value, I::CartesianIndex{N}
-) where {N}
-  set!(sparse_storage(a), I, value)
-  return a
-end
-
-function Base.similar(a::SparseArrayDOK, elt::Type, dims::Tuple{Vararg{Int}})
-  return SparseArrayDOK{elt}(undef, dims, getindex_zero_function(a))
-end
-
-# `SparseArraysBase` interface
-sparse_storage(a::SparseArrayDOK) = a.data
-
-function dropall!(a::SparseArrayDOK)
-  return empty!(sparse_storage(a))
-end
-
-SparseArrayDOK(a::AbstractArray) = SparseArrayDOK{eltype(a)}(a)
-
-SparseArrayDOK{T}(a::AbstractArray) where {T} = SparseArrayDOK{T,ndims(a)}(a)
-
-function SparseArrayDOK{T,N}(a::AbstractArray) where {T,N}
-  return sparse_convert(SparseArrayDOK{T,N}, a)
-end
-
-function Base.resize!(a::SparseArrayDOK{<:Any,N}, new_size::NTuple{N,Integer}) where {N}
-  a.dims[] = new_size
-  return a
-end
-
-function setindex_maybe_grow!(a::SparseArrayDOK{<:Any,N}, value, I::Vararg{Int,N}) where {N}
-  if any(I .> size(a))
-    resize!(a, max.(I, size(a)))
-  end
-  a[I...] = value
-  return a
-end
-
-function is_setindex!_expr(expr::Expr)
-  return is_assignment_expr(expr) && is_getindex_expr(first(expr.args))
-end
-is_setindex!_expr(x) = false
-
-is_getindex_expr(expr::Expr) = (expr.head === :ref)
-is_getindex_expr(x) = false
-
-is_assignment_expr(expr::Expr) = (expr.head === :(=))
-is_assignment_expr(expr) = false
-
-macro maybe_grow(expr)
-  if !is_setindex!_expr(expr)
-    error(
-      "@maybe_grow must be used with setindex! syntax (as @maybe_grow a[i,j,...] = value)"
-    )
-  end
-  @capture(expr, array_[indices__] = value_)
-  return :(setindex_maybe_grow!($(esc(array)), $(esc(value)), $(esc.(indices)...)))
-end
diff --git a/src/sparsearraydok/sparsematrixdok.jl b/src/sparsearraydok/sparsematrixdok.jl
deleted file mode 100644
index f568760..0000000
--- a/src/sparsearraydok/sparsematrixdok.jl
+++ /dev/null
@@ -1 +0,0 @@
-const SparseMatrixDOK{T} = SparseArrayDOK{T,2}
diff --git a/src/sparsearraydok/sparsevectordok.jl b/src/sparsearraydok/sparsevectordok.jl
deleted file mode 100644
index 7cc7df0..0000000
--- a/src/sparsearraydok/sparsevectordok.jl
+++ /dev/null
@@ -1 +0,0 @@
-const SparseVectorDOK{T} = SparseArrayDOK{T,1}
diff --git a/src/sparsearrayinterface.jl b/src/sparsearrayinterface.jl
new file mode 100644
index 0000000..577ff33
--- /dev/null
+++ b/src/sparsearrayinterface.jl
@@ -0,0 +1,11 @@
+using Derive: Derive
+
+struct SparseArrayInterface <: AbstractSparseArrayInterface end
+
+# Convenient shorthand to refer to the sparse interface.
+# Can turn a function into a sparse function with the syntax `sparse(f)`,
+# i.e. `sparse(map)(x -> 2x, randn(2, 2))` while use the sparse
+# version of `map`.
+# const sparse = SparseArrayInterface()
+
+Derive.interface(::Type{<:AbstractSparseArrayStyle}) = SparseArrayInterface()
diff --git a/src/sparsearrayinterface/SparseArraysBaseLinearAlgebraExt.jl b/src/sparsearrayinterface/SparseArraysBaseLinearAlgebraExt.jl
deleted file mode 100644
index 4e4ed25..0000000
--- a/src/sparsearrayinterface/SparseArraysBaseLinearAlgebraExt.jl
+++ /dev/null
@@ -1,78 +0,0 @@
-using LinearAlgebra: dot, mul!, norm
-
-sparse_norm(a::AbstractArray, p::Real=2) = norm(sparse_storage(a))
-
-function mul_indices(I1::CartesianIndex{2}, I2::CartesianIndex{2})
-  if I1[2] ≠ I2[1]
-    return nothing
-  end
-  return CartesianIndex(I1[1], I2[2])
-end
-
-# TODO: Is this needed? Maybe when multiplying vectors?
-function mul_indices(I1::CartesianIndex{1}, I2::CartesianIndex{1})
-  if I1 ≠ I2
-    return nothing
-  end
-  return CartesianIndex(I1)
-end
-
-function default_mul!!(
-  a_dest::AbstractMatrix,
-  a1::AbstractMatrix,
-  a2::AbstractMatrix,
-  α::Number=true,
-  β::Number=false,
-)
-  mul!(a_dest, a1, a2, α, β)
-  return a_dest
-end
-
-function default_mul!!(
-  a_dest::Number, a1::Number, a2::Number, α::Number=true, β::Number=false
-)
-  return a1 * a2 * α + a_dest * β
-end
-
-# a1 * a2 * α + a_dest * β
-function sparse_mul!(
-  a_dest::AbstractArray,
-  a1::AbstractArray,
-  a2::AbstractArray,
-  α::Number=true,
-  β::Number=false;
-  (mul!!)=(default_mul!!),
-)
-  for I1 in stored_indices(a1)
-    for I2 in stored_indices(a2)
-      I_dest = mul_indices(I1, I2)
-      if !isnothing(I_dest)
-        a_dest[I_dest] = mul!!(a_dest[I_dest], a1[I1], a2[I2], α, β)
-      end
-    end
-  end
-  return a_dest
-end
-
-function sparse_dot(a1::AbstractArray, a2::AbstractArray)
-  # This requires that `a1` and `a2` have the same shape.
-  # TODO: Generalize (Base supports dot products of
-  # arrays with the same length but different sizes).
-  size(a1) == size(a2) ||
-    throw(DimensionMismatch("Sizes $(size(a1)) and $(size(a2)) don't match."))
-  dot_dest = zero(Base.promote_op(dot, eltype(a1), eltype(a2)))
-  # TODO: First check if the number of stored elements (`stored_length`, to be renamed
-  # `stored_length`) is smaller in `a1` or `a2` and use whicheven one is smallar
-  # as the outer loop.
-  for I1 in stored_indices(a1)
-    # TODO: Overload and use `Base.isstored(a, I) = I in stored_indices(a)` instead.
-    # TODO: This assumes fast lookup of indices, which may not always be the case.
-    # It could be better to loop over `stored_indices(a2)` and check that
-    # `I1 == I2` instead (say using `mul_indices(I1, I2)`. We could have a trait
-    # `HasFastIsStored(a::AbstractArray)` to choose between the two.
-    if I1 in stored_indices(a2)
-      dot_dest += dot(a1[I1], a2[I1])
-    end
-  end
-  return dot_dest
-end
diff --git a/src/sparsearrayinterface/arraylayouts.jl b/src/sparsearrayinterface/arraylayouts.jl
deleted file mode 100644
index 4ff6cc1..0000000
--- a/src/sparsearrayinterface/arraylayouts.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-using ArrayLayouts: MemoryLayout
-
-abstract type AbstractSparseLayout <: MemoryLayout end
-
-struct SparseLayout <: AbstractSparseLayout end
diff --git a/src/sparsearrayinterface/base.jl b/src/sparsearrayinterface/base.jl
deleted file mode 100644
index 9a6fd24..0000000
--- a/src/sparsearrayinterface/base.jl
+++ /dev/null
@@ -1,126 +0,0 @@
-# This is used when a sparse output structure not matching
-# the input structure is needed, for example when reshaping
-# a DiagonalArray. Overload:
-#
-# sparse_similar(a::AbstractArray, elt::Type, dims::Tuple{Vararg{Int}})
-#
-# as needed.
-function sparse_similar(a::AbstractArray, elt::Type)
-  return similar(a, elt, size(a))
-end
-
-function sparse_similar(a::AbstractArray, dims::Tuple{Vararg{Int}})
-  return sparse_similar(a, eltype(a), dims)
-end
-
-function sparse_similar(a::AbstractArray)
-  return sparse_similar(a, eltype(a), size(a))
-end
-
-function sparse_reduce(op, a::AbstractArray; kwargs...)
-  return sparse_mapreduce(identity, op, a; kwargs...)
-end
-
-function sparse_all(a::AbstractArray)
-  return sparse_reduce(&, a; init=true)
-end
-
-function sparse_all(f, a::AbstractArray)
-  return sparse_mapreduce(f, &, a; init=true)
-end
-
-function sparse_iszero(a::AbstractArray)
-  return sparse_all(iszero, a)
-end
-
-function sparse_isreal(a::AbstractArray)
-  return sparse_all(isreal, a)
-end
-
-# This is equivalent to:
-#
-# sparse_map!(Returns(x), a, a)
-#
-# but we don't use that here since `sparse_fill!`
-# is used inside of `sparse_map!`.
-function sparse_fill!(a::AbstractArray, x)
-  if iszero(x)
-    # This checks that `x` is compatible
-    # with `eltype(a)`.
-    x = convert(eltype(a), x)
-    sparse_zero!(a)
-    return a
-  end
-  for I in eachindex(a)
-    a[I] = x
-  end
-  return a
-end
-
-# This could just call `sparse_fill!`
-# but it avoids a zero construction and check.
-function sparse_zero!(a::AbstractArray)
-  dropall!(a)
-  sparse_zerovector!(a)
-  return a
-end
-
-function sparse_zero(a::AbstractArray)
-  # Need to overload `similar` for custom types
-  a = similar(a)
-  sparse_zerovector!(a)
-  return a
-end
-
-# TODO: Is this a good definition?
-function sparse_zero(arraytype::Type{<:AbstractArray}, dims::Tuple{Vararg{Int}})
-  a = arraytype(undef, dims)
-  sparse_zerovector!(a)
-  return a
-end
-
-function sparse_one!(a::AbstractMatrix)
-  sparse_zero!(a)
-  m, n = size(a)
-  @assert m == n
-  for i in 1:m
-    a[i, i] = one(eltype(a))
-  end
-  return a
-end
-
-function sparse_one(a::AbstractMatrix)
-  a = sparse_zero(a)
-  sparse_one!(a)
-  return a
-end
-
-# TODO: Use `sparse_mapreduce(==, &, a1, a2)`?
-function sparse_isequal(a1::AbstractArray, a2::AbstractArray)
-  Is = collect(stored_indices(a1))
-  intersect!(Is, stored_indices(a2))
-  if !(length(Is) == stored_length(a1) == stored_length(a2))
-    return false
-  end
-  for I in Is
-    a1[I] == a2[I] || return false
-  end
-  return true
-end
-
-function sparse_reshape!(a_dest::AbstractArray, a_src::AbstractArray, dims)
-  @assert length(a_src) == prod(dims)
-  sparse_zero!(a_dest)
-  linear_inds = LinearIndices(a_src)
-  dest_cartesian_inds = CartesianIndices(dims)
-  for I in stored_indices(a_src)
-    a_dest[dest_cartesian_inds[linear_inds[I]]] = a_src[I]
-  end
-  return a_dest
-end
-
-function sparse_reshape(a::AbstractArray, dims)
-  a_reshaped = sparse_similar(a, dims)
-  sparse_reshape!(a_reshaped, a, dims)
-  return a_reshaped
-end
diff --git a/src/sparsearrayinterface/broadcast.jl b/src/sparsearrayinterface/broadcast.jl
deleted file mode 100644
index cb18aa2..0000000
--- a/src/sparsearrayinterface/broadcast.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-using Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted
-using BroadcastMapConversion: map_function, map_args
-
-struct SparseArrayStyle{N} <: AbstractArrayStyle{N} end
-
-# Define for new sparse array types.
-# function Broadcast.BroadcastStyle(arraytype::Type{<:MySparseArray})
-#   return SparseArrayStyle{ndims(arraytype)}()
-# end
-
-SparseArrayStyle(::Val{N}) where {N} = SparseArrayStyle{N}()
-SparseArrayStyle{M}(::Val{N}) where {M,N} = SparseArrayStyle{N}()
-
-Broadcast.BroadcastStyle(a::SparseArrayStyle, ::DefaultArrayStyle{0}) = a
-function Broadcast.BroadcastStyle(::SparseArrayStyle{N}, a::DefaultArrayStyle) where {N}
-  return BroadcastStyle(DefaultArrayStyle{N}(), a)
-end
-function Broadcast.BroadcastStyle(::SparseArrayStyle{N}, ::Broadcast.Style{Tuple}) where {N}
-  return DefaultArrayStyle{N}()
-end
-
-# TODO: Use `allocate_output`, share logic with `map`.
-function Base.similar(bc::Broadcasted{<:SparseArrayStyle}, elt::Type)
-  # TODO: Is this a good definition? Probably should check that
-  # they have consistent axes.
-  return similar(first(map_args(bc)), elt)
-end
-
-# Broadcasting implementation
-function Base.copyto!(
-  dest::AbstractArray{<:Any,N}, bc::Broadcasted{SparseArrayStyle{N}}
-) where {N}
-  # convert to map
-  # flatten and only keep the AbstractArray arguments
-  sparse_map!(map_function(bc), dest, map_args(bc)...)
-  return dest
-end
diff --git a/src/sparsearrayinterface/cat.jl b/src/sparsearrayinterface/cat.jl
deleted file mode 100644
index 9f2b317..0000000
--- a/src/sparsearrayinterface/cat.jl
+++ /dev/null
@@ -1,64 +0,0 @@
-unval(x) = x
-unval(::Val{x}) where {x} = x
-
-# TODO: Assert that `a1` and `a2` start at one.
-axis_cat(a1::AbstractUnitRange, a2::AbstractUnitRange) = Base.OneTo(length(a1) + length(a2))
-function axis_cat(
-  a1::AbstractUnitRange, a2::AbstractUnitRange, a_rest::AbstractUnitRange...
-)
-  return axis_cat(axis_cat(a1, a2), a_rest...)
-end
-function cat_axes(as::AbstractArray...; dims)
-  return ntuple(length(first(axes.(as)))) do dim
-    return if dim in unval(dims)
-      axis_cat(map(axes -> axes[dim], axes.(as))...)
-    else
-      axes(first(as))[dim]
-    end
-  end
-end
-
-function allocate_cat_output(as::AbstractArray...; dims)
-  eltype_dest = promote_type(eltype.(as)...)
-  axes_dest = cat_axes(as...; dims)
-  # TODO: Promote the block types of the inputs rather than using
-  # just the first input.
-  # TODO: Make this customizable with `cat_similar`.
-  # TODO: Base the zero element constructor on those of the inputs,
-  # for example block sparse arrays.
-  return similar(first(as), eltype_dest, axes_dest...)
-end
-
-# https://github.com/JuliaLang/julia/blob/v1.11.1/base/abstractarray.jl#L1748-L1857
-# https://docs.julialang.org/en/v1/base/arrays/#Concatenation-and-permutation
-# This is very similar to the `Base.cat` implementation but handles zero values better.
-function cat_offset!(
-  a_dest::AbstractArray, offsets, a1::AbstractArray, a_rest::AbstractArray...; dims
-)
-  inds = ntuple(ndims(a_dest)) do dim
-    dim in unval(dims) ? offsets[dim] .+ axes(a1, dim) : axes(a_dest, dim)
-  end
-  a_dest[inds...] = a1
-  new_offsets = ntuple(ndims(a_dest)) do dim
-    dim in unval(dims) ? offsets[dim] + size(a1, dim) : offsets[dim]
-  end
-  cat_offset!(a_dest, new_offsets, a_rest...; dims)
-  return a_dest
-end
-function cat_offset!(a_dest::AbstractArray, offsets; dims)
-  return a_dest
-end
-
-# TODO: Define a generic `cat!` function.
-function sparse_cat!(a_dest::AbstractArray, as::AbstractArray...; dims)
-  offsets = ntuple(zero, ndims(a_dest))
-  # TODO: Fill `a_dest` with zeros if needed.
-  cat_offset!(a_dest, offsets, as...; dims)
-  return a_dest
-end
-
-function sparse_cat(as::AbstractArray...; dims)
-  a_dest = allocate_cat_output(as...; dims)
-  sparse_cat!(a_dest, as...; dims)
-  return a_dest
-end
diff --git a/src/sparsearrayinterface/conversion.jl b/src/sparsearrayinterface/conversion.jl
deleted file mode 100644
index 57f1850..0000000
--- a/src/sparsearrayinterface/conversion.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-function sparse_convert(arraytype::Type{<:AbstractArray}, a::AbstractArray)
-  a_dest = sparse_zero(arraytype, size(a))
-  sparse_copyto!(a_dest, a)
-  return a_dest
-end
diff --git a/src/sparsearrayinterface/copyto.jl b/src/sparsearrayinterface/copyto.jl
deleted file mode 100644
index 218502f..0000000
--- a/src/sparsearrayinterface/copyto.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-function sparse_copy!(dest::AbstractArray, src::AbstractArray)
-  @assert axes(dest) == axes(src)
-  sparse_map!(identity, dest, src)
-  return dest
-end
-
-function sparse_copyto!(dest::AbstractArray, src::AbstractArray)
-  sparse_map!(identity, dest, src)
-  return dest
-end
-
-function sparse_permutedims!(dest::AbstractArray, src::AbstractArray, perm)
-  sparse_copyto!(dest, PermutedDimsArray(src, perm))
-  return dest
-end
diff --git a/src/sparsearrayinterface/densearray.jl b/src/sparsearrayinterface/densearray.jl
deleted file mode 100644
index 3b98c18..0000000
--- a/src/sparsearrayinterface/densearray.jl
+++ /dev/null
@@ -1,12 +0,0 @@
-# Generic functionality for converting to a
-# dense array, trying to preserve information
-# about the array (such as which device it is on).
-# TODO: Maybe call `densecopy`?
-# TODO: Make sure this actually preserves the device,
-# maybe use `TypeParameterAccessors.unwrap_array_type`.
-function densearray(a::AbstractArray)
-  # TODO: `set_ndims(unwrap_array_type(a), ndims(a))(a)`
-  # Maybe define `densetype(a) = set_ndims(unwrap_array_type(a), ndims(a))`.
-  # Or could use `unspecify_parameters(unwrap_array_type(a))(a)`.
-  return Array(a)
-end
diff --git a/src/sparsearrayinterface/indexing.jl b/src/sparsearrayinterface/indexing.jl
deleted file mode 100644
index 5f8c1ca..0000000
--- a/src/sparsearrayinterface/indexing.jl
+++ /dev/null
@@ -1,216 +0,0 @@
-using ArrayLayouts: ArrayLayouts
-
-# An index into the storage of the sparse array.
-struct StorageIndex{I}
-  i::I
-end
-index(i::StorageIndex) = i.i
-
-# Indicate if the index into the sparse array is
-# stored or not.
-abstract type MaybeStoredIndex{I} end
-
-# An index into a stored value of the sparse array.
-# Stores both the index into the outer array
-# as well as into the underlying storage.
-struct StoredIndex{Iouter,Istorage} <: MaybeStoredIndex{Iouter}
-  iouter::Iouter
-  istorage::StorageIndex{Istorage}
-end
-index(i::StoredIndex) = i.iouter
-StorageIndex(i::StoredIndex) = i.istorage
-
-stored_length(a::AbstractArray) = length(sparse_storage(a))
-
-struct NotStoredIndex{Iouter} <: MaybeStoredIndex{Iouter}
-  iouter::Iouter
-end
-index(i::NotStoredIndex) = i.iouter
-
-function MaybeStoredIndex(a::AbstractArray, I)
-  return MaybeStoredIndex(I, index_to_storage_index(a, I))
-end
-MaybeStoredIndex(I, I_storage) = StoredIndex(I, StorageIndex(I_storage))
-MaybeStoredIndex(I, I_storage::Nothing) = NotStoredIndex(I)
-
-# Convert the index into an index into the storage.
-# Return `NotStoredIndex(I)` if it isn't in the storage.
-storage_index(a::AbstractArray, I...) = MaybeStoredIndex(a, I...)
-
-function storage_indices(a::AbstractArray)
-  return eachindex(sparse_storage(a))
-end
-
-# Derived
-function index_to_storage_index(a::AbstractArray{<:Any,N}, I::Vararg{Int,N}) where {N}
-  return index_to_storage_index(a, CartesianIndex(I))
-end
-
-function sparse_getindex(a::AbstractArray, I::NotStoredIndex)
-  return getindex_notstored(a, index(I))
-end
-
-function sparse_getindex(a::AbstractArray, I::StoredIndex)
-  return sparse_getindex(a, StorageIndex(I))
-end
-
-function sparse_getindex(a::AbstractArray, I::StorageIndex)
-  return sparse_storage(a)[index(I)]
-end
-
-function sparse_getindex(a::AbstractArray{<:Any,N}, I::Vararg{Int,N}) where {N}
-  return sparse_getindex(a, CartesianIndex(I))
-end
-
-function sparse_getindex(a::AbstractArray{<:Any,N}, I::CartesianIndex{N}) where {N}
-  return _sparse_getindex(a, I)
-end
-
-# Ambiguity with linear indexing
-function sparse_getindex(a::AbstractArray{<:Any,1}, I::CartesianIndex{1})
-  return _sparse_getindex(a, I)
-end
-
-# Implementation of element access
-function _sparse_getindex(a::AbstractArray{<:Any,N}, I::CartesianIndex{N}) where {N}
-  @boundscheck checkbounds(a, I)
-  return sparse_getindex(a, storage_index(a, I))
-end
-
-# Handle trailing indices or linear indexing
-function sparse_getindex(a::AbstractArray, I::Vararg{Int})
-  return sparse_getindex(a, CartesianIndex(I))
-end
-
-# Fix ambiguity error.
-function sparse_getindex(a::AbstractArray{<:Any,0})
-  return sparse_getindex(a, CartesianIndex())
-end
-
-# Linear indexing
-function sparse_getindex(a::AbstractArray, I::CartesianIndex{1})
-  return sparse_getindex(a, CartesianIndices(a)[I])
-end
-
-# Handle trailing indices
-function sparse_getindex(a::AbstractArray, I::CartesianIndex)
-  t = Tuple(I)
-  length(t) < ndims(a) && error("Not enough indices passed")
-  I′ = ntuple(i -> t[i], ndims(a))
-  @assert all(i -> isone(I[i]), (ndims(a) + 1):length(I))
-  return _sparse_getindex(a, CartesianIndex(I′))
-end
-
-# Slicing
-function sparse_getindex(a::AbstractArray, I::AbstractVector...)
-  return copy(@view a[I...])
-end
-
-function ArrayLayouts.sub_materialize(::SparseLayout, a::AbstractArray, axes)
-  a_dest = similar(a, axes)
-  a_dest .= a
-  return a_dest
-end
-
-# Update a nonzero value
-function sparse_setindex!(a::AbstractArray, value, I::StorageIndex)
-  sparse_storage(a)[index(I)] = value
-  return a
-end
-
-# Implementation of element access
-function _sparse_setindex!(a::AbstractArray{<:Any,N}, value, I::CartesianIndex{N}) where {N}
-  @boundscheck checkbounds(a, I)
-  sparse_setindex!(a, value, storage_index(a, I))
-  return a
-end
-
-# Ambiguity with linear indexing
-function sparse_setindex!(a::AbstractArray{<:Any,1}, value, I::CartesianIndex{1})
-  _sparse_setindex!(a, value, I)
-  return a
-end
-
-# Handle trailing indices or linear indexing
-function sparse_setindex!(a::AbstractArray, value, I::Vararg{Int})
-  sparse_setindex!(a, value, CartesianIndex(I))
-  return a
-end
-
-# Fix ambiguity error
-function sparse_setindex!(a::AbstractArray, value)
-  sparse_setindex!(a, value, CartesianIndex())
-  return a
-end
-
-# Linear indexing
-function sparse_setindex!(a::AbstractArray, value, I::CartesianIndex{1})
-  sparse_setindex!(a, value, CartesianIndices(a)[I])
-  return a
-end
-
-# Slicing
-# TODO: Make this handle more general slicing operations,
-# base it off of `ArrayLayouts.sub_materialize`.
-function sparse_setindex!(a::AbstractArray, value, I::AbstractUnitRange...)
-  inds = CartesianIndices(I)
-  for i in stored_indices(value)
-    if i in CartesianIndices(inds)
-      a[inds[i]] = value[i]
-    end
-  end
-  return a
-end
-
-# Handle trailing indices
-function sparse_setindex!(a::AbstractArray, value, I::CartesianIndex)
-  t = Tuple(I)
-  length(t) < ndims(a) && error("Not enough indices passed")
-  I′ = ntuple(i -> t[i], ndims(a))
-  @assert all(i -> isone(I[i]), (ndims(a) + 1):length(I))
-  return _sparse_setindex!(a, value, CartesianIndex(I′))
-end
-
-function sparse_setindex!(a::AbstractArray, value, I::StoredIndex)
-  sparse_setindex!(a, value, StorageIndex(I))
-  return a
-end
-
-function sparse_setindex!(a::AbstractArray, value, I::NotStoredIndex)
-  setindex_notstored!(a, value, index(I))
-  return a
-end
-
-# isassigned
-function sparse_isassigned(a::AbstractArray{<:Any,N}, I::CartesianIndex{N}) where {N}
-  return sparse_isassigned(a, Tuple(I)...)
-end
-function sparse_isassigned(a::AbstractArray, I::Integer...)
-  # Check trailing dimensions are one. This is needed in generic
-  # AbstractArray show when `a isa AbstractVector`.
-  all(d -> isone(I[d]), (ndims(a) + 1):length(I)) || return false
-  return all(dim -> I[dim] ∈ axes(a, dim), 1:ndims(a))
-end
-
-# A set of indices into the storage of the sparse array.
-struct StorageIndices{I}
-  i::I
-end
-indices(i::StorageIndices) = i.i
-
-function sparse_getindex(a::AbstractArray, I::StorageIndices{Colon})
-  return sparse_storage(a)
-end
-
-function sparse_getindex(a::AbstractArray, I::StorageIndices)
-  return error("Not implemented")
-end
-
-function sparse_setindex!(a::AbstractArray, value, I::StorageIndices{Colon})
-  sparse_storage(a) .= value
-  return a
-end
-
-function sparse_setindex!(a::AbstractArray, value, I::StorageIndices)
-  return error("Not implemented")
-end
diff --git a/src/sparsearrayinterface/interface.jl b/src/sparsearrayinterface/interface.jl
deleted file mode 100644
index ed5f7cb..0000000
--- a/src/sparsearrayinterface/interface.jl
+++ /dev/null
@@ -1,18 +0,0 @@
-# Also look into:
-# https://juliaarrays.github.io/ArrayInterface.jl/stable/sparsearrays/
-
-# Minimal sparse array interface.
-# Data structure of the stored (generally nonzero) values.
-# By default assume it is dense, so all values are stored.
-sparse_storage(a::AbstractArray) = a
-
-# Minimal sparse array interface.
-# Map an index into the stored data to a CartesianIndex of the
-# outer array.
-storage_index_to_index(a::AbstractArray, I) = I
-
-# Minimal interface
-# Map a `CartesianIndex` to an index/key into the nonzero data structure
-# returned by `storage`.
-# Return `nothing` if the index corresponds to a structural zero (unstored) value.
-index_to_storage_index(a::AbstractArray{<:Any,N}, I::CartesianIndex{N}) where {N} = I
diff --git a/src/sparsearrayinterface/interface_optional.jl b/src/sparsearrayinterface/interface_optional.jl
deleted file mode 100644
index 710cd25..0000000
--- a/src/sparsearrayinterface/interface_optional.jl
+++ /dev/null
@@ -1,44 +0,0 @@
-# Optional interface.
-
-# Function for computing unstored zero values.
-getindex_zero_function(::AbstractArray) = Zero()
-
-# Change the function for computing unstored values
-set_getindex_zero_function(a::AbstractArray, f) = error("Not implemented")
-
-function getindex_notstored(a::AbstractArray, I)
-  return getindex_zero_function(a)(a, I)
-end
-
-# Optional interface.
-# Insert a new value into a location
-# where there is not a stored value.
-# Some types (like `Diagonal`) may not support this.
-function setindex_notstored!(a::AbstractArray, value, I)
-  iszero(value) && return a
-  return throw(ArgumentError("Can't set nonzero values of $(typeof(a))."))
-end
-
-# Optional interface.
-# Iterates over the indices of `a` where there are stored values.
-# Can overload for faster iteration when there is more structure,
-# for example with DiagonalArrays.
-function stored_indices(a::AbstractArray)
-  return Iterators.map(Inz -> storage_index_to_index(a, Inz), storage_indices(a))
-end
-
-# Empty the sparse storage if possible.
-# Array types should overload `Base.dataids` to opt-in
-# to aliasing detection with `Base.mightalias`
-# to avoid emptying an input array in the case of `sparse_map!`.
-# `dropall!` is used to zero out the output array.
-# See also `Base.unalias` and `Base.unaliascopy`.
-# Interface is inspired by Base `SparseArrays.droptol!`
-# and `SparseArrays.dropzeros!`, and is like
-# `dropall!(a) = SparseArrays.droptol!(a, Inf)`.
-dropall!(a::AbstractArray) = a
-
-# Overload
-function sparse_similar(a::AbstractArray, elt::Type, dims::Tuple{Vararg{Int}})
-  return similar(a, elt, dims)
-end
diff --git a/src/sparsearrayinterface/map.jl b/src/sparsearrayinterface/map.jl
deleted file mode 100644
index 0f9d9aa..0000000
--- a/src/sparsearrayinterface/map.jl
+++ /dev/null
@@ -1,120 +0,0 @@
-using Base.Broadcast: BroadcastStyle, combine_styles
-using Compat: allequal
-using LinearAlgebra: LinearAlgebra
-
-# Represents a value that isn't stored
-# Used to hijack dispatch
-struct NotStoredValue{Value}
-  value::Value
-end
-value(v::NotStoredValue) = v.value
-stored_length(::NotStoredValue) = false
-Base.:*(x::Number, y::NotStoredValue) = false
-Base.:*(x::NotStoredValue, y::Number) = false
-Base.:/(x::NotStoredValue, y::Number) = false
-Base.:+(::NotStoredValue, ::NotStoredValue...) = false
-Base.:-(::NotStoredValue, ::NotStoredValue...) = false
-Base.:+(x::Number, ::NotStoredValue...) = x
-Base.iszero(::NotStoredValue) = true
-Base.isreal(::NotStoredValue) = true
-Base.conj(x::NotStoredValue) = conj(value(x))
-Base.iterate(x::NotStoredValue) = (x, nothing)
-Base.mapreduce(f, op, x::NotStoredValue) = f(x)
-Base.zero(x::NotStoredValue) = zero(value(x))
-LinearAlgebra.norm(x::NotStoredValue, p::Real=2) = zero(value(x))
-
-notstored_index(a::AbstractArray) = NotStoredIndex(first(eachindex(a)))
-
-# Get some not-stored value
-function get_notstored(a::AbstractArray)
-  return sparse_getindex(a, notstored_index(a))
-end
-
-function apply_notstored(f, as::Vararg{AbstractArray})
-  return apply(f, NotStoredValue.(get_notstored.(as))...)
-end
-
-function apply(f, xs::Vararg{NotStoredValue})
-  return f(xs...)
-  #return try
-  #  return f(xs...)
-  #catch
-  #  f(value(x))
-  #end
-end
-
-# Test if the function preserves zero values and therefore
-# preserves the sparsity structure.
-function preserves_zero(f, as...)
-  # return iszero(f(map(get_notstored, as)...))
-  return iszero(apply_notstored(f, as...))
-end
-
-# Map a subset of indices
-function sparse_map_indices!(f, a_dest::AbstractArray, indices, as::AbstractArray...)
-  for I in indices
-    a_dest[I] = f(map(a -> a[I], as)...)
-  end
-  return a_dest
-end
-
-# Overload for custom `stored_indices` types.
-function promote_indices(I1, I2)
-  return union(I1, I2)
-end
-
-function promote_indices(I1, I2, Is...)
-  return promote_indices(promote_indices(I1, I2), Is...)
-end
-
-# Base case
-promote_indices(I) = I
-
-function promote_stored_indices(as::AbstractArray...)
-  return promote_indices(stored_indices.(as)...)
-end
-
-function sparse_map_stored!(f, a_dest::AbstractArray, as::AbstractArray...)
-  # Need to zero out the destination.
-  sparse_zero!(a_dest)
-  Is = promote_stored_indices(as...)
-  sparse_map_indices!(f, a_dest, Is, as...)
-  return a_dest
-end
-
-# Handle nonzero case, fill all values.
-function sparse_map_all!(f, a_dest::AbstractArray, as::AbstractArray...)
-  Is = eachindex(a_dest)
-  sparse_map_indices!(f, a_dest, Is, as...)
-  return a_dest
-end
-
-function sparse_map!(f, a_dest::AbstractArray, as::AbstractArray...)
-  return sparse_map!(combine_styles(as...), f, a_dest, as...)
-end
-
-function sparse_map!(::BroadcastStyle, f, a_dest::AbstractArray, as::AbstractArray...)
-  @assert allequal(axes.((a_dest, as...)))
-  if preserves_zero(f, as...)
-    # Remove aliases to avoid overwriting inputs.
-    as = map(a -> Base.unalias(a_dest, a), as)
-    sparse_map_stored!(f, a_dest, as...)
-  else
-    sparse_map_all!(f, a_dest, as...)
-  end
-  return a_dest
-end
-
-# `f::typeof(norm)`, `op::typeof(max)` used by `norm`.
-function reduce_init(f, op, a)
-  return f(zero(eltype(a)))
-end
-
-# TODO: Generalize to multiple arguements.
-# TODO: Define `sparse_mapreducedim!`.
-function sparse_mapreduce(f, op, a::AbstractArray; init=reduce_init(f, op, a), kwargs...)
-  output = mapreduce(f, op, sparse_storage(a); init, kwargs...)
-  f_notstored = apply_notstored(f, a)
-  @assert isequal(op(output, eltype(output)(f_notstored)), output)
-  return output
-end
diff --git a/src/sparsearrayinterface/vectorinterface.jl b/src/sparsearrayinterface/vectorinterface.jl
deleted file mode 100644
index 787ec83..0000000
--- a/src/sparsearrayinterface/vectorinterface.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-using VectorInterface: zerovector!
-
-##################################################
-# TODO: Move to `DictionariesVectorInterfaceExt`.
-using VectorInterface: VectorInterface, zerovector!, zerovector!!
-using Dictionaries: AbstractDictionary
-
-function VectorInterface.zerovector!(x::AbstractDictionary{<:Number})
-  return fill!(x, zero(scalartype(x)))
-end
-function VectorInterface.zerovector!(x::AbstractDictionary)
-  T = eltype(x)
-  for I in eachindex(x)
-    if isbitstype(T) || isassigned(x, I)
-      x[I] = zerovector!!(x[I])
-    else
-      x[I] = zero(eltype(x))
-    end
-  end
-  return x
-end
-##################################################
-
-function sparse_zerovector!(a::AbstractArray)
-  dropall!(a)
-  zerovector!(sparse_storage(a))
-  return a
-end
diff --git a/src/sparsearrayinterface/wrappers.jl b/src/sparsearrayinterface/wrappers.jl
deleted file mode 100644
index 27c493f..0000000
--- a/src/sparsearrayinterface/wrappers.jl
+++ /dev/null
@@ -1,51 +0,0 @@
-using NestedPermutedDimsArrays: NestedPermutedDimsArray
-
-## PermutedDimsArray
-
-const AnyPermutedDimsArray{T,N,perm,iperm,P} = Union{
-  PermutedDimsArray{T,N,perm,iperm,P},NestedPermutedDimsArray{T,N,perm,iperm,P}
-}
-
-# TODO: Use `TypeParameterAccessors`.
-perm(::AnyPermutedDimsArray{<:Any,<:Any,Perm}) where {Perm} = Perm
-iperm(::AnyPermutedDimsArray{<:Any,<:Any,<:Any,IPerm}) where {IPerm} = IPerm
-
-# TODO: Use `Base.PermutedDimsArrays.genperm` or
-# https://github.com/jipolanco/StaticPermutations.jl?
-genperm(v, perm) = map(j -> v[j], perm)
-genperm(v::CartesianIndex, perm) = CartesianIndex(map(j -> Tuple(v)[j], perm))
-
-function storage_index_to_index(a::AnyPermutedDimsArray, I)
-  return genperm(storage_index_to_index(parent(a), I), perm(a))
-end
-
-function index_to_storage_index(
-  a::AnyPermutedDimsArray{<:Any,N}, I::CartesianIndex{N}
-) where {N}
-  return index_to_storage_index(parent(a), genperm(I, perm(a)))
-end
-
-# TODO: Add `getindex_zero_function` definition?
-
-## SubArray
-
-function map_index(
-  indices::Tuple{Vararg{Any,N}}, cartesian_index::CartesianIndex{N}
-) where {N}
-  index = Tuple(cartesian_index)
-  new_index = ntuple(length(indices)) do i
-    findfirst(==(index[i]), indices[i])
-  end
-  any(isnothing, new_index) && return nothing
-  return CartesianIndex(new_index)
-end
-
-function storage_index_to_index(a::SubArray, I)
-  return storage_index_to_index(parent(a), I)
-end
-
-function index_to_storage_index(a::SubArray{<:Any,N}, I::CartesianIndex{N}) where {N}
-  return index_to_storage_index(parent(a), I)
-end
-
-getindex_zero_function(a::SubArray) = getindex_zero_function(parent(a))
diff --git a/src/sparsearrayinterface/zero.jl b/src/sparsearrayinterface/zero.jl
deleted file mode 100644
index 72899b4..0000000
--- a/src/sparsearrayinterface/zero.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-# Represents a zero value and an index
-# TODO: Rename `GetIndexZero`?
-struct Zero end
-(f::Zero)(a::AbstractArray, I) = f(eltype(a), I)
-(::Zero)(type::Type, I) = zero(type)
diff --git a/src/wrappers.jl b/src/wrappers.jl
new file mode 100644
index 0000000..a7123a3
--- /dev/null
+++ b/src/wrappers.jl
@@ -0,0 +1,137 @@
+parentvalue_to_value(a::AbstractArray, value) = value
+value_to_parentvalue(a::AbstractArray, value) = value
+eachstoredparentindex(a::AbstractArray) = eachstoredindex(parent(a))
+storedparentvalues(a::AbstractArray) = storedvalues(parent(a))
+parentindex_to_index(a::AbstractArray, I::CartesianIndex) = error()
+function parentindex_to_index(a::AbstractArray, I::Int...)
+  return Tuple(parentindex_to_index(a, CartesianIndex(I)))
+end
+index_to_parentindex(a::AbstractArray, I::CartesianIndex) = error()
+function index_to_parentindex(a::AbstractArray, I::Int...)
+  return Tuple(index_to_parentindex(a, CartesianIndex(I)))
+end
+
+function cartesianindex_reverse(I::CartesianIndex)
+  return CartesianIndex(reverse(Tuple(I)))
+end
+tuple_oneto(n) = ntuple(identity, n)
+
+# TODO: Use `Base.PermutedDimsArrays.genperm` or
+# https://github.com/jipolanco/StaticPermutations.jl?
+genperm(v, perm) = map(j -> v[j], perm)
+
+using LinearAlgebra: Adjoint
+function parentindex_to_index(a::Adjoint, I::CartesianIndex)
+  return cartesianindex_reverse(I)
+end
+function index_to_parentindex(a::Adjoint, I::CartesianIndex)
+  return cartesianindex_reverse(I)
+end
+function parentvalue_to_value(a::Adjoint, value)
+  return adjoint(value)
+end
+function value_to_parentvalue(a::Adjoint, value)
+  return adjoint(value)
+end
+
+perm(::PermutedDimsArray{<:Any,<:Any,p}) where {p} = p
+iperm(::PermutedDimsArray{<:Any,<:Any,<:Any,ip}) where {ip} = ip
+function index_to_parentindex(a::PermutedDimsArray, I::CartesianIndex)
+  return CartesianIndex(genperm(I, iperm(a)))
+end
+function parentindex_to_index(a::PermutedDimsArray, I::CartesianIndex)
+  return CartesianIndex(genperm(I, perm(a)))
+end
+
+using Base: ReshapedArray
+function parentindex_to_index(a::ReshapedArray, I::CartesianIndex)
+  return CartesianIndices(size(a))[LinearIndices(parent(a))[I]]
+end
+function index_to_parentindex(a::ReshapedArray, I::CartesianIndex)
+  return CartesianIndices(parent(a))[LinearIndices(size(a))[I]]
+end
+
+function eachstoredparentindex(a::SubArray)
+  return filter(eachstoredindex(parent(a))) do I
+    return all(d -> I[d] ∈ parentindices(a)[d], 1:ndims(parent(a)))
+  end
+end
+function index_to_parentindex(a::SubArray, I::CartesianIndex)
+  return CartesianIndex(Base.reindex(parentindices(a), Tuple(I)))
+end
+function parentindex_to_index(a::SubArray, I::CartesianIndex)
+  nonscalardims = filter(tuple_oneto(ndims(parent(a)))) do d
+    return !(parentindices(a)[d] isa Real)
+  end
+  return CartesianIndex(
+    map(nonscalardims) do d
+      return findfirst(==(I[d]), parentindices(a)[d])
+    end,
+  )
+end
+## TODO: Use this and something similar for `Dictionary` to make a faster
+## implementation of `storedvalues(::SubArray)`.
+## function valuesview(d::Dict, keys)
+##   return @view d.vals[[Base.ht_keyindex(d, key) for key in keys]]
+## end
+function storedparentvalues(a::SubArray)
+  # We use `StoredValues` rather than `@view`/`SubArray` so that
+  # it gets interpreted as a dense array.
+  return StoredValues(parent(a), collect(eachstoredparentindex(a)))
+end
+
+using LinearAlgebra: Transpose
+function parentindex_to_index(a::Transpose, I::CartesianIndex)
+  return cartesianindex_reverse(I)
+end
+function index_to_parentindex(a::Transpose, I::CartesianIndex)
+  return cartesianindex_reverse(I)
+end
+function parentvalue_to_value(a::Transpose, value)
+  return transpose(value)
+end
+function value_to_parentvalue(a::Transpose, value)
+  return transpose(value)
+end
+
+# TODO: Turn these into `AbstractWrappedSparseArrayInterface` functions?
+for type in (:Adjoint, :PermutedDimsArray, :ReshapedArray, :SubArray, :Transpose)
+  @eval begin
+    @interface ::AbstractSparseArrayInterface storedvalues(a::$type) = storedparentvalues(a)
+    @interface ::AbstractSparseArrayInterface function isstored(a::$type, I::Int...)
+      return isstored(parent(a), index_to_parentindex(a, I...)...)
+    end
+    @interface ::AbstractSparseArrayInterface function eachstoredindex(a::$type)
+      # TODO: Make lazy with `Iterators.map`.
+      return map(collect(eachstoredparentindex(a))) do I
+        return parentindex_to_index(a, I)
+      end
+    end
+    @interface ::AbstractSparseArrayInterface function getstoredindex(a::$type, I::Int...)
+      return parentvalue_to_value(
+        a, getstoredindex(parent(a), index_to_parentindex(a, I...)...)
+      )
+    end
+    @interface ::AbstractSparseArrayInterface function getunstoredindex(a::$type, I::Int...)
+      return parentvalue_to_value(
+        a, getunstoredindex(parent(a), index_to_parentindex(a, I...)...)
+      )
+    end
+    @interface ::AbstractSparseArrayInterface function setstoredindex!(
+      a::$type, value, I::Int...
+    )
+      setstoredindex!(
+        parent(a), value_to_parentvalue(a, value), index_to_parentindex(a, I...)...
+      )
+      return a
+    end
+    @interface ::AbstractSparseArrayInterface function setunstoredindex!(
+      a::$type, value, I::Int...
+    )
+      setunstoredindex!(
+        parent(a), value_to_parentvalue(a, value), index_to_parentindex(a, I...)...
+      )
+      return a
+    end
+  end
+end
diff --git a/test/Project.toml b/test/Project.toml
index 74d11d0..9a4e6ac 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,12 +1,11 @@
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+Derive = "a07dfc7f-7d04-4eb5-84cc-a97f051f655a"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NestedPermutedDimsArrays = "2c2a8ec4-3cfc-4276-aa3e-1307b4294e58"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/SparseArraysBaseTestUtils/AbstractSparseArrays.jl b/test/SparseArraysBaseTestUtils/AbstractSparseArrays.jl
deleted file mode 100644
index ccf8e32..0000000
--- a/test/SparseArraysBaseTestUtils/AbstractSparseArrays.jl
+++ /dev/null
@@ -1,74 +0,0 @@
-module AbstractSparseArrays
-using ArrayLayouts: ArrayLayouts, MatMulMatAdd, MemoryLayout, MulAdd
-using SparseArraysBase: SparseArraysBase, AbstractSparseArray, Zero
-
-struct SparseArray{T,N,Zero} <: AbstractSparseArray{T,N}
-  data::Vector{T}
-  dims::Tuple{Vararg{Int,N}}
-  index_to_dataindex::Dict{CartesianIndex{N},Int}
-  dataindex_to_index::Vector{CartesianIndex{N}}
-  zero::Zero
-end
-function SparseArray{T,N}(dims::Tuple{Vararg{Int,N}}; zero=Zero()) where {T,N}
-  return SparseArray{T,N,typeof(zero)}(
-    T[], dims, Dict{CartesianIndex{N},Int}(), Vector{CartesianIndex{N}}(), zero
-  )
-end
-function SparseArray{T,N}(dims::Vararg{Int,N}; kwargs...) where {T,N}
-  return SparseArray{T,N}(dims; kwargs...)
-end
-function SparseArray{T}(dims::Tuple{Vararg{Int}}; kwargs...) where {T}
-  return SparseArray{T,length(dims)}(dims; kwargs...)
-end
-function SparseArray{T}(::UndefInitializer, dims::Tuple{Vararg{Int}}; kwargs...) where {T}
-  return SparseArray{T}(dims; kwargs...)
-end
-SparseArray{T}(dims::Vararg{Int}; kwargs...) where {T} = SparseArray{T}(dims; kwargs...)
-
-# ArrayLayouts interface
-struct SparseLayout <: MemoryLayout end
-ArrayLayouts.MemoryLayout(::Type{<:SparseArray}) = SparseLayout()
-function Base.similar(::MulAdd{<:SparseLayout,<:SparseLayout}, elt::Type, axes)
-  return similar(SparseArray{elt}, axes)
-end
-function ArrayLayouts.materialize!(
-  m::MatMulMatAdd{<:SparseLayout,<:SparseLayout,<:SparseLayout}
-)
-  α, a1, a2, β, a_dest = m.α, m.A, m.B, m.β, m.C
-  SparseArraysBase.sparse_mul!(a_dest, a1, a2, α, β)
-  return a_dest
-end
-
-# AbstractArray interface
-Base.size(a::SparseArray) = a.dims
-function Base.similar(a::SparseArray, elt::Type, dims::Tuple{Vararg{Int}})
-  return SparseArray{elt}(dims)
-end
-
-# Minimal interface
-SparseArraysBase.getindex_zero_function(a::SparseArray) = a.zero
-SparseArraysBase.sparse_storage(a::SparseArray) = a.data
-function SparseArraysBase.index_to_storage_index(
-  a::SparseArray{<:Any,N}, I::CartesianIndex{N}
-) where {N}
-  return get(a.index_to_dataindex, I, nothing)
-end
-SparseArraysBase.storage_index_to_index(a::SparseArray, I) = a.dataindex_to_index[I]
-function SparseArraysBase.setindex_notstored!(
-  a::SparseArray{<:Any,N}, value, I::CartesianIndex{N}
-) where {N}
-  push!(a.data, value)
-  push!(a.dataindex_to_index, I)
-  a.index_to_dataindex[I] = length(a.data)
-  return a
-end
-
-# Empty the storage, helps with efficiency in `map!` to drop
-# zeros.
-function SparseArraysBase.dropall!(a::SparseArray)
-  empty!(a.data)
-  empty!(a.index_to_dataindex)
-  empty!(a.dataindex_to_index)
-  return a
-end
-end
diff --git a/test/SparseArraysBaseTestUtils/DiagonalArrays.jl b/test/SparseArraysBaseTestUtils/DiagonalArrays.jl
deleted file mode 100644
index 96b062c..0000000
--- a/test/SparseArraysBaseTestUtils/DiagonalArrays.jl
+++ /dev/null
@@ -1,95 +0,0 @@
-module DiagonalArrays
-using SparseArraysBase: SparseArraysBase
-
-struct DiagonalArray{T,N} <: AbstractArray{T,N}
-  data::Vector{T}
-  dims::Tuple{Vararg{Int,N}}
-end
-function DiagonalArray{T,N}(::UndefInitializer, dims::Tuple{Vararg{Int,N}}) where {T,N}
-  return DiagonalArray{T,N}(Vector{T}(undef, minimum(dims)), dims)
-end
-function DiagonalArray{T,N}(::UndefInitializer, dims::Vararg{Int,N}) where {T,N}
-  return DiagonalArray{T,N}(undef, dims)
-end
-function DiagonalArray{T}(::UndefInitializer, dims::Tuple{Vararg{Int}}) where {T}
-  return DiagonalArray{T,length(dims)}(undef, dims)
-end
-function DiagonalArray{T}(::UndefInitializer, dims::Vararg{Int}) where {T}
-  return DiagonalArray{T}(undef, dims)
-end
-
-# AbstractArray interface
-Base.size(a::DiagonalArray) = a.dims
-function Base.getindex(a::DiagonalArray, I...)
-  return SparseArraysBase.sparse_getindex(a, I...)
-end
-function Base.setindex!(a::DiagonalArray, I...)
-  return SparseArraysBase.sparse_setindex!(a, I...)
-end
-function Base.similar(a::DiagonalArray, elt::Type, dims::Tuple{Vararg{Int}})
-  return DiagonalArray{elt}(undef, dims)
-end
-
-# Minimal interface
-SparseArraysBase.sparse_storage(a::DiagonalArray) = a.data
-function SparseArraysBase.index_to_storage_index(
-  a::DiagonalArray{<:Any,N}, I::CartesianIndex{N}
-) where {N}
-  !allequal(Tuple(I)) && return nothing
-  return first(Tuple(I))
-end
-function SparseArraysBase.storage_index_to_index(a::DiagonalArray, I)
-  return CartesianIndex(ntuple(Returns(I), ndims(a)))
-end
-function SparseArraysBase.sparse_similar(
-  a::DiagonalArray, elt::Type, dims::Tuple{Vararg{Int}}
-)
-  return Array{elt}(undef, dims)
-end
-function SparseArraysBase.sparse_similar(a::DiagonalArray, elt::Type, dims::Tuple{Int})
-  return similar(a, elt, dims)
-end
-
-# Broadcasting
-function Broadcast.BroadcastStyle(arraytype::Type{<:DiagonalArray})
-  return SparseArraysBase.SparseArrayStyle{ndims(arraytype)}()
-end
-
-# Base
-function Base.iszero(a::DiagonalArray)
-  return SparseArraysBase.sparse_iszero(a)
-end
-function Base.isreal(a::DiagonalArray)
-  return SparseArraysBase.sparse_isreal(a)
-end
-function Base.zero(a::DiagonalArray)
-  return SparseArraysBase.sparse_zero(a)
-end
-function Base.one(a::DiagonalArray)
-  return SparseArraysBase.sparse_one(a)
-end
-function Base.:(==)(a1::DiagonalArray, a2::DiagonalArray)
-  return SparseArraysBase.sparse_isequal(a1, a2)
-end
-function Base.reshape(a::DiagonalArray, dims::Tuple{Vararg{Int}})
-  return SparseArraysBase.sparse_reshape(a, dims)
-end
-
-# Map
-function Base.map!(f, dest::AbstractArray, src::DiagonalArray)
-  SparseArraysBase.sparse_map!(f, dest, src)
-  return dest
-end
-function Base.copy!(dest::AbstractArray, src::DiagonalArray)
-  SparseArraysBase.sparse_copy!(dest, src)
-  return dest
-end
-function Base.copyto!(dest::AbstractArray, src::DiagonalArray)
-  SparseArraysBase.sparse_copyto!(dest, src)
-  return dest
-end
-function Base.permutedims!(dest::AbstractArray, src::DiagonalArray, perm)
-  SparseArraysBase.sparse_permutedims!(dest, src, perm)
-  return dest
-end
-end
diff --git a/test/SparseArraysBaseTestUtils/SparseArrays.jl b/test/SparseArraysBaseTestUtils/SparseArrays.jl
deleted file mode 100644
index 827ce01..0000000
--- a/test/SparseArraysBaseTestUtils/SparseArrays.jl
+++ /dev/null
@@ -1,166 +0,0 @@
-module SparseArrays
-using LinearAlgebra: LinearAlgebra
-using SparseArraysBase: SparseArraysBase, Zero
-
-struct SparseArray{T,N,Zero} <: AbstractArray{T,N}
-  data::Vector{T}
-  dims::Tuple{Vararg{Int,N}}
-  index_to_dataindex::Dict{CartesianIndex{N},Int}
-  dataindex_to_index::Vector{CartesianIndex{N}}
-  zero::Zero
-end
-function SparseArray{T,N}(dims::Tuple{Vararg{Int,N}}; zero=Zero()) where {T,N}
-  return SparseArray{T,N,typeof(zero)}(
-    T[], dims, Dict{CartesianIndex{N},Int}(), Vector{CartesianIndex{N}}(), zero
-  )
-end
-function SparseArray{T,N}(dims::Vararg{Int,N}; kwargs...) where {T,N}
-  return SparseArray{T,N}(dims; kwargs...)
-end
-function SparseArray{T}(dims::Tuple{Vararg{Int}}; kwargs...) where {T}
-  return SparseArray{T,length(dims)}(dims; kwargs...)
-end
-function SparseArray{T}(::UndefInitializer, dims::Tuple{Vararg{Int}}; kwargs...) where {T}
-  return SparseArray{T}(dims; kwargs...)
-end
-SparseArray{T}(dims::Vararg{Int}; kwargs...) where {T} = SparseArray{T}(dims; kwargs...)
-
-# LinearAlgebra interface
-function LinearAlgebra.mul!(
-  a_dest::AbstractMatrix,
-  a1::SparseArray{<:Any,2},
-  a2::SparseArray{<:Any,2},
-  α::Number,
-  β::Number,
-)
-  SparseArraysBase.sparse_mul!(a_dest, a1, a2, α, β)
-  return a_dest
-end
-
-function LinearAlgebra.dot(a1::SparseArray, a2::SparseArray)
-  return SparseArraysBase.sparse_dot(a1, a2)
-end
-
-# AbstractArray interface
-Base.size(a::SparseArray) = a.dims
-function Base.similar(a::SparseArray, elt::Type, dims::Tuple{Vararg{Int}})
-  return SparseArray{elt}(dims)
-end
-
-function Base.getindex(a::SparseArray, I...)
-  return SparseArraysBase.sparse_getindex(a, I...)
-end
-function Base.setindex!(a::SparseArray, value, I...)
-  return SparseArraysBase.sparse_setindex!(a, value, I...)
-end
-function Base.fill!(a::SparseArray, value)
-  return SparseArraysBase.sparse_fill!(a, value)
-end
-
-# Minimal interface
-SparseArraysBase.getindex_zero_function(a::SparseArray) = a.zero
-SparseArraysBase.sparse_storage(a::SparseArray) = a.data
-function SparseArraysBase.index_to_storage_index(
-  a::SparseArray{<:Any,N}, I::CartesianIndex{N}
-) where {N}
-  return get(a.index_to_dataindex, I, nothing)
-end
-SparseArraysBase.storage_index_to_index(a::SparseArray, I) = a.dataindex_to_index[I]
-function SparseArraysBase.setindex_notstored!(
-  a::SparseArray{<:Any,N}, value, I::CartesianIndex{N}
-) where {N}
-  push!(a.data, value)
-  push!(a.dataindex_to_index, I)
-  a.index_to_dataindex[I] = length(a.data)
-  return a
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-using SparseArraysBase: perm, stored_indices
-function SparseArraysBase.stored_indices(
-  a::PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:SparseArray}
-)
-  return Iterators.map(
-    I -> CartesianIndex(map(i -> I[i], perm(a))), stored_indices(parent(a))
-  )
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-using SparseArraysBase: sparse_storage
-function SparseArraysBase.sparse_storage(
-  a::PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:SparseArray}
-)
-  return sparse_storage(parent(a))
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-using NestedPermutedDimsArrays: NestedPermutedDimsArray
-function SparseArraysBase.stored_indices(
-  a::NestedPermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:SparseArray}
-)
-  return Iterators.map(
-    I -> CartesianIndex(map(i -> I[i], perm(a))), stored_indices(parent(a))
-  )
-end
-
-# TODO: Make this into a generic definition of all `AbstractArray`?
-using NestedPermutedDimsArrays: NestedPermutedDimsArray
-using SparseArraysBase: sparse_storage
-function SparseArraysBase.sparse_storage(
-  a::NestedPermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:SparseArray}
-)
-  return sparse_storage(parent(a))
-end
-
-# Empty the storage, helps with efficiency in `map!` to drop
-# zeros.
-function SparseArraysBase.dropall!(a::SparseArray)
-  empty!(a.data)
-  empty!(a.index_to_dataindex)
-  empty!(a.dataindex_to_index)
-  return a
-end
-
-# Broadcasting
-function Broadcast.BroadcastStyle(arraytype::Type{<:SparseArray})
-  return SparseArraysBase.SparseArrayStyle{ndims(arraytype)}()
-end
-
-# Map
-function Base.map!(f, dest::AbstractArray, src::SparseArray)
-  SparseArraysBase.sparse_map!(f, dest, src)
-  return dest
-end
-function Base.copy!(dest::AbstractArray, src::SparseArray)
-  SparseArraysBase.sparse_copy!(dest, src)
-  return dest
-end
-function Base.copyto!(dest::AbstractArray, src::SparseArray)
-  SparseArraysBase.sparse_copyto!(dest, src)
-  return dest
-end
-function Base.permutedims!(dest::AbstractArray, src::SparseArray, perm)
-  SparseArraysBase.sparse_permutedims!(dest, src, perm)
-  return dest
-end
-
-# Base
-function Base.:(==)(a1::SparseArray, a2::SparseArray)
-  return SparseArraysBase.sparse_isequal(a1, a2)
-end
-function Base.reshape(a::SparseArray, dims::Tuple{Vararg{Int}})
-  return SparseArraysBase.sparse_reshape(a, dims)
-end
-function Base.iszero(a::SparseArray)
-  return SparseArraysBase.sparse_iszero(a)
-end
-function Base.isreal(a::SparseArray)
-  return SparseArraysBase.sparse_isreal(a)
-end
-function Base.zero(a::SparseArray)
-  return SparseArraysBase.sparse_zero(a)
-end
-function Base.one(a::SparseArray)
-  return SparseArraysBase.sparse_one(a)
-end
-end
diff --git a/test/SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl b/test/SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl
deleted file mode 100644
index 9067f9d..0000000
--- a/test/SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-module SparseArraysBaseTestUtils
-include("AbstractSparseArrays.jl")
-include("DiagonalArrays.jl")
-include("SparseArrays.jl")
-end
diff --git a/test/basics/test_basics.jl b/test/basics/test_basics.jl
new file mode 100644
index 0000000..d968acd
--- /dev/null
+++ b/test/basics/test_basics.jl
@@ -0,0 +1,62 @@
+using Adapt: adapt
+using JLArrays: JLArray, @allowscalar
+using SparseArraysBase:
+  SparseArraysBase,
+  eachstoredindex,
+  getstoredindex,
+  getunstoredindex,
+  isstored,
+  setstoredindex!,
+  setunstoredindex!,
+  storedlength,
+  storedpairs,
+  storedvalues
+using Test: @test, @testset
+
+elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
+arrayts = (Array, JLArray)
+@testset "SparseArraysBase (arraytype=$arrayt, eltype=$elt)" for arrayt in arrayts,
+  elt in elts
+
+  dev(x) = adapt(arrayt, x)
+
+  n = 2
+  a = dev(randn(elt, n, n))
+  @test storedlength(a) == length(a)
+  for indexstyle in (IndexLinear(), IndexCartesian())
+    for I in eachindex(indexstyle, a)
+      @test isstored(a, I)
+    end
+  end
+  @test eachstoredindex(a) == eachindex(a)
+  # TODO: We should be specializing these for dense/strided arrays,
+  # probably we can have a trait for that. It could be based
+  # on the `ArrayLayouts.MemoryLayout`.
+  @allowscalar @test storedvalues(a) == a
+  @allowscalar @test storedpairs(a) == collect(pairs(vec(a)))
+  @allowscalar for I in eachindex(a)
+    @test getstoredindex(a, I) == a[I]
+    @test iszero(getunstoredindex(a, I))
+  end
+  @allowscalar for I in eachindex(IndexCartesian(), a)
+    @test getstoredindex(a, I) == a[I]
+    @test iszero(getunstoredindex(a, I))
+  end
+
+  a = dev(randn(elt, n, n))
+  for I in ((1, 2), (CartesianIndex(1, 2),))
+    b = copy(a)
+    value = randn(elt)
+    @allowscalar setstoredindex!(b, value, I...)
+    @allowscalar b[I...] == value
+  end
+
+  # TODO: Should `setunstoredindex!` error by default
+  # if the value at that index is already stored?
+  a = dev(randn(elt, n, n))
+  for I in ((1, 2), (CartesianIndex(1, 2),))
+    b = copy(a)
+    value = randn(elt)
+    @test_throws ErrorException setunstoredindex!(b, value, I...)
+  end
+end
diff --git a/test/basics/test_sparsearraydok.jl b/test/basics/test_sparsearraydok.jl
new file mode 100644
index 0000000..f2302eb
--- /dev/null
+++ b/test/basics/test_sparsearraydok.jl
@@ -0,0 +1,171 @@
+using Adapt: adapt
+using ArrayLayouts: zero!
+using JLArrays: JLArray, @allowscalar
+using SparseArraysBase:
+  SparseArraysBase,
+  SparseArrayDOK,
+  eachstoredindex,
+  getstoredindex,
+  getunstoredindex,
+  isstored,
+  setstoredindex!,
+  setunstoredindex!,
+  storedlength,
+  storedpairs,
+  storedvalues
+using Test: @test, @testset
+
+elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
+# arrayts = (Array, JLArray)
+arrayts = (Array,)
+@testset "SparseArrayDOK (arraytype=$arrayt, eltype=$elt)" for arrayt in arrayts,
+  elt in elts
+
+  dev(x) = adapt(arrayt, x)
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  @test a isa SparseArrayDOK{elt,2}
+  @test size(a) == (2, 2)
+  @test a[1, 1] == 0
+  @test a[1, 1, 1] == 0
+  @test a[1, 2] == 12
+  @test a[1, 2, 1] == 12
+  @test storedlength(a) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  for b in (similar(a, Float32, (3, 3)), similar(a, Float32, Base.OneTo.((3, 3))))
+    @test b isa SparseArrayDOK{Float32,2}
+    @test b == zeros(Float32, 3, 3)
+    @test size(b) == (3, 3)
+  end
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = similar(a)
+  bc = Broadcast.Broadcasted(x -> 2x, (a,))
+  copyto!(b, bc)
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == [0 24; 0 0]
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(3, 3, 3)
+  a[1, 2, 3] = 123
+  b = permutedims(a, (2, 3, 1))
+  @test b isa SparseArrayDOK{elt,3}
+  @test b[2, 3, 1] == 123
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = copy(a')
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == [0 0; 12 0]
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = map(x -> 2x, a)
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == [0 24; 0 0]
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = a * a'
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == [144 0; 0 0]
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = a .+ 2 .* a'
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == [0 12; 24 0]
+  @test storedlength(b) == 2
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = a[1:2, 2]
+  @test b isa SparseArrayDOK{elt,1}
+  @test b == [12, 0]
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  @test iszero(a)
+  a[2, 1] = 21
+  a[1, 2] = 12
+  @test !iszero(a)
+  @test isreal(a)
+  @test sum(a) == 33
+  @test mapreduce(x -> 2x, +, a) == 66
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = similar(a)
+  copyto!(b, a)
+  @test b isa SparseArrayDOK{elt,2}
+  @test b == a
+  @test b[1, 2] == 12
+  @test storedlength(b) == 1
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a .= 2
+  @test storedlength(a) == length(a)
+  for I in eachindex(a)
+    @test a[I] == 2
+  end
+
+  a = SparseArrayDOK{elt}(2, 2)
+  fill!(a, 2)
+  @test storedlength(a) == length(a)
+  for I in eachindex(a)
+    @test a[I] == 2
+  end
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  zero!(a)
+  @test iszero(a)
+  @test iszero(storedlength(a))
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = zero(a)
+  @test b isa SparseArrayDOK{elt,2}
+  @test iszero(b)
+  @test iszero(storedlength(b))
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = SparseArrayDOK{elt}(4, 4)
+  b[2:3, 2:3] .= a
+  @test isone(storedlength(b))
+  @test b[2, 3] == 12
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = SparseArrayDOK{elt}(4, 4)
+  b[2:3, 2:3] = a
+  @test isone(storedlength(b))
+  @test b[2, 3] == 12
+
+  a = SparseArrayDOK{elt}(2, 2)
+  a[1, 2] = 12
+  b = SparseArrayDOK{elt}(4, 4)
+  c = @view b[2:3, 2:3]
+  c .= a
+  @test isone(storedlength(b))
+  @test b[2, 3] == 12
+
+  a1 = SparseArrayDOK{elt}(2, 2)
+  a1[1, 2] = 12
+  a2 = SparseArrayDOK{elt}(2, 2)
+  a2[2, 1] = 21
+  b = cat(a1, a2; dims=(1, 2))
+  @test b isa SparseArrayDOK{elt,2}
+  @test storedlength(b) == 2
+  @test b[1, 2] == 12
+  @test b[4, 3] == 21
+end
diff --git a/test/test_abstractsparsearray.jl b/test/test_abstractsparsearray.jl
deleted file mode 100644
index 7cc7c2f..0000000
--- a/test/test_abstractsparsearray.jl
+++ /dev/null
@@ -1,439 +0,0 @@
-@eval module $(gensym())
-using LinearAlgebra: dot, mul!, norm
-using SparseArraysBase: SparseArraysBase
-using NestedPermutedDimsArrays: NestedPermutedDimsArray
-include("SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl")
-using .SparseArraysBaseTestUtils.AbstractSparseArrays: AbstractSparseArrays
-using .SparseArraysBaseTestUtils.SparseArrays: SparseArrays
-using Test: @test, @testset
-@testset "AbstractSparseArray (arraytype=$SparseArray, eltype=$elt)" for SparseArray in (
-    AbstractSparseArrays.SparseArray, SparseArrays.SparseArray
-  ),
-  elt in (Float32, ComplexF32, Float64, ComplexF64)
-
-  a = SparseArray{elt}(2, 3)
-  @test size(a) == (2, 3)
-  @test axes(a) == (1:2, 1:3)
-  @test SparseArraysBase.sparse_storage(a) == elt[]
-  @test iszero(SparseArraysBase.stored_length(a))
-  @test collect(SparseArraysBase.stored_indices(a)) == CartesianIndex{2}[]
-  @test iszero(a)
-  @test iszero(norm(a))
-  for I in eachindex(a)
-    @test iszero(a)
-  end
-  for I in CartesianIndices(a)
-    @test isassigned(a, Tuple(I)...)
-    @test isassigned(a, I)
-  end
-  @test !isassigned(a, 0, 1)
-  @test !isassigned(a, CartesianIndex(0, 1))
-  @test !isassigned(a, 1, 4)
-  @test !isassigned(a, CartesianIndex(1, 4))
-
-  a = SparseArray{elt}(2, 3)
-  fill!(a, 0)
-  @test size(a) == (2, 3)
-  @test iszero(a)
-  @test iszero(SparseArraysBase.stored_length(a))
-
-  a_dense = SparseArraysBase.densearray(a)
-  @test a_dense == a
-  @test a_dense isa Array{elt,ndims(a)}
-
-  a = SparseArray{elt}(2, 3)
-  fill!(a, 2)
-  @test size(a) == (2, 3)
-  @test !iszero(a)
-  @test SparseArraysBase.stored_length(a) == length(a)
-  for I in eachindex(a)
-    @test a[I] == 2
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  @test a[1, 2] == 12
-  @test a[3] == 12 # linear indexing
-  @test size(a) == (2, 3)
-  @test axes(a) == (1:2, 1:3)
-  @test a[SparseArraysBase.StorageIndex(1)] == 12
-  @test SparseArraysBase.sparse_storage(a) == elt[12]
-  @test isone(SparseArraysBase.stored_length(a))
-  @test collect(SparseArraysBase.stored_indices(a)) == [CartesianIndex(1, 2)]
-  @test !iszero(a)
-  @test !iszero(norm(a))
-  for I in eachindex(a)
-    if I == CartesianIndex(1, 2)
-      @test a[I] == 12
-    else
-      @test iszero(a[I])
-    end
-  end
-  for I in CartesianIndices(a)
-    @test isassigned(a, Tuple(I)...)
-    @test isassigned(a, I)
-  end
-  @test !isassigned(a, 0, 1)
-  @test !isassigned(a, CartesianIndex(0, 1))
-  @test !isassigned(a, 1, 4)
-  @test !isassigned(a, CartesianIndex(1, 4))
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  a = map(x -> 2x, a)
-  for I in eachindex(a)
-    if I == CartesianIndex(1, 2)
-      @test a[I] == 2 * 12
-    else
-      @test iszero(a[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 2, 2)
-  a[1, 2, 2] = 122
-  a_r = reshape(a, 2, 4)
-  @test a_r[1, 4] == a[1, 2, 2] == 122
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  a = zero(a)
-  @test size(a) == (2, 3)
-  @test iszero(SparseArraysBase.stored_length(a))
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = SparseArray{elt}(2, 3)
-  b[2, 1] = 21
-  @test a == a
-  @test a ≠ b
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  @test isreal(a)
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = randn(elt)
-  b = copy(a)
-  conj!(b)
-  for I in eachindex(a)
-    @test conj(a[I]) == b[I]
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = randn(elt)
-  b = conj(a)
-  for I in eachindex(a)
-    @test conj(a[I]) == b[I]
-  end
-
-  if !(elt <: Real)
-    a = SparseArray{elt}(2, 3)
-    a[1, 2] = 12 + 12im
-    @test !isreal(a)
-  end
-
-  a = SparseArray{elt}(2, 2)
-  a[1, 2] = 12
-  a = one(a)
-  @test size(a) == (2, 2)
-  @test isone(a[1, 1])
-  @test isone(a[2, 2])
-  @test iszero(a[1, 2])
-  @test iszero(a[2, 1])
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  a = zero(a)
-  @test size(a) == (2, 3)
-  @test iszero(SparseArraysBase.stored_length(a))
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  a = copy(a)
-  @test size(a) == (2, 3)
-  @test axes(a) == (1:2, 1:3)
-  @test SparseArraysBase.sparse_storage(a) == elt[12]
-  @test isone(SparseArraysBase.stored_length(a))
-  @test SparseArraysBase.storage_indices(a) == 1:1
-  @test collect(SparseArraysBase.stored_indices(a)) == [CartesianIndex(1, 2)]
-  @test !iszero(a)
-  @test !iszero(norm(a))
-  for I in eachindex(a)
-    if I == CartesianIndex(1, 2)
-      @test a[I] == 12
-    else
-      @test iszero(a[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  a = 2 * a
-  @test size(a) == (2, 3)
-  @test axes(a) == (1:2, 1:3)
-  @test SparseArraysBase.sparse_storage(a) == elt[24]
-  @test isone(SparseArraysBase.stored_length(a))
-  @test collect(SparseArraysBase.stored_indices(a)) == [CartesianIndex(1, 2)]
-  @test !iszero(a)
-  @test !iszero(norm(a))
-  for I in eachindex(a)
-    if I == CartesianIndex(1, 2)
-      @test a[I] == 24
-    else
-      @test iszero(a[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = SparseArray{elt}(2, 3)
-  b[2, 1] = 21
-  c = a + b
-  @test size(c) == (2, 3)
-  @test axes(c) == (1:2, 1:3)
-  @test SparseArraysBase.sparse_storage(c) == elt[12, 21]
-  @test SparseArraysBase.stored_length(c) == 2
-  @test collect(SparseArraysBase.stored_indices(c)) ==
-    [CartesianIndex(1, 2), CartesianIndex(2, 1)]
-  @test !iszero(c)
-  @test !iszero(norm(c))
-  for I in eachindex(c)
-    if I == CartesianIndex(1, 2)
-      @test c[I] == 12
-    elseif I == CartesianIndex(2, 1)
-      @test c[I] == 21
-    else
-      @test iszero(c[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = permutedims(a, (2, 1))
-  @test size(b) == (3, 2)
-  @test axes(b) == (1:3, 1:2)
-  @test SparseArraysBase.sparse_storage(b) == elt[12]
-  @test SparseArraysBase.stored_length(b) == 1
-  @test collect(SparseArraysBase.stored_indices(b)) == [CartesianIndex(2, 1)]
-  @test !iszero(b)
-  @test !iszero(norm(b))
-  for I in eachindex(b)
-    if I == CartesianIndex(2, 1)
-      @test b[I] == 12
-    else
-      @test iszero(b[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = PermutedDimsArray(a, (2, 1))
-  @test size(b) == (3, 2)
-  @test axes(b) == (1:3, 1:2)
-  @test SparseArraysBase.sparse_storage(b) == elt[12]
-  @test SparseArraysBase.stored_length(b) == 1
-  @test collect(SparseArraysBase.stored_indices(b)) == [CartesianIndex(2, 1)]
-  @test !iszero(b)
-  @test !iszero(norm(b))
-  for I in eachindex(b)
-    if I == CartesianIndex(2, 1)
-      @test b[I] == 12
-    else
-      @test iszero(b[I])
-    end
-  end
-
-  a = SparseArray{Matrix{elt}}(
-    2, 3; zero=(a, I) -> (z = similar(eltype(a), 2, 3); fill!(z, false); z)
-  )
-  a[1, 2] = randn(elt, 2, 3)
-  b = NestedPermutedDimsArray(a, (2, 1))
-  @test size(b) == (3, 2)
-  @test axes(b) == (1:3, 1:2)
-  @test SparseArraysBase.sparse_storage(b) == [a[1, 2]]
-  @test SparseArraysBase.stored_length(b) == 1
-  @test collect(SparseArraysBase.stored_indices(b)) == [CartesianIndex(2, 1)]
-  @test !iszero(b)
-  @test !iszero(norm(b))
-  for I in eachindex(b)
-    if I == CartesianIndex(2, 1)
-      @test b[I] == permutedims(a[1, 2], (2, 1))
-    else
-      @test iszero(b[I])
-    end
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = randn(elt, 2, 3)
-  b .= a
-  @test a == b
-  for I in eachindex(a)
-    @test a[I] == b[I]
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = randn(elt, 2, 3)
-  b .= 2 .* a
-  @test 2 * a == b
-  for I in eachindex(a)
-    @test 2 * a[I] == b[I]
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = randn(elt, 2, 3)
-  b .= 2 .+ a
-  @test 2 .+ a == b
-  for I in eachindex(a)
-    @test 2 + a[I] == b[I]
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = randn(elt, 2, 3)
-  map!(x -> 2x, b, a)
-  @test 2 * a == b
-  for I in eachindex(a)
-    @test 2 * a[I] == b[I]
-  end
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = zeros(elt, 2, 3)
-  b[2, 1] = 21
-  @test Array(a) == a
-  @test a + b == Array(a) + b
-  @test b + a == Array(a) + b
-  @test b .+ 2 .* a == 2 * Array(a) + b
-  @test a .+ 2 .* b == Array(a) + 2b
-  @test a + b isa Matrix{elt}
-  @test b + a isa Matrix{elt}
-  @test SparseArraysBase.stored_length(a + b) == length(a)
-
-  a = SparseArray{elt}(2, 3)
-  a[1, 2] = 12
-  b = zeros(elt, 2, 3)
-  b[2, 1] = 21
-  a′ = copy(a)
-  a′ .+= b
-  @test a′ == a + b
-  # TODO: Should this be:
-  # ```julia
-  # @test SparseArraysBase.stored_length(a′) == 2
-  # ```
-  # ? I.e. should it only store the nonzero values?
-  @test SparseArraysBase.stored_length(a′) == 6
-
-  # Matrix multiplication
-  a1 = SparseArray{elt}(2, 3)
-  a1[1, 2] = 12
-  a1[2, 1] = 21
-  a2 = SparseArray{elt}(3, 4)
-  a2[1, 1] = 11
-  a2[2, 2] = 22
-  a_dest = a1 * a2
-  @test Array(a_dest) ≈ Array(a1) * Array(a2)
-  @test a_dest isa SparseArray{elt}
-  @test SparseArraysBase.stored_length(a_dest) == 2
-
-  # Dot product
-  a1 = SparseArray{elt}(4)
-  a1[1] = randn()
-  a1[3] = randn()
-  a2 = SparseArray{elt}(4)
-  a2[2] = randn()
-  a2[3] = randn()
-  a_dest = a1' * a2
-  @test a_dest isa elt
-  @test a_dest ≈ Array(a1)' * Array(a2)
-  @test a_dest ≈ dot(a1, a2)
-
-  # In-place matrix multiplication
-  a1 = SparseArray{elt}(2, 3)
-  a1[1, 2] = 12
-  a1[2, 1] = 21
-  a2 = SparseArray{elt}(3, 4)
-  a2[1, 1] = 11
-  a2[2, 2] = 22
-  a_dest = SparseArray{elt}(2, 4)
-  mul!(a_dest, a1, a2)
-  @test Array(a_dest) ≈ Array(a1) * Array(a2)
-  @test a_dest isa SparseArray{elt}
-  @test SparseArraysBase.stored_length(a_dest) == 2
-
-  # In-place matrix multiplication
-  a1 = SparseArray{elt}(2, 3)
-  a1[1, 2] = 12
-  a1[2, 1] = 21
-  a2 = SparseArray{elt}(3, 4)
-  a2[1, 1] = 11
-  a2[2, 2] = 22
-  a_dest = SparseArray{elt}(2, 4)
-  a_dest[1, 2] = 12
-  a_dest[2, 1] = 21
-  α = elt(2)
-  β = elt(3)
-  a_dest′ = copy(a_dest)
-  mul!(a_dest, a1, a2, α, β)
-  @test Array(a_dest) ≈ Array(a1) * Array(a2) * α + Array(a_dest′) * β
-  @test a_dest isa SparseArray{elt}
-  @test SparseArraysBase.stored_length(a_dest) == 2
-
-  # cat
-  a1 = SparseArray{elt}(2, 3)
-  a1[1, 2] = 12
-  a1[2, 1] = 21
-  a2 = SparseArray{elt}(2, 3)
-  a2[1, 1] = 11
-  a2[2, 2] = 22
-
-  a_dest = cat(a1, a2; dims=1)
-  @test size(a_dest) == (4, 3)
-  @test SparseArraysBase.stored_length(a_dest) == 4
-  @test a_dest[1, 2] == a1[1, 2]
-  @test a_dest[2, 1] == a1[2, 1]
-  @test a_dest[3, 1] == a2[1, 1]
-  @test a_dest[4, 2] == a2[2, 2]
-
-  a_dest = cat(a1, a2; dims=2)
-  @test size(a_dest) == (2, 6)
-  @test SparseArraysBase.stored_length(a_dest) == 4
-  @test a_dest[1, 2] == a1[1, 2]
-  @test a_dest[2, 1] == a1[2, 1]
-  @test a_dest[1, 4] == a2[1, 1]
-  @test a_dest[2, 5] == a2[2, 2]
-
-  a_dest = cat(a1, a2; dims=(1, 2))
-  @test size(a_dest) == (4, 6)
-  @test SparseArraysBase.stored_length(a_dest) == 4
-  @test a_dest[1, 2] == a1[1, 2]
-  @test a_dest[2, 1] == a1[2, 1]
-  @test a_dest[3, 4] == a2[1, 1]
-  @test a_dest[4, 5] == a2[2, 2]
-
-  ## # Sparse matrix of matrix multiplication
-  ## TODO: Make this work, seems to require
-  ## a custom zero constructor.
-  ## a1 = SparseArray{Matrix{elt}}(2, 3)
-  ## a1[1, 1] = zeros(elt, (2, 3))
-  ## a1[1, 2] = randn(elt, (2, 3))
-  ## a1[2, 1] = randn(elt, (2, 3))
-  ## a1[2, 2] = zeros(elt, (2, 3))
-  ## a2 = SparseArray{Matrix{elt}}(3, 4)
-  ## a2[1, 1] = randn(elt, (3, 4))
-  ## a2[1, 2] = zeros(elt, (3, 4))
-  ## a2[2, 2] = randn(elt, (3, 4))
-  ## a2[2, 2] = zeros(elt, (3, 4))
-  ## a_dest = SparseArray{Matrix{elt}}(2, 4)
-  ## a_dest[1, 1] = zeros(elt, (3, 4))
-  ## a_dest[1, 2] = zeros(elt, (3, 4))
-  ## a_dest[2, 2] = zeros(elt, (3, 4))
-  ## a_dest[2, 2] = zeros(elt, (3, 4))
-  ## mul!(a_dest, a1, a2)
-  ## @test Array(a_dest) ≈ Array(a1) * Array(a2)
-  ## @test a_dest isa SparseArray{Matrix{elt}}
-  ## @test SparseArraysBase.stored_length(a_dest) == 2
-end
-end
diff --git a/test/test_array.jl b/test/test_array.jl
deleted file mode 100644
index 8306eef..0000000
--- a/test/test_array.jl
+++ /dev/null
@@ -1,13 +0,0 @@
-@eval module $(gensym())
-using SparseArraysBase: SparseArraysBase
-include("SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl")
-using Test: @test, @testset
-@testset "Array (eltype=$elt)" for elt in (Float32, ComplexF32, Float64, ComplexF64)
-  a = randn(2, 3)
-  @test SparseArraysBase.sparse_storage(a) == a
-  @test SparseArraysBase.index_to_storage_index(a, CartesianIndex(1, 2)) ==
-    CartesianIndex(1, 2)
-  @test SparseArraysBase.storage_index_to_index(a, CartesianIndex(1, 2)) ==
-    CartesianIndex(1, 2)
-end
-end
diff --git a/test/test_diagonalarray.jl b/test/test_diagonalarray.jl
deleted file mode 100644
index 59bc470..0000000
--- a/test/test_diagonalarray.jl
+++ /dev/null
@@ -1,76 +0,0 @@
-@eval module $(gensym())
-using LinearAlgebra: norm
-using SparseArraysBase: SparseArraysBase
-include("SparseArraysBaseTestUtils/SparseArraysBaseTestUtils.jl")
-using .SparseArraysBaseTestUtils.DiagonalArrays: DiagonalArray
-using Test: @test, @testset, @test_throws
-@testset "DiagonalArray (eltype=$elt)" for elt in (Float32, ComplexF32, Float64, ComplexF64)
-  # TODO: Test `fill!`.
-
-  # Test
-  a = DiagonalArray{elt}(undef, 2, 3)
-  @test size(a) == (2, 3)
-  a[1, 1] = 11
-  a[2, 2] = 22
-  @test a[1, 1] == 11
-  @test a[2, 2] == 22
-  @test_throws ArgumentError a[1, 2] = 12
-  @test SparseArraysBase.storage_indices(a) == 1:2
-  @test collect(SparseArraysBase.stored_indices(a)) ==
-    [CartesianIndex(1, 1), CartesianIndex(2, 2)]
-  a[1, 2] = 0
-  @test a[1, 1] == 11
-  @test a[2, 2] == 22
-
-  a_dense = SparseArraysBase.densearray(a)
-  @test a_dense == a
-  @test a_dense isa Array{elt,ndims(a)}
-
-  b = similar(a)
-  @test b isa DiagonalArray
-  @test size(b) == (2, 3)
-
-  a = DiagonalArray(elt[1, 2, 3], (3, 3))
-  @test size(a) == (3, 3)
-  @test a[1, 1] == 1
-  @test a[2, 2] == 2
-  @test a[3, 3] == 3
-  @test a[SparseArraysBase.StorageIndex(1)] == 1
-  @test a[SparseArraysBase.StorageIndex(2)] == 2
-  @test a[SparseArraysBase.StorageIndex(3)] == 3
-  @test iszero(a[1, 2])
-
-  a = DiagonalArray(elt[1, 2, 3], (3, 3))
-  a = 2 * a
-  @test size(a) == (3, 3)
-  @test a[1, 1] == 2
-  @test a[2, 2] == 4
-  @test a[3, 3] == 6
-  @test iszero(a[1, 2])
-
-  a = DiagonalArray(elt[1, 2, 3], (3, 3))
-  a_r = reshape(a, 9)
-  @test a_r isa DiagonalArray{elt,1}
-  for I in LinearIndices(a)
-    @test a[I] == a_r[I]
-  end
-
-  # This needs `Base.reshape` with a custom destination
-  # calling `SparseArraysBase.sparse_reshape!`
-  # in order to specify an appropriate output
-  # type to work.
-  a = DiagonalArray(elt[1, 2], (2, 2, 2))
-  a_r = reshape(a, 2, 4)
-  @test a_r isa Matrix{elt}
-  for I in LinearIndices(a)
-    @test a[I] == a_r[I]
-  end
-
-  # Matrix multiplication!
-  a1 = DiagonalArray(elt[1, 2], (2, 2))
-  a2 = DiagonalArray(elt[2, 3], (2, 2))
-  a_dest = a1 * a2
-  @test Array(a_dest) ≈ Array(a1) * Array(a2)
-  @test a_dest isa DiagonalArray{elt}
-end
-end
diff --git a/test/test_sparsearraydok.jl b/test/test_sparsearraydok.jl
deleted file mode 100644
index eb43f28..0000000
--- a/test/test_sparsearraydok.jl
+++ /dev/null
@@ -1,138 +0,0 @@
-@eval module $(gensym())
-
-# TODO: Test:
-# zero (PermutedDimsArray)
-# Custom zero type
-# Slicing
-
-using Dictionaries: Dictionary
-using Test: @test, @testset, @test_broken
-using SparseArraysBase: SparseArraysBase, SparseArrayDOK, SparseMatrixDOK, @maybe_grow
-using SparseArraysBase: storage_indices, stored_length
-using SparseArrays: SparseMatrixCSC, nnz
-@testset "SparseArrayDOK (eltype=$elt)" for elt in
-                                            (Float32, ComplexF32, Float64, ComplexF64)
-  @testset "Basics" begin
-    a = SparseArrayDOK{elt}(3, 4)
-    @test a == SparseArrayDOK{elt}((3, 4))
-    @test a == SparseArrayDOK{elt}(undef, 3, 4)
-    @test a == SparseArrayDOK{elt}(undef, (3, 4))
-    @test iszero(a)
-    @test iszero(nnz(a))
-    @test stored_length(a) == nnz(a)
-    @test size(a) == (3, 4)
-    @test eltype(a) == elt
-    for I in eachindex(a)
-      @test iszero(a[I])
-      @test a[I] isa elt
-    end
-    @test isempty(storage_indices(a))
-
-    x12 = randn(elt)
-    x23 = randn(elt)
-    b = copy(a)
-    @test b isa SparseArrayDOK{elt}
-    @test iszero(b)
-    b[1, 2] = x12
-    b[2, 3] = x23
-    @test iszero(a)
-    @test !iszero(b)
-    @test b[1, 2] == x12
-    @test b[2, 3] == x23
-    @test iszero(stored_length(a))
-    @test stored_length(b) == 2
-  end
-  @testset "map/broadcast" begin
-    a = SparseArrayDOK{elt}(3, 4)
-    a[1, 1] = 11
-    a[3, 4] = 34
-    @test stored_length(a) == 2
-    b = 2 * a
-    @test stored_length(b) == 2
-    @test b[1, 1] == 2 * 11
-    @test b[3, 4] == 2 * 34
-  end
-  @testset "reshape" begin
-    a = SparseArrayDOK{elt}(2, 2, 2)
-    a[1, 2, 2] = 122
-    b = reshape(a, 2, 4)
-    @test b[1, 4] == 122
-  end
-  @testset "Matrix multiplication" begin
-    a1 = SparseArrayDOK{elt}(2, 3)
-    a1[1, 2] = 12
-    a1[2, 1] = 21
-    a2 = SparseArrayDOK{elt}(3, 4)
-    a2[1, 1] = 11
-    a2[2, 2] = 22
-    a2[3, 3] = 33
-    a_dest = a1 * a2
-    # TODO: Use `densearray` to make generic to GPU.
-    @test Array(a_dest) ≈ Array(a1) * Array(a2)
-    # TODO: Make this work with `ArrayLayouts`.
-    @test stored_length(a_dest) == 2
-    @test a_dest isa SparseMatrixDOK{elt}
-
-    a2 = randn(elt, (3, 4))
-    a_dest = a1 * a2
-    # TODO: Use `densearray` to make generic to GPU.
-    @test Array(a_dest) ≈ Array(a1) * Array(a2)
-    @test stored_length(a_dest) == 8
-    @test a_dest isa Matrix{elt}
-  end
-  @testset "SparseMatrixCSC" begin
-    a = SparseArrayDOK{elt}(2, 2)
-    a[1, 2] = 12
-    for (type, a′) in ((SparseMatrixCSC, a), (SparseArrayDOK, SparseMatrixCSC(a)))
-      b = type(a′)
-      @test b isa type{elt}
-      @test b[1, 2] == 12
-      @test isone(nnz(b))
-      for I in eachindex(b)
-        if I ≠ CartesianIndex(1, 2)
-          @test iszero(b[I])
-        end
-      end
-    end
-  end
-  @testset "Maybe Grow Feature" begin
-    a = SparseArrayDOK{elt,2}((0, 0))
-    SparseArraysBase.setindex_maybe_grow!(a, 230, 2, 3)
-    @test size(a) == (2, 3)
-    @test a[2, 3] == 230
-    # Test @maybe_grow macro
-    @maybe_grow a[5, 5] = 550
-    @test size(a) == (5, 5)
-    @test a[2, 3] == 230
-    @test a[5, 5] == 550
-    # Test that size remains same
-    # if we set at an index smaller than
-    # the maximum size:
-    @maybe_grow a[3, 4] = 340
-    @test size(a) == (5, 5)
-    @test a[2, 3] == 230
-    @test a[5, 5] == 550
-    @test a[3, 4] == 340
-    # Test vector case
-    v = SparseArrayDOK{elt,1}((0,))
-    @maybe_grow v[5] = 50
-    @test size(v) == (5,)
-    @test v[5] == 50
-    # Test setting from a variable (to test macro escaping)
-    i = 6
-    val = 60
-    @maybe_grow v[i] = val
-    @test v[i] == val
-    i, j = 1, 2
-    val = 120
-    @maybe_grow a[i, j] = val
-    @test a[i, j] == val
-  end
-  @testset "Test Lower Level Constructor" begin
-    d = Dictionary{CartesianIndex{2},elt}()
-    a = SparseArrayDOK(d, (2, 2), zero(elt))
-    a[1, 2] = 12.0
-    @test a[1, 2] == 12.0
-  end
-end
-end