diff --git a/Project.toml b/Project.toml
index cf3035bf..0a2cd552 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.2.28"
+version = "0.3.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/README.md b/README.md
index e2f415e4..72a1fb93 100644
--- a/README.md
+++ b/README.md
@@ -34,129 +34,17 @@ julia> Pkg.add("BlockSparseArrays")
 ## Examples
 
 ````julia
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
+using BlockArrays: Block
 using BlockSparseArrays: BlockSparseArray, blockstoredlength
-using Test: @test, @test_broken
-
-function main()
-  # Block dimensions
-  i1 = [2, 3]
-  i2 = [2, 3]
-
-  i_axes = (blockedrange(i1), blockedrange(i2))
-
-  function block_size(axes, block)
-    return length.(getindex.(axes, Block.(block.n)))
-  end
-
-  # Data
-  nz_blocks = Block.([(1, 1), (2, 2)])
-  nz_block_sizes = [block_size(i_axes, nz_block) for nz_block in nz_blocks]
-  nz_block_lengths = prod.(nz_block_sizes)
-
-  # Blocks with contiguous underlying data
-  d_data = BlockedVector(randn(sum(nz_block_lengths)), nz_block_lengths)
-  d_blocks = [
-    reshape(@view(d_data[Block(i)]), block_size(i_axes, nz_blocks[i])) for
-    i in 1:length(nz_blocks)
-  ]
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  # Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  # Access a block
-  @test b[Block(1, 1)] == d_blocks[1]
-
-  # Access a zero block, returns a zero matrix
-  @test b[Block(1, 2)] == zeros(2, 3)
-
-  # Set a zero block
-  a₁₂ = randn(2, 3)
-  b[Block(1, 2)] = a₁₂
-  @test b[Block(1, 2)] == a₁₂
-
-  # Matrix multiplication
-  @test b * b ≈ Array(b) * Array(b)
-
-  permuted_b = permutedims(b, (2, 1))
-  @test permuted_b isa BlockSparseArray
-  @test permuted_b == permutedims(Array(b), (2, 1))
-
-  @test b + b ≈ Array(b) + Array(b)
-  @test b + b isa BlockSparseArray
-  # TODO: Fix this, broken.
-  @test_broken blockstoredlength(b + b) == 2
-
-  scaled_b = 2b
-  @test scaled_b ≈ 2Array(b)
-  @test scaled_b isa BlockSparseArray
-
-  # TODO: Fix this, broken.
-  @test_broken reshape(b, ([4, 6, 6, 9],)) isa BlockSparseArray{<:Any,1}
-
-  return nothing
-end
-
-main()
-````
-
-# BlockSparseArrays.jl and BlockArrays.jl interface
-
-````julia
-using BlockArrays: BlockArrays, Block
-using BlockSparseArrays: BlockSparseArray
-
-i1 = [2, 3]
-i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
-B[Block(1, 1)] = randn(2, 2)
-B[Block(2, 2)] = randn(3, 3)
-
-# Minimal interface
-
-# Specifies the block structure
-@show collect.(BlockArrays.blockaxes(axes(B, 1)))
-
-# Index range of a block
-@show axes(B, 1)[Block(1)]
-
-# Last index of each block
-@show BlockArrays.blocklasts(axes(B, 1))
-
-# Find the block containing the index
-@show BlockArrays.findblock(axes(B, 1), 3)
-
-# Retrieve a block
-@show B[Block(1, 1)]
-@show BlockArrays.viewblock(B, Block(1, 1))
-
-# Check block bounds
-@show BlockArrays.blockcheckbounds(B, 2, 2)
-@show BlockArrays.blockcheckbounds(B, Block(2, 2))
-
-# Derived interface
-
-# Specifies the block structure
-@show collect(Iterators.product(BlockArrays.blockaxes(B)...))
-
-# Iterate over block views
-@show sum.(BlockArrays.eachblock(B))
-
-# Reshape into 1-d
-# TODO: Fix this, broken.
-# @show BlockArrays.blockvec(B)[Block(1)]
-
-# Array-of-array view
-@show BlockArrays.blocks(B)[1, 1] == B[Block(1, 1)]
-
-# Access an index within a block
-@show B[Block(1, 1)[1, 1]] == B[1, 1]
+using Test: @test
+
+a = BlockSparseArray{Float64}(undef, [2, 3], [2, 3])
+a[Block(1, 2)] = randn(2, 3)
+a[Block(2, 1)] = randn(3, 2)
+@test blockstoredlength(a) == 2
+b = a .+ 2 .* a'
+@test Array(b) ≈ Array(a) + 2 * Array(a')
+@test blockstoredlength(b) == 2
 ````
 
 ---
diff --git a/docs/Project.toml b/docs/Project.toml
index 6dca2c9c..1c7114e2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 Documenter = "1"
 Literate = "2"
diff --git a/examples/Project.toml b/examples/Project.toml
index 48948bf1..cf876f1c 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 Test = "1"
diff --git a/examples/README.jl b/examples/README.jl
index 2d00e5b0..ba0467b7 100644
--- a/examples/README.jl
+++ b/examples/README.jl
@@ -39,124 +39,14 @@ julia> Pkg.add("BlockSparseArrays")
 
 # ## Examples
 
-using BlockArrays: BlockArrays, BlockedVector, Block, blockedrange
+using BlockArrays: Block
 using BlockSparseArrays: BlockSparseArray, blockstoredlength
-using Test: @test, @test_broken
-
-function main()
-  ## Block dimensions
-  i1 = [2, 3]
-  i2 = [2, 3]
-
-  i_axes = (blockedrange(i1), blockedrange(i2))
-
-  function block_size(axes, block)
-    return length.(getindex.(axes, Block.(block.n)))
-  end
-
-  ## Data
-  nz_blocks = Block.([(1, 1), (2, 2)])
-  nz_block_sizes = [block_size(i_axes, nz_block) for nz_block in nz_blocks]
-  nz_block_lengths = prod.(nz_block_sizes)
-
-  ## Blocks with contiguous underlying data
-  d_data = BlockedVector(randn(sum(nz_block_lengths)), nz_block_lengths)
-  d_blocks = [
-    reshape(@view(d_data[Block(i)]), block_size(i_axes, nz_blocks[i])) for
-    i in 1:length(nz_blocks)
-  ]
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  ## Blocks with discontiguous underlying data
-  d_blocks = randn.(nz_block_sizes)
-  b = BlockSparseArray(nz_blocks, d_blocks, i_axes)
-
-  @test blockstoredlength(b) == 2
-
-  ## Access a block
-  @test b[Block(1, 1)] == d_blocks[1]
-
-  ## Access a zero block, returns a zero matrix
-  @test b[Block(1, 2)] == zeros(2, 3)
-
-  ## Set a zero block
-  a₁₂ = randn(2, 3)
-  b[Block(1, 2)] = a₁₂
-  @test b[Block(1, 2)] == a₁₂
-
-  ## Matrix multiplication
-  @test b * b ≈ Array(b) * Array(b)
-
-  permuted_b = permutedims(b, (2, 1))
-  @test permuted_b isa BlockSparseArray
-  @test permuted_b == permutedims(Array(b), (2, 1))
-
-  @test b + b ≈ Array(b) + Array(b)
-  @test b + b isa BlockSparseArray
-  ## TODO: Fix this, broken.
-  @test_broken blockstoredlength(b + b) == 2
-
-  scaled_b = 2b
-  @test scaled_b ≈ 2Array(b)
-  @test scaled_b isa BlockSparseArray
-
-  ## TODO: Fix this, broken.
-  @test_broken reshape(b, ([4, 6, 6, 9],)) isa BlockSparseArray{<:Any,1}
-
-  return nothing
-end
-
-main()
-
-# # BlockSparseArrays.jl and BlockArrays.jl interface
-
-using BlockArrays: BlockArrays, Block
-using BlockSparseArrays: BlockSparseArray
-
-i1 = [2, 3]
-i2 = [2, 3]
-B = BlockSparseArray{Float64}(i1, i2)
-B[Block(1, 1)] = randn(2, 2)
-B[Block(2, 2)] = randn(3, 3)
-
-## Minimal interface
-
-## Specifies the block structure
-@show collect.(BlockArrays.blockaxes(axes(B, 1)))
-
-## Index range of a block
-@show axes(B, 1)[Block(1)]
-
-## Last index of each block
-@show BlockArrays.blocklasts(axes(B, 1))
-
-## Find the block containing the index
-@show BlockArrays.findblock(axes(B, 1), 3)
-
-## Retrieve a block
-@show B[Block(1, 1)]
-@show BlockArrays.viewblock(B, Block(1, 1))
-
-## Check block bounds
-@show BlockArrays.blockcheckbounds(B, 2, 2)
-@show BlockArrays.blockcheckbounds(B, Block(2, 2))
-
-## Derived interface
-
-## Specifies the block structure
-@show collect(Iterators.product(BlockArrays.blockaxes(B)...))
-
-## Iterate over block views
-@show sum.(BlockArrays.eachblock(B))
-
-## Reshape into 1-d
-## TODO: Fix this, broken.
-## @show BlockArrays.blockvec(B)[Block(1)]
-
-## Array-of-array view
-@show BlockArrays.blocks(B)[1, 1] == B[Block(1, 1)]
-
-## Access an index within a block
-@show B[Block(1, 1)[1, 1]] == B[1, 1]
+using Test: @test
+
+a = BlockSparseArray{Float64}(undef, [2, 3], [2, 3])
+a[Block(1, 2)] = randn(2, 3)
+a[Block(2, 1)] = randn(3, 2)
+@test blockstoredlength(a) == 2
+b = a .+ 2 .* a'
+@test Array(b) ≈ Array(a) + 2 * Array(a')
+@test blockstoredlength(b) == 2
diff --git a/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl b/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl
index 376f53e8..49762edf 100644
--- a/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl
+++ b/ext/BlockSparseArraysGradedUnitRangesExt/BlockSparseArraysGradedUnitRangesExt.jl
@@ -17,7 +17,7 @@ function similar_blocksparse(
   return BlockSparseArray{
     elt,length(axes),similartype(unwrap_array_type(blocktype(a)), elt, axes)
   }(
-    axes
+    undef, axes
   )
 end
 
diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index d58a9b75..a350dfff 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -5,7 +5,8 @@ export BlockSparseArray,
   BlockSparseVector,
   blockstoredlength,
   eachblockstoredindex,
-  eachstoredblock
+  eachstoredblock,
+  sparsemortar
 
 # factorizations
 include("factorizations/svd.jl")
@@ -38,7 +39,6 @@ include("abstractblocksparsearray/cat.jl")
 include("abstractblocksparsearray/adapt.jl")
 
 # functions specifically for BlockSparseArray
-include("blocksparsearray/defaults.jl")
 include("blocksparsearray/blocksparsearray.jl")
 include("blocksparsearray/blockdiagonalarray.jl")
 
diff --git a/src/abstractblocksparsearray/arraylayouts.jl b/src/abstractblocksparsearray/arraylayouts.jl
index e02eb979..8b771208 100644
--- a/src/abstractblocksparsearray/arraylayouts.jl
+++ b/src/abstractblocksparsearray/arraylayouts.jl
@@ -38,7 +38,7 @@ function ArrayLayouts.sub_materialize(layout::BlockLayout{<:SparseLayout}, a, ax
   # TODO: Define `blocktype`/`blockstype` for `SubArray` wrapping `BlockSparseArray`.
   # TODO: Use `similar`?
   blocktype_a = blocktype(parent(a))
-  a_dest = BlockSparseArray{eltype(a),length(axes),blocktype_a}(axes)
+  a_dest = BlockSparseArray{eltype(a),length(axes),blocktype_a}(undef, axes)
   a_dest .= a
   return a_dest
 end
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 7d0f281b..05345ac3 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -1,7 +1,8 @@
-# type alias for block-diagonal
-using LinearAlgebra: Diagonal
+using BlockArrays: blockedrange
 using DiagonalArrays: DiagonalArrays, diagonal
+using LinearAlgebra: Diagonal
 
+# type alias for block-diagonal
 const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{
   T,A,Diagonal{A,V},Axes
 }
@@ -12,8 +13,8 @@ const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
 end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
-  return BlockSparseArray(
-    Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
+  return _BlockSparseArray(
+    Diagonal(blocks), blockedrange.((size.(blocks, 1), size.(blocks, 2)))
   )
 end
 
diff --git a/src/blocksparsearray/blocksparsearray.jl b/src/blocksparsearray/blocksparsearray.jl
index 44b8d774..9d7e570e 100644
--- a/src/blocksparsearray/blocksparsearray.jl
+++ b/src/blocksparsearray/blocksparsearray.jl
@@ -1,168 +1,225 @@
-using BlockArrays: BlockArrays, Block, BlockedUnitRange, blockedrange, blocklength
+using BlockArrays:
+  BlockArrays,
+  Block,
+  BlockedUnitRange,
+  UndefBlocksInitializer,
+  blockedrange,
+  blocklength,
+  undef_blocks
 using DerivableInterfaces: @interface
 using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
 
-# TODO: Delete this.
-## using BlockArrays: blocks
+"""
+    SparseArrayDOK{T}(undef_blocks, axes)
+    SparseArrayDOK{T,N}(undef_blocks, axes)
+
+Construct the block structure of an undefined BlockSparseArray that will have
+blocked axes `axes`.
+
+Note that `undef_blocks` is defined in
+[BlockArrays.jl](https://juliaarrays.github.io/BlockArrays.jl/stable/lib/public/#BlockArrays.undef_blocks)
+and should be imported from that package to use it as an input to this constructor.
+"""
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef, blocklength.(ax), GetUnstoredBlock(ax))
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer, ax::Vararg{AbstractUnitRange{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, blockedrange.(dims))
+end
+function SparseArraysBase.SparseArrayDOK{T,N}(
+  ::UndefBlocksInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, (dim1, dim_rest...))
+end
+
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, ax::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer, ax::Vararg{AbstractUnitRange{<:Integer},N}
+) where {T,N}
+  return SparseArrayDOK{T,N}(undef_blocks, ax)
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
+) where {T}
+  return SparseArrayDOK{T}(undef_blocks, blockedrange.(dims))
+end
+function SparseArraysBase.SparseArrayDOK{T}(
+  ::UndefBlocksInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
+) where {T}
+  return SparseArrayDOK{T}(undef_blocks, (dim1, dim_rest...))
+end
+
+function _BlockSparseArray end
 
 struct BlockSparseArray{
   T,
   N,
   A<:AbstractArray{T,N},
   Blocks<:AbstractArray{A,N},
-  Axes<:Tuple{Vararg{AbstractUnitRange,N}},
+  Axes<:Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
 } <: AbstractBlockSparseArray{T,N}
   blocks::Blocks
   axes::Axes
+  global @inline function _BlockSparseArray(
+    blocks::AbstractArray{<:AbstractArray{T,N},N},
+    axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
+  ) where {T,N}
+    Base.require_one_based_indexing(axes...)
+    Base.require_one_based_indexing(blocks)
+    return new{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
+  end
 end
 
 # TODO: Can this definition be shortened?
-const BlockSparseMatrix{T,A<:AbstractMatrix{T},Blocks<:AbstractMatrix{A},Axes<:Tuple{AbstractUnitRange,AbstractUnitRange}} = BlockSparseArray{
+const BlockSparseMatrix{T,A<:AbstractMatrix{T},Blocks<:AbstractMatrix{A},Axes<:Tuple{AbstractUnitRange{<:Integer},AbstractUnitRange{<:Integer}}} = BlockSparseArray{
   T,2,A,Blocks,Axes
 }
 
 # TODO: Can this definition be shortened?
-const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:Tuple{AbstractUnitRange}} = BlockSparseArray{
+const BlockSparseVector{T,A<:AbstractVector{T},Blocks<:AbstractVector{A},Axes<:Tuple{AbstractUnitRange{<:Integer}}} = BlockSparseArray{
   T,1,A,Blocks,Axes
 }
 
-function BlockSparseArray(
-  block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  blocks = default_blocks(block_data, axes)
-  return BlockSparseArray(blocks, axes)
-end
-
-function BlockSparseArray(
-  block_indices::Vector{<:Block{N}},
-  block_data::Vector{<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return BlockSparseArray(Dictionary(block_indices, block_data), axes)
-end
+"""
+    sparsemortar(blocks::AbstractArray{<:AbstractArray{T,N},N}, axes) -> ::BlockSparseArray{T,N}
 
-function BlockSparseArray{T,N,A,Blocks}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
-) where {T,N,A<:AbstractArray{T,N},Blocks<:AbstractArray{A,N}}
-  return BlockSparseArray{T,N,A,Blocks,typeof(axes)}(blocks, axes)
+Construct a block sparse array from a sparse array of arrays and specified blocked axes.
+The block sizes must be commensurate with the blocks of the axes.
+"""
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}},
+) where {T,N}
+  return _BlockSparseArray(blocks, axes)
 end
-
-function BlockSparseArray{T,N,A}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
-) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A,typeof(blocks)}(blocks, axes)
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  axes::Vararg{AbstractUnitRange{<:Integer},N},
+) where {T,N}
+  return sparsemortar(blocks, axes)
 end
-
-function BlockSparseArray{T,N}(
-  blocks::AbstractArray{<:AbstractArray{T,N},N}, axes::Tuple{Vararg{AbstractUnitRange,N}}
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N}
-  return BlockSparseArray{T,N,eltype(blocks),typeof(blocks),typeof(axes)}(blocks, axes)
+  return sparsemortar(blocks, blockedrange.(dims))
 end
-
-function BlockSparseArray{T,N}(
-  block_data::Dictionary{Block{N,Int},<:AbstractArray{T,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
+function sparsemortar(
+  blocks::AbstractArray{<:AbstractArray{T,N},N},
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N}
-  blocks = default_blocks(block_data, axes)
-  return BlockSparseArray{T,N}(blocks, axes)
+  return sparsemortar(blocks, (dim1, dim_rest...))
 end
 
+@doc """
+    BlockSparseArray{T}(undef, dims)
+    BlockSparseArray{T,N}(undef, dims)
+    BlockSparseArray{T,N,A}(undef, dims)
+
+Construct an uninitialized N-dimensional BlockSparseArray containing elements of type T. `dims` should be a list
+of block lengths in each dimension or a list of blocked ranges representing the axes.
+""" BlockSparseArray
+
 function BlockSparseArray{T,N,A}(
-  axes::Tuple{Vararg{AbstractUnitRange,N}}
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
 ) where {T,N,A<:AbstractArray{T,N}}
-  blocks = default_blocks(A, axes)
-  return BlockSparseArray{T,N,A}(blocks, axes)
+  return _BlockSparseArray(SparseArrayDOK{A}(undef_blocks, axes), axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  axes::Vararg{AbstractUnitRange,N}
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer},N}
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(axes)
+  return BlockSparseArray{T,N,A}(undef, axes)
 end
 
 function BlockSparseArray{T,N,A}(
-  dims::Tuple{Vararg{Vector{Int},N}}
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(blockedrange.(dims))
-end
-
-# Fix ambiguity error.
-function BlockSparseArray{T,0,A}(axes::Tuple{}) where {T,A<:AbstractArray{T,0}}
-  blocks = default_blocks(A, axes)
-  return BlockSparseArray{T,0,A}(blocks, axes)
+  return BlockSparseArray{T,N,A}(undef, blockedrange.(dims))
 end
 
 function BlockSparseArray{T,N,A}(
-  dims::Vararg{Vector{Int},N}
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
 ) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(dims)
-end
-
-function BlockSparseArray{T,N}(axes::Tuple{Vararg{AbstractUnitRange,N}}) where {T,N}
-  return BlockSparseArray{T,N,default_arraytype(T, axes)}(axes)
-end
-
-function BlockSparseArray{T,N}(axes::Vararg{AbstractUnitRange,N}) where {T,N}
-  return BlockSparseArray{T,N}(axes)
-end
-
-function BlockSparseArray{T,0}(axes::Tuple{}) where {T}
-  return BlockSparseArray{T,0,default_arraytype(T, axes)}(axes)
-end
-
-function BlockSparseArray{T,N}(dims::Tuple{Vararg{Vector{Int},N}}) where {T,N}
-  return BlockSparseArray{T,N}(blockedrange.(dims))
-end
-
-function BlockSparseArray{T,N}(dims::Vararg{Vector{Int},N}) where {T,N}
-  return BlockSparseArray{T,N}(dims)
+  return BlockSparseArray{T,N,A}(undef, (dim1, dim_rest...))
 end
 
-function BlockSparseArray{T}(dims::Tuple{Vararg{Vector{Int}}}) where {T}
-  return BlockSparseArray{T,length(dims)}(dims)
-end
-
-function BlockSparseArray{T}(axes::Tuple{Vararg{AbstractUnitRange}}) where {T}
-  return BlockSparseArray{T,length(axes)}(axes)
-end
-
-function BlockSparseArray{T}(axes::Tuple{}) where {T}
-  return BlockSparseArray{T,length(axes)}(axes)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
+) where {T,N}
+  return BlockSparseArray{T,N,Array{T,N}}(undef, axes)
 end
 
-function BlockSparseArray{T}(dims::Vararg{Vector{Int}}) where {T}
-  return BlockSparseArray{T}(dims)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer},N}
+) where {T,N}
+  return BlockSparseArray{T,N}(undef, axes)
 end
 
-function BlockSparseArray{T}(axes::Vararg{AbstractUnitRange}) where {T}
-  return BlockSparseArray{T}(axes)
+function BlockSparseArray{T,N}(
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
+) where {T,N}
+  return BlockSparseArray{T,N}(undef, blockedrange.(dims))
 end
 
-function BlockSparseArray{T}() where {T}
-  return BlockSparseArray{T}(())
+function BlockSparseArray{T,N}(
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
+) where {T,N}
+  return BlockSparseArray{T,N}(undef, (dim1, dim_rest...))
 end
 
-# undef
-function BlockSparseArray{T,N,A,Blocks}(
-  ::UndefInitializer, args...
-) where {T,N,A<:AbstractArray{T,N},Blocks<:AbstractArray{A,N}}
-  return BlockSparseArray{T,N,A,Blocks}(args...)
+function BlockSparseArray{T}(
+  ::UndefInitializer,
+  dims::Tuple{AbstractVector{<:Integer},Vararg{AbstractVector{<:Integer}}},
+) where {T}
+  return BlockSparseArray{T,length(dims)}(undef, dims)
 end
 
-function BlockSparseArray{T,N,A}(
-  ::UndefInitializer, args...
-) where {T,N,A<:AbstractArray{T,N}}
-  return BlockSparseArray{T,N,A}(args...)
+function BlockSparseArray{T}(
+  ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer}}}
+) where {T}
+  return BlockSparseArray{T,length(axes)}(undef, axes)
 end
 
