diff --git a/.gitignore b/.gitignore
index ca28dd4d..3d39b6dc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,7 +4,8 @@
 docs/build/
 docs/site/
 *.jld
-benchmark/*.md
+benchmark/*.json
+benchmark/results_*.md
 src/.DS_Store
 Manifest.toml
 Manifest-v*.*.toml
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 00000000..f740db19
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,9 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[compat]
+BenchmarkTools = "1"
+julia = "1.10"
diff --git a/benchmark/README.md b/benchmark/README.md
new file mode 100644
index 00000000..7cdff9e0
--- /dev/null
+++ b/benchmark/README.md
@@ -0,0 +1,23 @@
+# Benchmarks
+
+Set up the benchmark environment from the repository root:
+
+```julia
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+```
+
+Run and save the suite under a descriptive name:
+
+```julia
+julia --project=benchmark -L benchmark/runbenchmarks.jl -e 'run_benchmarks("main")'
+```
+
+Generate a standalone report or compare two saved runs:
+
+```julia
+julia --project=benchmark -L benchmark/runbenchmarks.jl -e 'generate_report("main")'
+julia --project=benchmark -L benchmark/runbenchmarks.jl -e 'generate_report("main", "branch")'
+```
+
+Benchmark parameters and results are stored as ignored JSON files under
+`benchmark/`.
diff --git a/benchmark/generate_report.jl b/benchmark/generate_report.jl
index be976b93..5786f4a6 100644
--- a/benchmark/generate_report.jl
+++ b/benchmark/generate_report.jl
@@ -40,7 +40,11 @@ function printreport(io::IO, results; iscomparisonjob::Bool = false)
     return nothing
 end
 
-idrepr(id) = (str = repr(id); str[searchindex(str, '['):end])
+function idrepr(id)
+    str = repr(id)
+    firstbracket = findfirst(==('['), str)
+    return isnothing(firstbracket) ? str : str[firstbracket:end]
+end
 intpercent(p) = string(ceil(Int, p * 100), "%")
 resultrow(ids, t::BenchmarkTools.Trial) = resultrow(ids, minimum(t))
 
diff --git a/benchmark/runbenchmarks.jl b/benchmark/runbenchmarks.jl
index bfdaa202..55d3e3d7 100644
--- a/benchmark/runbenchmarks.jl
+++ b/benchmark/runbenchmarks.jl
@@ -1,7 +1,6 @@
 using BlockArrays
 using BenchmarkTools
-using FileIO
-using JLD
+using LinearAlgebra
 
 include("generate_report.jl")
 
@@ -10,6 +9,10 @@ const SUITE = BenchmarkGroup()
 g = addgroup!(SUITE, "indexing")
 # g_block = addgroup!(SUITE, "blockindexing")
 g_size = addgroup!(SUITE, "size")
+g_metadata = addgroup!(SUITE, "metadata")
+g_product = addgroup!(SUITE, "product")
+g_broadcast = addgroup!(SUITE, "broadcast")
+g_reduction = addgroup!(SUITE, "reduction")
 
 for n = (5,)
     for BT in (BlockArray, BlockedArray)
@@ -30,31 +33,68 @@ for n = (5,)
     end
 end
 
+blockkron_vector = Ref(BlockKron(1:10_000, 1:4, 1:3))
+blockkron_matrix = Ref(BlockKron(
+    reshape(1:800_000, 1_000, 800), reshape(1:20, 4, 5), reshape(1:6, 3, 2),
+))
+g_metadata["BlockKron", "axes", "vector"] = @benchmarkable axes($blockkron_vector[])
+g_metadata["BlockKron", "axes", "matrix"] = @benchmarkable axes($blockkron_matrix[])
+
+blocklasts_a = collect(2:2:20_000)
+blocklasts_b = collect(3:3:30_000)
+g_metadata["sortedunion", "vectors"] =
+    @benchmarkable BlockArrays.sortedunion($blocklasts_a, $blocklasts_b)
+
+khatri_block_sizes = fill(2, 10)
+khatri_a = BlockArray(randn(20, 20), khatri_block_sizes, khatri_block_sizes)
+khatri_b = BlockArray(randn(20, 20), khatri_block_sizes, khatri_block_sizes)
+g_product["khatri_rao", "10x10", "2x2 blocks"] =
+    @benchmarkable khatri_rao($khatri_a, $khatri_b)
+
+broadcast_block_sizes = fill(4, 32)
+broadcast_blocked_a = BlockedArray(randn(128, 128), broadcast_block_sizes, broadcast_block_sizes)
+broadcast_blocked_b = BlockedArray(randn(128, 128), broadcast_block_sizes, broadcast_block_sizes)
+broadcast_blocked_dest = similar(broadcast_blocked_a)
+g_broadcast["in-place", "BlockedArray", "matching"] =
+    @benchmarkable $broadcast_blocked_dest .= $broadcast_blocked_a .+ $broadcast_blocked_b
+
+broadcast_block_a = BlockArray(randn(128, 128), broadcast_block_sizes, broadcast_block_sizes)
+broadcast_block_b = BlockArray(randn(128, 128), broadcast_block_sizes, broadcast_block_sizes)
+broadcast_block_dest = similar(broadcast_block_a)
+g_broadcast["in-place", "BlockArray", "matching"] =
+    @benchmarkable $broadcast_block_dest .= $broadcast_block_a .+ $broadcast_block_b
+
+reduction_block_sizes = fill(4, 32)
+reduction_block_array = BlockArray(randn(128, 128), reduction_block_sizes, reduction_block_sizes)
+g_reduction["sum", "BlockArray", "matrix"] = @benchmarkable sum($reduction_block_array)
+g_reduction["sum(abs2)", "BlockArray", "matrix"] =
+    @benchmarkable sum(abs2, $reduction_block_array)
+g_reduction["norm", "BlockArray", "matrix"] = @benchmarkable norm($reduction_block_array)
+
 
 function run_benchmarks(name, tagfilter = @tagged ALL)
-    paramspath = joinpath(dirname(@__FILE__), "params.jld")
+    paramspath = joinpath(@__DIR__, "params.json")
     if !isfile(paramspath)
         println("Tuning benchmarks...")
         tune!(SUITE, verbose=true)
-        JLD.save(paramspath, "SUITE", params(SUITE))
+        BenchmarkTools.save(paramspath, params(SUITE))
     end
-    loadparams!(SUITE, JLD.load(paramspath, "SUITE"), :evals, :samples)
+    loadparams!(SUITE, only(BenchmarkTools.load(paramspath)), :evals, :samples)
     results = run(SUITE[tagfilter], verbose = true, seconds = 2)
-    JLD.save(joinpath(dirname(@__FILE__), name * ".jld"), "results", results)
+    BenchmarkTools.save(joinpath(@__DIR__, name * ".json"), results)
 end
 
 function generate_report(v1, v2)
-    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
-    v2_res = load(joinpath(dirname(@__FILE__), v2 * ".jld"), "results")
-    open(joinpath(dirname(@__FILE__), "results_$(v1)_$(v2).md"), "w") do f
+    v1_res = only(BenchmarkTools.load(joinpath(@__DIR__, v1 * ".json")))
+    v2_res = only(BenchmarkTools.load(joinpath(@__DIR__, v2 * ".json")))
+    open(joinpath(@__DIR__, "results_$(v1)_$(v2).md"), "w") do f
         printreport(f, judge(minimum(v1_res), minimum(v2_res)); iscomparisonjob = true)
     end
 end
 
 function generate_report(v1)
-    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
-    open(joinpath(dirname(@__FILE__), "results_$(v1).md"), "w") do f
+    v1_res = only(BenchmarkTools.load(joinpath(@__DIR__, v1 * ".json")))
+    open(joinpath(@__DIR__, "results_$(v1).md"), "w") do f
         printreport(f, minimum(v1_res); iscomparisonjob = false)
     end
 end
-
diff --git a/src/blockaxis.jl b/src/blockaxis.jl
index 1acfd5d4..c71ae4dd 100644
--- a/src/blockaxis.jl
+++ b/src/blockaxis.jl
@@ -735,6 +735,34 @@ Base.BroadcastStyle(::Type{<:AbstractBlockedUnitRange{<:Any,R}}) where R = _broa
 ###
 
 const OneToCumsum{T<:Integer} = RangeCumsum{T,Base.OneTo{T}}
+
+# Base has no public vector merge that exploits sorted inputs, so merge block
+# boundaries directly instead of hashing and then sorting them with `union`.
+function sortedunion(a::StridedVector{<:Integer}, b::StridedVector{<:Integer})
+    T = promote_type(eltype(a), eltype(b))
+    result = Vector{T}(undef, length(a) + length(b))
+    ia, ib = firstindex(a), firstindex(b)
+    lasta, lastb = lastindex(a), lastindex(b)
+    nresult = 0
+
+    @inbounds while ia <= lasta || ib <= lastb
+        value = if ib > lastb || (ia <= lasta && !isless(b[ib], a[ia]))
+            value = a[ia]
+            ia += 1
+            value
+        else
+            value = b[ib]
+            ib += 1
+            value
+        end
+        if iszero(nresult) || !isequal(result[nresult], value)
+            nresult += 1
+            result[nresult] = value
+        end
+    end
+    return resize!(result, nresult)
+end
+
 sortedunion(a::OneToCumsum, ::OneToCumsum) = a
 function sortedunion(a::RangeCumsum{<:Any,<:AbstractRange}, b::RangeCumsum{<:Any,<:AbstractRange})
     @assert a == b
diff --git a/src/blockbroadcast.jl b/src/blockbroadcast.jl
index c1180268..bb74900d 100644
--- a/src/blockbroadcast.jl
+++ b/src/blockbroadcast.jl
@@ -33,7 +33,6 @@ BroadcastStyle(::BlockStyle{M}, ::BlockedStyle{N}) where {M,N} = BlockStyle(Val(
 BroadcastStyle(::BlockedStyle{M}, ::BlockStyle{N}) where {M,N} = BlockStyle(Val(max(M,N)))
 
 
-# sortedunion can assume inputs are already sorted so this could be improved
 maybeinplacesort!(v::StridedVector) = sort!(v)
 maybeinplacesort!(v) = sort(v)
 sortedunion(a,b) = maybeinplacesort!(union(a,b))
@@ -42,6 +41,16 @@ sortedunion(a::AbstractUnitRange, b::AbstractUnitRange) = min(first(a),first(b))
 combine_blockaxes(a, b) = _BlockedUnitRange(sortedunion(blocklasts(a), blocklasts(b)))
 combine_blockaxes(a::BlockedOneTo, b::BlockedOneTo) = BlockedOneTo(sortedunion(blocklasts(a), blocklasts(b)))
 
+function combine_blockaxes(a::BlockedUnitRange{T,Vector{T}},
+                           b::BlockedUnitRange{T,Vector{T}}) where T
+    return blockisequal(a, b) ? a : _BlockedUnitRange(sortedunion(blocklasts(a), blocklasts(b)))
+end
+
+function combine_blockaxes(a::BlockedOneTo{T,CS}, b::BlockedOneTo{T,CS}) where
+                           {T<:Integer,CS<:Union{Vector{T},Base.OneTo{T}}}
+    return blockisequal(a, b) ? a : BlockedOneTo(sortedunion(blocklasts(a), blocklasts(b)))
+end
+
 Base.Broadcast.axistype(a::AbstractBlockedUnitRange, b::AbstractBlockedUnitRange) = length(b) == 1 ? a : combine_blockaxes(a, b)
 Base.Broadcast.axistype(a::AbstractBlockedUnitRange, b) = length(b) == 1 ? a : combine_blockaxes(a, b)
 Base.Broadcast.axistype(a, b::AbstractBlockedUnitRange) = length(b) == 1 ? a : combine_blockaxes(a, b)
@@ -176,31 +185,26 @@ copyto!(dest::AbstractArray,
     _generic_blockbroadcast_copyto!(dest, bc)
 
 # type-stable version of _bview.(args, K)
-__bview(args::Tuple{}, K) = ()
-__bview(args::Tuple, K) = tuple(_bview(args[1],K), __bview(tail(args), K)...)
+__bview(args::Tuple{}, K...) = ()
+__bview(args::Tuple, K...) = tuple(_bview(args[1], K...), __bview(tail(args), K...)...)
 
-function _fast_blockbradcast_copyto!(dest, bc)
-    @inbounds for K in blockaxes(bc)[1]
-        broadcast!(bc.f, view(dest,K), __bview(bc.args, K)...)
+function _fast_blockbroadcast_copyto!(dest, bc)
+    @inbounds for K in Iterators.product(blockaxes(bc)...)
+        broadcast!(bc.f, view(dest, K...), __bview(bc.args, K...)...)
     end
     dest
 end
 
-_hasscalarlikevec() = false
-_hasscalarlikevec(a, b...) = _hasscalarlikevec(b...)
-_hasscalarlikevec(a::AbstractVector, b...) = size(a,1) == 1 || _hasscalarlikevec(b...)
-
-blockisequalorscalar(ax, ::Number) = true
-blockisequalorscalar(ax, a) = blockisequal(ax, Base.axes1(a))
+blockisequalorscalar(ax::Tuple, ::Number) = true
+blockisequalorscalar(ax::Tuple, a) = blockisequal(ax, axes(a))
 
-function copyto!(dest::AbstractVector,
-        bc::Broadcasted{<:AbstractBlockStyle{1}, <:Any, <:Any, Args}) where {Args <: Tuple}
-    _hasscalarlikevec(bc.args...) && return _generic_blockbroadcast_copyto!(dest, bc)
-    ax = axes(dest,1)
+function copyto!(dest::AbstractArray{<:Any,N},
+        bc::Broadcasted{<:AbstractBlockStyle{N},<:Any,<:Any,Args}) where {N,Args<:Tuple}
+    ax = axes(dest)
     for a in bc.args
         blockisequalorscalar(ax, a) || return _generic_blockbroadcast_copyto!(dest, bc)
     end
-    return _fast_blockbradcast_copyto!(dest, bc)
+    return _fast_blockbroadcast_copyto!(dest, bc)
 end
 @inline function Broadcast.instantiate(bc::Broadcasted{Style}) where {Style <:BlockStyle}
     bcf = Broadcast.instantiate(Broadcast.flatten(Broadcasted{Nothing}(bc.f, bc.args, bc.axes)))
@@ -214,6 +218,13 @@ _removeblocks(a::Adjoint) = _removeblocks(parent(a))'
 _removeblocks(a::Transpose) = transpose(_removeblocks(parent(a)))
 _removeblocks(a::SubArray{<:Any,N,<:BlockedArray}) where N = view(_removeblocks(parent(a)), map(_removeblocks, parentindices(a))...)
 _removeblocks(a) = a
+
+function copyto!(dest::BlockedArray{<:Any,N},
+        bc::Broadcasted{BlockedStyle{N},<:Any,<:Any,Args}) where {N,Args<:Tuple}
+    copyto!(parent(dest), _removeblocks(bc))
+    return dest
+end
+
 copy(bc::Broadcasted{BlockedStyle{N}}) where N = BlockedArray(Broadcast.materialize(_removeblocks(bc)), axes(bc))
 
 for op in (:+, :-, :*)
diff --git a/src/blockproduct.jl b/src/blockproduct.jl
index 19504822..baa982b5 100644
--- a/src/blockproduct.jl
+++ b/src/blockproduct.jl
@@ -6,23 +6,14 @@ References
 * Khatri, C. G., and Rao, C. Radhakrishna (1968) Solutions to Some Functional Equations and Their Applications to Characterization of Probability Distributions. Sankhya: Indian J. Statistics, Series A 30, 167–180.
 """
 function khatri_rao(A::AbstractBlockMatrix, B::AbstractBlockMatrix)
-    #
-    Ablksize = blocksize(A)
-    Bblksize = blocksize(B)
-
-    @assert Ablksize == Bblksize "A and B must have the same blocksize"
-
-    kblk = []
-    for iblk in blockaxes(A,1)
-        kblk_j = []
-        for _jblk in blockaxes(A,2)
-            Ablk = A[iblk, _jblk]
-            Bblk = B[iblk, _jblk]
-            push!(kblk_j, kron(Ablk, Bblk))
-        end
-        push!(kblk, tuple(kblk_j...))
+    @assert blocksize(A) == blocksize(B) "A and B must have the same blocksize"
+
+    product = Iterators.product(blockaxes(A)...)
+    result_blocks = map(product) do block_index
+        K, J = block_index
+        kron(view(A, K, J), view(B, K, J))
     end
-    mortar(kblk...)
+    return mortar(result_blocks)
 end
 
 function khatri_rao(A::AbstractMatrix, B::AbstractMatrix)
@@ -51,15 +42,11 @@ size(K::BlockKron, j::Int) = prod(size.(K.args, j))
 size(a::BlockKron{<:Any,1}) = (size(a,1),)
 size(a::BlockKron{<:Any,2}) = (size(a,1), size(a,2))
 
-function axes(K::BlockKron{<:Any,1})
-    A,B = K.args
-    (blockedrange(fill(prod(size.(tail(K.args),1)), size(K.args[1],1))),)
-end
-
-function axes(K::BlockKron{<:Any,2})
-    A,B = K.args
-    blockedrange.((fill(prod(size.(tail(K.args),1)), size(K.args[1],1)),
-                   fill(prod(size.(tail(K.args),2)), size(K.args[1],2))))
+function axes(K::BlockKron{<:Any,N}) where N
+    ntuple(Val(N)) do dim
+        blocklength = prod(size.(tail(K.args), dim))
+        blockedrange(Fill(blocklength, size(K.args[1], dim)))
+    end
 end
 
 kron_getindex((A,)::Tuple{AbstractVector}, k::Integer) = A[k]
diff --git a/src/blockreduce.jl b/src/blockreduce.jl
index f1c55363..873a662a 100644
--- a/src/blockreduce.jl
+++ b/src/blockreduce.jl
@@ -17,5 +17,25 @@ Base.mapfoldl(f::F, op::OP, B::BlockedArray; kw...) where {F, OP} =
 Base.mapreduce(f::F, op::OP, B::BlockedArray; kw...) where {F, OP} =
     mapreduce(f, op, B.blocks; kw...)
 
+function Base.mapreduce(f::F, op::typeof(Base.add_sum), B::BlockArray;
+                        dims=:, kw...) where F
+    if dims isa Colon && !isempty(B)
+        nonemptyblocks = Iterators.filter(!isempty, B.blocks)
+        return mapreduce(block -> mapreduce(f, op, block), op, nonemptyblocks; kw...)
+    end
+    return invoke(mapreduce, Tuple{Any,Any,AbstractArray}, f, op, B; dims, kw...)
+end
+Base.mapreduce(f::F, op::typeof(Base.add_sum), B::BlockVector;
+               dims=:, kw...) where F =
+    invoke(mapreduce, Tuple{F,typeof(op),BlockArray}, f, op, B; dims, kw...)
+
+_norm_hypot(x, y) = (ismissing(x) || isnan(x)) ? x :
+                    (ismissing(y) || isnan(y)) ? y : hypot(x, y)
+
+function LinearAlgebra.norm2(B::BlockArray)
+    isempty(B.blocks) && return float(norm(zero(eltype(B))))
+    return mapreduce(norm, _norm_hypot, B.blocks)
+end
+
 # support sum, need to return something analogous to Base.OneTo(1) but same type
 Base.reduced_index(::BR) where BR<:AbstractBlockedUnitRange = convert(BR, Base.OneTo(1))
diff --git a/test/test_blockbroadcast.jl b/test/test_blockbroadcast.jl
index cbfdd97c..307ea258 100644
--- a/test/test_blockbroadcast.jl
+++ b/test/test_blockbroadcast.jl
@@ -28,6 +28,15 @@ using StaticArrays
         @test axes(A + A) == axes(A .+ A) == axes(A)
         @test axes(A .+ 1) == axes(A)
 
+        dest = similar(A)
+        @test (dest .= A .+ 2 .* A) === dest
+        @test Matrix(dest) ≈ Matrix(A) .+ 2 .* Matrix(A)
+
+        A3 = BlockArray(randn(4, 4, 4), [2, 2], [1, 3], [3, 1])
+        dest3 = similar(A3)
+        @test (dest3 .= A3 .+ A3) === dest3
+        @test Array(dest3) ≈ 2 .* Array(A3)
+
         @testset "mismatched ndims" begin
             u = BlockArray(randn(5), [2,3])
             dest = zeros(size(u)..., 1)
@@ -81,6 +90,19 @@ using StaticArrays
             @test dest ≈ x + 2y
         end
 
+        @testset "matrix in-place broadcast" begin
+            x = BlockedMatrix(randn(6, 6), [2, 4], [3, 3])
+            y = BlockedMatrix(randn(6, 6), [1, 2, 3], [2, 4])
+            dest = similar(x)
+            @test (dest .= x .+ 2 .* y) === dest
+            @test parent(dest) ≈ parent(x) .+ 2 .* parent(y)
+            @test axes(x .+ x) === axes(x)
+
+            expected = parent(x) .+ parent(y)
+            x .+= y
+            @test parent(x) ≈ expected
+        end
+
         @testset "0-dim nested in-place broadcast" begin
             x = BlockedArray(randn(()))
             y = BlockedArray(randn(()))
@@ -127,6 +149,17 @@ using StaticArrays
         B = BlockArray(randn(6,6), fill(2,3), fill(3,2))
 
         @test blocksize(A+B) == (5,3)
+
+        dest = similar(A)
+        @test (dest .= A .+ B) === dest
+        @test Matrix(dest) ≈ Matrix(A) .+ Matrix(B)
+    end
+
+    @testset "sorted block boundary union" begin
+        @test BlockArrays.sortedunion([1, 3, 3, 7], [2, 3, 8]) == [1, 2, 3, 7, 8]
+        @test BlockArrays.sortedunion(Int[], Int[]) == Int[]
+        @test BlockArrays.sortedunion(Int32[1, 3], Int64[2, 3]) == [1, 2, 3]
+        @test eltype(BlockArrays.sortedunion(Int32[1], Int64[2])) == Int64
     end
 
     @testset "UnitRange" begin
@@ -169,10 +202,18 @@ using StaticArrays
     @testset "special axes" begin
         A = BlockArray(randn(6), Ones{Int}(6))
         B = BlockArray(randn(6), Ones{Int}(6))
+        @test (@inferred BlockArrays.combine_blockaxes(axes(A, 1), axes(B, 1))) === axes(A, 1)
         @test axes(A+B,1) === axes(A,1)
 
         C = BlockArray(randn(6), (BlockArrays._BlockedUnitRange(1,2:6),))
         @test axes(A+C,1) === BlockArrays._BlockedUnitRange(1,1:6)
+
+        a32 = blockedrange(Int32[1, 2, 3])
+        a64 = blockedrange(Int64[1, 2, 3])
+        @test blockisequal(@inferred(BlockArrays.combine_blockaxes(a32, a64)), a64)
+
+        lazy = blockedrange(Fill(4, 32))
+        @test blockisequal(@inferred(BlockArrays.combine_blockaxes(lazy, lazy)), lazy)
     end
 
     @testset "Views" begin
diff --git a/test/test_blocklinalg.jl b/test/test_blocklinalg.jl
index d10f0d70..23cbe060 100644
--- a/test/test_blocklinalg.jl
+++ b/test/test_blocklinalg.jl
@@ -35,6 +35,28 @@ bview(a, b) = Base.invoke(view, Tuple{AbstractArray,Any}, a, b)
     @test a .^ 2 == 4
     end
 
+    @testset "norm" begin
+        data = randn(7, 9)
+        A = BlockArray(data, [1, 2, 4], [3, 1, 5])
+        @test @inferred(norm(A)) ≈ norm(data)
+        @test @inferred(LinearAlgebra.norm2(A)) ≈ norm(data)
+        @test norm(A, 1) ≈ norm(data, 1)
+        @test norm(A, 3) ≈ norm(data, 3)
+
+        extremes = BlockArray([1.0e300 1.0e-300], [1], [1, 1])
+        @test norm(extremes) == norm(Matrix(extremes)) == 1.0e300
+
+        nonfinite = BlockArray([NaN, Inf], [1, 1])
+        @test isnan(norm(nonfinite)) == isnan(norm(Vector(nonfinite)))
+
+        withmissing = BlockArray(Union{Missing,Float64}[1.0, missing], [1, 1])
+        @test norm(withmissing) === norm(Vector(withmissing)) === missing
+
+        emptyblocks = BlockArray(zeros(0, 0), Int[], Int[])
+        @test norm(emptyblocks) === 0.0
+        @test LinearAlgebra.norm2(emptyblocks) === 0.0
+    end
+
     @testset "BlockArray scalar * matrix" begin
         A = BlockArray{Float64}(randn(6,6), fill(2,3), 1:3)
         @test 2A == A*2 == 2Matrix(A)
diff --git a/test/test_blockproduct.jl b/test/test_blockproduct.jl
index 13203f11..0dd0286c 100644
--- a/test/test_blockproduct.jl
+++ b/test/test_blockproduct.jl
@@ -27,10 +27,11 @@ using BlockArrays, Test
         A = BlockArray(ones(m, n), mi, ni)
         B = BlockArray(ones(p, q), pi, qi)
 
-        AB = khatri_rao(A, B)
+        AB = @inferred khatri_rao(A, B)
 
         @test blocksize(AB) == blocksize(A)
         @test blocksize(AB) == blocksize(B)
+        @test khatri_rao(BlockedArray(A), BlockedArray(B)) == AB
 
         #Test: Size of resulting blocks
         for i in blockaxes(AB,1)
@@ -116,6 +117,8 @@ using BlockArrays, Test
         c = 6:8
         k̄ = BlockKron(a,b,c)
         @test k̄ == blockkron(a,b,c) == kron(a,b,c)
+        @test @inferred(axes(k̄)) == (blockedrange(fill(length(b) * length(c), length(a))),)
+        @test blocklasts(axes(k̄, 1)) isa FirstStepRangeLen
         @test k̄[Block(1)][Block(1)] == a[1]*b[1]*c
         @test k̄[Block(1)][Block(2)] == a[1]*b[2]*c
         @test k̄[Block(2)][Block(3)] == a[2]*b[3]*c
@@ -133,6 +136,11 @@ using BlockArrays, Test
         C = randn(2,5)
         K̄ = BlockKron(A,B,C)
         @test K̄ ≈ blockkron(A,B,C) ≈ kron(A,B,C)
+        @test @inferred(axes(K̄)) == (
+            blockedrange(fill(size(B, 1) * size(C, 1), size(A, 1))),
+            blockedrange(fill(size(B, 2) * size(C, 2), size(A, 2))),
+        )
+        @test all(ax -> blocklasts(ax) isa FirstStepRangeLen, axes(K̄))
         @test K̄[Block(1,1)][Block(1,1)] ≈ A[1,1]*B[1,1]*C
         @test K̄[Block(2,3)][Block(3,4)] ≈ A[2,3]*B[3,4]*C
 
diff --git a/test/test_blockreduce.jl b/test/test_blockreduce.jl
index 5af3afab..56b1fd7a 100644
--- a/test/test_blockreduce.jl
+++ b/test/test_blockreduce.jl
@@ -2,6 +2,11 @@ module TestBlockReduce
 
 using BlockArrays, Test
 
+struct NoZero
+    value::Int
+end
+Base.:+(a::NoZero, b::NoZero) = NoZero(a.value + b.value)
+
 @testset "foldl" begin
     x = mortar([rand(3), rand(2)])
     @test foldl(push!, x; init = []) == collect(x)
@@ -21,11 +26,38 @@ end
 @testset "sum (#141)" begin
     data = reshape(collect(1:20), 4, 5)
     A = BlockArray(data, [1,3], [2,3])
-    @test sum(A) == sum(data)
+    @test @inferred(sum(A)) == sum(data)
+    @test @inferred(mapreduce(identity, Base.add_sum, A)) == sum(data)
+    @test mapreduce(identity, Base.add_sum, A; init=10) == sum(data; init=10)
+    @test mapreduce(identity, Base.add_sum, A; dims=1) == sum(data; dims=1)
+    @test @inferred(sum(abs2, A)) == sum(abs2, data)
+    @test @inferred(mapreduce(abs2, Base.add_sum, A)) == sum(abs2, data)
+    @test mapreduce(abs2, Base.add_sum, A; init=10) == sum(abs2, data; init=10)
+    @test mapreduce(abs2, Base.add_sum, A; dims=2) == sum(abs2, data; dims=2)
+    @test sum(A; init=10) == sum(data; init=10)
     @test sum(A; dims=1) == sum(data; dims=1)
     @test sum(A; dims=2) == sum(data; dims=2)
     @test blockisequal(axes(A,2), axes(sum(A; dims=1),2))
     @test blockisequal(axes(A,1), axes(sum(A; dims=2),1))
+
+    smallints = BlockArray(fill(Int8(1), 4, 4), [2, 2], [1, 3])
+    @test sum(smallints) === sum(Matrix(smallints)) === 16
+
+    emptyblocks = BlockArray(zeros(0, 0), Int[], Int[])
+    @test sum(emptyblocks) === 0.0
+    @test sum(abs2, emptyblocks) === 0.0
+    @test sum(emptyblocks; init=10.0) === 10.0
+
+    zeroblocks = BlockArray(collect(1:6), [0, 2, 0, 4])
+    @test sum(zeroblocks) == sum(Vector(zeroblocks))
+    @test sum(abs2, zeroblocks) == sum(abs2, Vector(zeroblocks))
+    @test sum(NoZero, zeroblocks).value == sum(NoZero, Vector(zeroblocks)).value
+
+    v = BlockArray(collect(1:6), [2, 4])
+    @test @inferred(sum(v)) == 21
+    @test @inferred(mapreduce(identity, Base.add_sum, v)) == 21
+    @test @inferred(sum(abs2, v)) == 91
+    @test @inferred(mapreduce(abs2, Base.add_sum, v)) == 91
 end
 
 end # module