-function BlockSparseArray{T,N}(::UndefInitializer, args...) where {T,N}
-  return BlockSparseArray{T,N}(args...)
+function BlockSparseArray{T}(
+  ::UndefInitializer,
+  dim1::AbstractVector{<:Integer},
+  dim_rest::AbstractVector{<:Integer}...,
+) where {T}
+  return BlockSparseArray{T}(undef, (dim1, dim_rest...))
 end
 
-function BlockSparseArray{T}(::UndefInitializer, args...) where {T}
-  return BlockSparseArray{T}(args...)
+function BlockSparseArray{T}(
+  ::UndefInitializer, axes::Vararg{AbstractUnitRange{<:Integer}}
+) where {T}
+  return BlockSparseArray{T}(undef, axes)
 end
 
 # Base `AbstractArray` interface
diff --git a/src/blocksparsearray/defaults.jl b/src/blocksparsearray/defaults.jl
deleted file mode 100644
index 788a9750..00000000
--- a/src/blocksparsearray/defaults.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-using BlockArrays: Block
-using Dictionaries: Dictionary
-using SparseArraysBase: SparseArrayDOK
-
-# Construct the sparse structure storing the blocks
-function default_blockdata(
-  block_data::Dictionary{<:CartesianIndex{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return error()
-end
-
-function default_blocks(
-  block_indices::Vector{<:Block{N}},
-  block_data::Vector{<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return default_blocks(Dictionary(block_indices, block_data), axes)
-end
-
-function default_blocks(
-  block_data::Dictionary{<:Block{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return default_blocks(blocks_to_cartesianindices(block_data), axes)
-end
-
-function default_arraytype(elt::Type, axes::Tuple{Vararg{AbstractUnitRange}})
-  return Array{elt,length(axes)}
-end
-
-function default_blocks(blocktype::Type, axes::Tuple{Vararg{AbstractUnitRange}})
-  block_data = Dictionary{Block{length(axes),Int},blocktype}()
-  return default_blocks(block_data, axes)
-end
-
-function default_blocks(
-  block_data::Dictionary{<:CartesianIndex{N},<:AbstractArray{<:Any,N}},
-  axes::Tuple{Vararg{AbstractUnitRange,N}},
-) where {N}
-  return SparseArrayDOK(block_data, blocklength.(axes), GetUnstoredBlock(axes))
-end
diff --git a/src/blocksparsearrayinterface/map.jl b/src/blocksparsearrayinterface/map.jl
index b97bff89..2df71338 100644
--- a/src/blocksparsearrayinterface/map.jl
+++ b/src/blocksparsearrayinterface/map.jl
@@ -94,23 +94,21 @@ function map_zero_dim! end
   return a_dest
 end
 
-# TODO: Decide what to do with these.
+# TODO: Do we need this function or can we just use `map`?
+# Probably it should be a special version of `map` where we
+# specify the function preserves zeros, i.e.
+# `map(f, a; preserves_zeros=true)` or `@preserves_zeros map(f, a)`.
 function map_stored_blocks(f, a::AbstractArray)
-  # TODO: Implement this as:
-  # ```julia
-  # mapped_blocks = SparseArraysInterface.map_stored(f, blocks(a))
-  # BlockSparseArray(mapped_blocks, axes(a))
-  # ```
-  # TODO: `block_stored_indices` should output `Indices` storing
-  # the stored Blocks, not a `Dictionary` from cartesian indices
-  # to Blocks.
-  bs = collect(eachblockstoredindex(a))
-  ds = map(b -> f(@view(a[b])), bs)
-  # We manually specify the block type using `Base.promote_op`
-  # since `a[b]` may not be inferrable. For example, if `blocktype(a)`
-  # is `Diagonal{Float64,Vector{Float64}}`, the non-stored blocks are `Matrix{Float64}`
-  # since they can't necessarily by `Diagonal` if there are rectangular blocks.
-  mapped_blocks = Dictionary{eltype(bs),eltype(ds)}(bs, ds)
-  # TODO: Use `similartype(typeof(a), eltype(eltype(mapped_blocks)))(...)`.
-  return BlockSparseArray(mapped_blocks, axes(a))
+  block_stored_indices = collect(eachblockstoredindex(a))
+  if isempty(block_stored_indices)
+    blocktype_a′ = Base.promote_op(f, blocktype(a))
+    return BlockSparseArray{eltype(blocktype_a′),ndims(a),blocktype_a′}(undef, axes(a))
+  end
+  stored_blocks = map(B -> f(@view!(a[B])), block_stored_indices)
+  blocktype_a′ = eltype(stored_blocks)
+  a′ = BlockSparseArray{eltype(blocktype_a′),ndims(a),blocktype_a′}(undef, axes(a))
+  for (B, b) in zip(block_stored_indices, stored_blocks)
+    a′[B] = b
+  end
+  return a′
 end
diff --git a/test/Project.toml b/test/Project.toml
index 4c93abdc..817fd103 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -26,7 +26,7 @@ Adapt = "4"
 Aqua = "0.8"
 ArrayLayouts = "1"
 BlockArrays = "1"
-BlockSparseArrays = "0.2"
+BlockSparseArrays = "0.3"
 DiagonalArrays = "0.2"
 GPUArraysCore = "0.2"
 GradedUnitRanges = "0.1"
diff --git a/test/test_basics.jl b/test/test_basics.jl
index da93d228..e39e2667 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -17,7 +17,8 @@ using BlockArrays:
   blocklengths,
   blocksize,
   blocksizes,
-  mortar
+  mortar,
+  undef_blocks
 using BlockSparseArrays:
   @view!,
   BlockSparseArray,
@@ -30,6 +31,7 @@ using BlockSparseArrays:
   eachstoredblock,
   blockstype,
   blocktype,
+  sparsemortar,
   view!
 using GPUArraysCore: @allowscalar
 using JLArrays: JLArray
@@ -48,20 +50,20 @@ arrayts = (Array, JLArray)
   dev(a) = adapt(arrayt, a)
   @testset "Broken" begin
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @test_broken a[:, 4]
 
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @test_broken a[:, [2, 4]]
     @test_broken a[[3, 5], [2, 4]]
 
     # TODO: Fix this and turn it into a proper test.
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a[Block(1, 1)] = dev(randn(elt, 2, 2))
     a[Block(2, 2)] = dev(randn(elt, 3, 3))
     @allowscalar @test a[2:4, 4] == Array(a)[2:4, 4]
@@ -85,10 +87,6 @@ arrayts = (Array, JLArray)
       ## BlockSparseMatrix{elt,Matrix{elt},SparseMatrixDOK{Matrix{elt}}}, # TODO
     )
       for args in (
-        bs,
-        (bs,),
-        blockedrange.(bs),
-        (blockedrange.(bs),),
         (undef, bs),
         (undef, bs...),
         (undef, blockedrange.(bs)),
@@ -117,10 +115,6 @@ arrayts = (Array, JLArray)
       ## BlockSparseVector{elt,Vector{elt},SparseVectorDOK{Vector{elt}}}, # TODO
     )
       for args in (
-        bs,
-        (bs,),
-        blockedrange.(bs),
-        (blockedrange.(bs),),
         (undef, bs),
         (undef, bs...),
         (undef, blockedrange.(bs)),
@@ -146,7 +140,7 @@ arrayts = (Array, JLArray)
     # TODO: This is difficult to determine just from type information.
     @test_broken blocktype(typeof(a)) <: SubArray{elt,2,arrayt{elt,2}}
 
-    a = BlockSparseMatrix{elt,arrayt{elt,2}}([1, 1], [1, 1])
+    a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [1, 1], [1, 1])
     @test (@constinferred blockstype(a)) <: SparseMatrixDOK{arrayt{elt,2}}
     @test (@constinferred blockstype(typeof(a))) <: SparseMatrixDOK{arrayt{elt,2}}
     @test (@constinferred blocktype(a)) <: arrayt{elt,2}
@@ -165,11 +159,29 @@ arrayts = (Array, JLArray)
     @test (@constinferred blocktype(a)) <: SubArray{elt,2,arrayt{elt,2}}
     # TODO: This is difficult to determine just from type information.
     @test_broken blocktype(typeof(a)) <: SubArray{elt,2,arrayt{elt,2}}
+
+    # sparsemortar
+    for ax in (
+      ([2, 3], [2, 3]),
+      (([2, 3], [2, 3]),),
+      blockedrange.(([2, 3], [2, 3])),
+      (blockedrange.(([2, 3], [2, 3])),),
+    )
+      blocks = SparseArrayDOK{arrayt{elt,2}}(undef_blocks, ax...)
+      blocks[2, 1] = arrayt(randn(elt, 3, 2))
+      blocks[1, 2] = arrayt(randn(elt, 2, 3))
+      a = sparsemortar(blocks, ax...)
+      @test a isa BlockSparseArray{elt,2,arrayt{elt,2}}
+      @test iszero(a[Block(1, 1)])
+      @test a[Block(2, 1)] == blocks[2, 1]
+      @test a[Block(1, 2)] == blocks[1, 2]
+      @test iszero(a[Block(2, 2)])
+    end
   end
   @testset "Basics" begin
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     @allowscalar @test a == dev(
-      BlockSparseArray{elt}(blockedrange([2, 3]), blockedrange([2, 3]))
+      BlockSparseArray{elt}(undef, blockedrange([2, 3]), blockedrange([2, 3]))
     )
     @test eltype(a) === elt
     @test axes(a) == (1:5, 1:5)
@@ -182,7 +194,7 @@ arrayts = (Array, JLArray)
     @allowscalar @test all(I -> iszero(a[I]), eachindex(a))
     @test_throws DimensionMismatch a[Block(1, 1)] = randn(elt, 2, 3)
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[3, 3] = 33
     @test eltype(a) === elt
     @test axes(a) == (1:5, 1:5)
@@ -201,7 +213,7 @@ arrayts = (Array, JLArray)
       end
     end
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(2, 1)] = randn(elt, 3, 2)
     a[Block(1, 2)] = randn(elt, 2, 3)
     @test issetequal(eachstoredblock(a), [a[Block(2, 1)], a[Block(1, 2)]])
@@ -211,7 +223,7 @@ arrayts = (Array, JLArray)
     @test isnan(norm(a))
 
     # Empty constructor
-    for a in (dev(BlockSparseArray{elt}()), dev(BlockSparseArray{elt}(undef)))
+    for a in (dev(BlockSparseArray{elt}(undef)),)
       @test size(a) == ()
       @test isone(length(a))
       @test blocksize(a) == ()
@@ -247,7 +259,7 @@ arrayts = (Array, JLArray)
     end
 
     @testset "Transpose" begin
-      a = dev(BlockSparseArray{elt}([2, 2], [3, 3, 1]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 2], [3, 3, 1]))
       a[Block(1, 1)] = dev(randn(elt, 2, 3))
       a[Block(2, 3)] = dev(randn(elt, 2, 1))
 
@@ -271,7 +283,7 @@ arrayts = (Array, JLArray)
     end
 
     @testset "Adjoint" begin
-      a = dev(BlockSparseArray{elt}([2, 2], [3, 3, 1]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 2], [3, 3, 1]))
       a[Block(1, 1)] = dev(randn(elt, 2, 3))
       a[Block(2, 3)] = dev(randn(elt, 2, 1))
 
@@ -305,12 +317,12 @@ arrayts = (Array, JLArray)
 
     # TODO: Broken on GPU.
     if arrayt ≠ Array
-      a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
       @test_broken a[Block(1, 2)] .= 2
     end
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 2)] .= 2
     @test eltype(a) == elt
     @test all(==(2), a[Block(1, 2)])
@@ -322,12 +334,12 @@ arrayts = (Array, JLArray)
 
     # TODO: Broken on GPU.
     if arrayt ≠ Array
-      a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+      a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
       @test_broken a[Block(1, 2)] .= 0
     end
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 2)] .= 0
     @test eltype(a) == elt
     @test iszero(a[Block(1, 1)])
@@ -367,7 +379,7 @@ arrayts = (Array, JLArray)
     @test size(b) == size(a)
     @test blocksize(b) == blocksize(a)
 
-    a = dev(BlockSparseArray{elt}([2, 3], [3, 4]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [3, 4]))
     b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
     c = @view b[Block(1, 1)]
     @test iszero(a)
@@ -468,7 +480,7 @@ arrayts = (Array, JLArray)
     @test blockstoredlength(b) == 2
     @test storedlength(b) == 2 * 4 + 3 * 3
 
-    a = dev(BlockSparseArray{elt}([1, 1, 1], [1, 2, 3], [2, 2, 1], [1, 2, 1]))
+    a = dev(BlockSparseArray{elt}(undef, [1, 1, 1], [1, 2, 3], [2, 2, 1], [1, 2, 1]))
     a[Block(3, 2, 2, 3)] = dev(randn(elt, 1, 2, 2, 1))
     perm = (2, 3, 4, 1)
     for b in (PermutedDimsArray(a, perm), permutedims(a, perm))
@@ -623,7 +635,7 @@ arrayts = (Array, JLArray)
     b[Block(1, 1)] .= 1
     @test b[Block(1, 1)] == trues(blocksizes(b)[1, 1])
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[Block(2, 2)]
     @test size(b) == (3, 4)
     for i in parentindices(b)
@@ -632,7 +644,7 @@ arrayts = (Array, JLArray)
     @test parentindices(b)[1] == 1:3
     @test parentindices(b)[2] == 1:4
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[Block(2, 2)[1:2, 2:2]]
     @test size(b) == (2, 1)
     for i in parentindices(b)
@@ -655,7 +667,7 @@ arrayts = (Array, JLArray)
     @test @view(a[Block(2, 2)])[1:1, 1:2] == x
     @test a[3:3, 4:5] == x
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @views a[Block(2, 2)][1:2, 2:3]
     @test b isa SubArray{<:Any,<:Any,<:BlockView}
     for i in parentindices(b)
@@ -667,7 +679,7 @@ arrayts = (Array, JLArray)
     @test a[Block(2, 2)[1:2, 2:3]] == b
     @test blockstoredlength(a) == 1
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     @views for b in [Block(1, 1), Block(2, 2)]
       a[b] = randn(elt, size(a[b]))
     end
@@ -681,7 +693,7 @@ arrayts = (Array, JLArray)
       end
     end
 
-    a = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     @views for b in [Block(1, 1), Block(2, 2)]
       # TODO: Use `blocksizes(a)[Int.(Tuple(b))...]` once available.
       a[b] = dev(randn(elt, size(a[b])))
@@ -768,7 +780,7 @@ arrayts = (Array, JLArray)
     @test !isassigned(a, Block(1)[1], Block(0)[1])
     @test !isassigned(a, Block(3)[3], Block(2)[4])
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     fill!(a, 0)
@@ -782,7 +794,7 @@ arrayts = (Array, JLArray)
     @test iszero(a)
     @test iszero(blockstoredlength(a))
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     zero!(a)
@@ -796,7 +808,7 @@ arrayts = (Array, JLArray)
     @test iszero(a)
     @test iszero(blockstoredlength(a))
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     @test iszero(a)
     @test iszero(blockstoredlength(a))
     a .= 0
@@ -811,7 +823,7 @@ arrayts = (Array, JLArray)
     @test iszero(blockstoredlength(a))
 
     # TODO: Broken on GPU.
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     for I in (Block.(1:2), [Block(1), Block(2)])
       b = @view a[I, I]
       x = randn(elt, 3, 4)
@@ -831,14 +843,14 @@ arrayts = (Array, JLArray)
     end
 
     function f1()
-      a = BlockSparseArray{elt}([2, 3], [3, 4])
+      a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
       b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
       x = randn(elt, 3, 4)
       b[Block(1, 1)] .= x
       return (; a, b, x)
     end
     function f2()
-      a = BlockSparseArray{elt}([2, 3], [3, 4])
+      a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
       b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
       x = randn(elt, 3, 4)
       b[Block(1, 1)] = x
@@ -861,7 +873,7 @@ arrayts = (Array, JLArray)
       @test_throws DimensionMismatch b[Block(1, 1)] .= randn(2, 3)
     end
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @views a[[Block(2), Block(1)], [Block(2), Block(1)]][Block(2, 1)]
     @test iszero(b)
     @test size(b) == (2, 4)
@@ -871,7 +883,7 @@ arrayts = (Array, JLArray)
     @test a[Block(1, 2)] == x
     @test blockstoredlength(a) == 1
 
-    a = BlockSparseArray{elt}([4, 3, 2], [4, 3, 2])
+    a = BlockSparseArray{elt}(undef, [4, 3, 2], [4, 3, 2])
     @views for B in [Block(1, 1), Block(2, 2), Block(3, 3)]
       a[B] = randn(elt, size(a[B]))
     end
@@ -897,13 +909,13 @@ arrayts = (Array, JLArray)
     @test c[Block(2, 2)] == x
     @test a[Block(1, 1)[1:3, 1:3]] == x
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     b = @view a[[Block(2), Block(1)], [Block(2), Block(1)]]
     for index in parentindices(@view(b[Block(1, 1)]))
       @test index isa Base.OneTo{Int}
     end
 
-    a = BlockSparseArray{elt}([2, 3], [3, 4])
+    a = BlockSparseArray{elt}(undef, [2, 3], [3, 4])
     a[Block(1, 1)] = randn(elt, 2, 3)
     b = @view a[Block(1, 1)[1:2, 1:1]]
     @test b isa SubArray{elt,2,Matrix{elt}}
@@ -911,7 +923,7 @@ arrayts = (Array, JLArray)
       @test i isa UnitRange{Int}
     end
 
-    a = BlockSparseArray{elt}([2, 2, 2, 2], [2, 2, 2, 2])
+    a = BlockSparseArray{elt}(undef, [2, 2, 2, 2], [2, 2, 2, 2])
     @views for I in [Block(1, 1), Block(2, 2), Block(3, 3), Block(4, 4)]
       a[I] = randn(elt, size(a[I]))
     end
@@ -932,7 +944,7 @@ arrayts = (Array, JLArray)
       @test c == a[Block.(3:4), Block.(3:4)]
     end
 
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(1, 1)] = randn(elt, 2, 2)
     a[Block(2, 2)] = randn(elt, 3, 3)
     for I in (mortar([Block(1)[2:2], Block(2)[2:3]]), [Block(1)[2:2], Block(2)[2:3]])
@@ -941,7 +953,7 @@ arrayts = (Array, JLArray)
     end
 
     # Merge and permute blocks.
-    a = BlockSparseArray{elt}([2, 2, 2, 2], [2, 2, 2, 2])
+    a = BlockSparseArray{elt}(undef, [2, 2, 2, 2], [2, 2, 2, 2])
     @views for I in [Block(1, 1), Block(2, 2), Block(3, 3), Block(4, 4)]
       a[I] = randn(elt, size(a[I]))
     end
@@ -966,7 +978,7 @@ arrayts = (Array, JLArray)
 
     # TODO: Add more tests of this, it may
     # only be working accidentally.
-    a = BlockSparseArray{elt}([2, 3], [2, 3])
+    a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
     a[Block(1, 1)] = randn(elt, 2, 2)
     a[Block(2, 2)] = randn(elt, 3, 3)
     @test a[2:4, 4] == Array(a)[2:4, 4]
@@ -975,7 +987,7 @@ arrayts = (Array, JLArray)
   end
   @testset "view!" begin
     for blk in ((Block(2, 2),), (Block(2), Block(2)))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = view!(a, blk...)
       x = randn(elt, 3, 3)
       b .= x
@@ -986,7 +998,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2),), (Block(2), Block(2)))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = @view! a[blk...]
       x = randn(elt, 3, 3)
       b .= x
@@ -997,7 +1009,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2)[2:3, 1:2],), (Block(2)[2:3], Block(2)[1:2]))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = view!(a, blk...)
       x = randn(elt, 2, 2)
       b .= x
@@ -1008,7 +1020,7 @@ arrayts = (Array, JLArray)
       @test @view!(a[blk...]) == x
     end
     for blk in ((Block(2, 2)[2:3, 1:2],), (Block(2)[2:3], Block(2)[1:2]))
-      a = BlockSparseArray{elt}([2, 3], [2, 3])
+      a = BlockSparseArray{elt}(undef, [2, 3], [2, 3])
       b = @view! a[blk...]
       x = randn(elt, 2, 2)
       b .= x
@@ -1020,9 +1032,9 @@ arrayts = (Array, JLArray)
     end
   end
   @testset "LinearAlgebra" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
     a_dest = a1 * a2
     @allowscalar @test Array(a_dest) ≈ Array(a1) * Array(a2)
@@ -1030,10 +1042,10 @@ arrayts = (Array, JLArray)
     @test blockstoredlength(a_dest) == 1
   end
   @testset "Matrix multiplication" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 2)] = dev(randn(elt, size(@view(a1[Block(1, 2)]))))
     a1[Block(2, 1)] = dev(randn(elt, size(@view(a1[Block(2, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 2)] = dev(randn(elt, size(@view(a2[Block(1, 2)]))))
     a2[Block(2, 1)] = dev(randn(elt, size(@view(a2[Block(2, 1)]))))
     for (a1′, a2′) in ((a1, a2), (a1', a2), (a1, a2'), (a1', a2'))
@@ -1042,19 +1054,19 @@ arrayts = (Array, JLArray)
     end
   end
   @testset "Dot product" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3, 4]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3, 4]))
     a1[Block(1)] = dev(randn(elt, size(@view(a1[Block(1)]))))
     a1[Block(3)] = dev(randn(elt, size(@view(a1[Block(3)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3, 4]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3, 4]))
     a2[Block(2)] = dev(randn(elt, size(@view(a1[Block(2)]))))
     a2[Block(3)] = dev(randn(elt, size(@view(a1[Block(3)]))))
     @test a1' * a2 ≈ Array(a1)' * Array(a2)
     @test dot(a1, a2) ≈ a1' * a2
   end
   @testset "cat" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(2, 1)] = dev(randn(elt, size(@view(a1[Block(2, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 2)] = dev(randn(elt, size(@view(a2[Block(1, 2)]))))
 
     a_dest = cat(a1, a2; dims=1)
@@ -1079,9 +1091,9 @@ arrayts = (Array, JLArray)
     @test a_dest[Block(3, 4)] == a2[Block(1, 2)]
   end
   @testset "TensorAlgebra" begin
-    a1 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a1 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a1[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
-    a2 = dev(BlockSparseArray{elt}([2, 3], [2, 3]))
+    a2 = dev(BlockSparseArray{elt}(undef, [2, 3], [2, 3]))
     a2[Block(1, 1)] = dev(randn(elt, size(@view(a1[Block(1, 1)]))))
     # TODO: Make this work, requires customization of `TensorAlgebra.fusedims` and
     # `TensorAlgebra.splitdims` in terms of `BlockSparseArrays.blockreshape`,
@@ -1107,7 +1119,7 @@ arrayts = (Array, JLArray)
     matrixt_elt = arrayt{elt,2}
     arrayt_elt = arrayt{elt,3}
 
-    a = BlockSparseVector{elt,arrayt{elt,1}}([2, 2])
+    a = BlockSparseVector{elt,arrayt{elt,1}}(undef, [2, 2])
     res = sprint(summary, a)
     function ref_vec(elt, arrayt, prefix="")
       return "2-blocked 4-element $(prefix)BlockSparseVector{$(elt), $(arrayt), …, …}"
@@ -1116,7 +1128,7 @@ arrayts = (Array, JLArray)
     @test (res == ref_vec(elt, vectort_elt)) ||
       (res == ref_vec(elt, vectort_elt, "BlockSparseArrays."))
 
-    a = BlockSparseMatrix{elt,arrayt{elt,2}}([2, 2], [2, 2])
+    a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [2, 2], [2, 2])
     res = sprint(summary, a)
     function ref_mat(elt, arrayt, prefix="")
       return "2×2-blocked 4×4 $(prefix)BlockSparseMatrix{$(elt), $(arrayt), …, …}"
@@ -1125,7 +1137,7 @@ arrayts = (Array, JLArray)
     @test (res == ref_mat(elt, matrixt_elt)) ||
       (res == ref_mat(elt, matrixt_elt, "BlockSparseArrays."))
 
-    a = BlockSparseArray{elt,3,arrayt{elt,3}}([2, 2], [2, 2], [2, 2])
+    a = BlockSparseArray{elt,3,arrayt{elt,3}}(undef, [2, 2], [2, 2], [2, 2])
     res = sprint(summary, a)
     function ref_arr(elt, arrayt, prefix="")
       return "2×2×2-blocked 4×4×4 $(prefix)BlockSparseArray{$(elt), 3, $(arrayt), …, …}"
@@ -1136,7 +1148,7 @@ arrayts = (Array, JLArray)
     if elt === Float64
       # Not testing other element types since they change the
       # spacing so it isn't easy to make the test general.
-      a = BlockSparseMatrix{elt,arrayt{elt,2}}([2, 2], [2, 2])
+      a = BlockSparseMatrix{elt,arrayt{elt,2}}(undef, [2, 2], [2, 2])
       @allowscalar a[1, 2] = 12
       @test sprint(show, "text/plain", a) ==
         "$(summary(a)):\n $(zero(eltype(a)))  $(eltype(a)(12))  │   ⋅    ⋅ \n $(zero(eltype(a)))   $(zero(eltype(a)))  │   ⋅    ⋅ \n ───────────┼──────────\n  ⋅     ⋅   │   ⋅    ⋅ \n  ⋅     ⋅   │   ⋅    ⋅ "
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 0e0d20ef..7121c797 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -9,6 +9,7 @@ using Test: @test, @testset
     :blockstoredlength,
     :eachblockstoredindex,
     :eachstoredblock,
+    :sparsemortar,
   ]
   @test issetequal(names(BlockSparseArrays), exports)
 end
diff --git a/test/test_gradedunitrangesext.jl b/test/test_gradedunitrangesext.jl
index 4988028c..fbbd0296 100644
--- a/test/test_gradedunitrangesext.jl
+++ b/test/test_gradedunitrangesext.jl
@@ -20,7 +20,7 @@ using TensorAlgebra: fusedims, splitdims
 using LinearAlgebra: adjoint
 using Random: randn!
 function randn_blockdiagonal(elt::Type, axes::Tuple)
-  a = BlockSparseArray{elt}(axes)
+  a = BlockSparseArray{elt}(undef, axes)
   blockdiaglength = minimum(blocksize(a))
   for i in 1:blockdiaglength
     b = Block(ntuple(Returns(i), ndims(a)))
@@ -135,7 +135,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
   @testset "dual axes" begin
     r = gradedrange([U1(0) => 2, U1(1) => 2])
     for ax in ((r, r), (dual(r), r), (r, dual(r)), (dual(r), dual(r)))
-      a = BlockSparseArray{elt}(ax...)
+      a = BlockSparseArray{elt}(undef, ax...)
       @views for b in [Block(1, 1), Block(2, 2)]
         a[b] = randn(elt, size(a[b]))
       end
@@ -178,7 +178,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "GradedOneTo" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -199,7 +199,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "GradedUnitRange" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])[1:3]
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -226,7 +226,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
     # Test case when all axes are dual.
     @testset "dual GradedOneTo" begin
       r = gradedrange([U1(-1) => 2, U1(1) => 2])
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -251,7 +251,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "dual GradedUnitRange" begin
       r = gradedrange([U1(0) => 2, U1(1) => 2])[1:3]
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -277,7 +277,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
 
     @testset "dual BlockedUnitRange" begin    # self dual
       r = blockedrange([2, 2])
-      a = BlockSparseArray{elt}(dual(r), dual(r))
+      a = BlockSparseArray{elt}(undef, dual(r), dual(r))
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -302,7 +302,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
       gradedrange([U1(0) => 2, U1(1) => 2]),
       gradedrange([U1(0) => 2, U1(1) => 2])[begin:end],
     )
-      a = BlockSparseArray{elt}(r, r)
+      a = BlockSparseArray{elt}(undef, r, r)
       @views for i in [Block(1, 1), Block(2, 2)]
         a[i] = randn(elt, size(a[i]))
       end
@@ -330,16 +330,16 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
   end
   @testset "Matrix multiplication" begin
     r = gradedrange([U1(0) => 2, U1(1) => 3])
-    a1 = BlockSparseArray{elt}(dual(r), r)
+    a1 = BlockSparseArray{elt}(undef, dual(r), r)
     a1[Block(1, 2)] = randn(elt, size(@view(a1[Block(1, 2)])))
     a1[Block(2, 1)] = randn(elt, size(@view(a1[Block(2, 1)])))
-    a2 = BlockSparseArray{elt}(dual(r), r)
+    a2 = BlockSparseArray{elt}(undef, dual(r), r)
     a2[Block(1, 2)] = randn(elt, size(@view(a2[Block(1, 2)])))
     a2[Block(2, 1)] = randn(elt, size(@view(a2[Block(2, 1)])))
     @test Array(a1 * a2) ≈ Array(a1) * Array(a2)
     @test Array(a1' * a2') ≈ Array(a1') * Array(a2')
 
-    a2 = BlockSparseArray{elt}(r, dual(r))
+    a2 = BlockSparseArray{elt}(undef, r, dual(r))
     a2[Block(1, 2)] = randn(elt, size(@view(a2[Block(1, 2)])))
     a2[Block(2, 1)] = randn(elt, size(@view(a2[Block(2, 1)])))
     @test Array(a1' * a2) ≈ Array(a1') * Array(a2)
diff --git a/test/test_svd.jl b/test/test_svd.jl
index 8d8dfe79..73efa160 100644
--- a/test/test_svd.jl
+++ b/test/test_svd.jl
@@ -47,7 +47,7 @@ end
 # -----------
 @testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in
                                               Iterators.product(blockszs, eltypes)
-  a = BlockSparseArray{T}(m, n)
+  a = BlockSparseArray{T}(undef, m, n)
 
   # test empty matrix
   usv_empty = svd(a)
diff --git a/test/test_tensoralgebraext.jl b/test/test_tensoralgebraext.jl
index 470efaad..3fe0b869 100644
--- a/test/test_tensoralgebraext.jl
+++ b/test/test_tensoralgebraext.jl
@@ -9,7 +9,7 @@ using SymmetrySectors: U1
 using TensorAlgebra: contract
 
 function randn_blockdiagonal(elt::Type, axes::Tuple)
-  a = BlockSparseArray{elt}(axes)
+  a = BlockSparseArray{elt}(undef, axes)
   blockdiaglength = minimum(blocksize(a))
   for i in 1:blockdiaglength
     b = Block(ntuple(Returns(i), ndims(a)))