diff --git a/Project.toml b/Project.toml
index 8a5f85c..28d36e7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "KroneckerArrays"
 uuid = "05d0b138-81bc-4ff7-84be-08becefb1ccc"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.30"
+version = "0.1.31"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -27,14 +27,14 @@ KroneckerArraysTensorProductsExt = "TensorProducts"
 [compat]
 Adapt = "4.3"
 BlockArrays = "1.6"
-BlockSparseArrays = "0.9"
+BlockSparseArrays = "0.9, 0.10.3"
 DerivableInterfaces = "0.5.3"
 DiagonalArrays = "0.3.11"
 FillArrays = "1.13"
 GPUArraysCore = "0.2"
 LinearAlgebra = "1.10"
 MapBroadcast = "0.1.10"
-MatrixAlgebraKit = "0.2"
+MatrixAlgebraKit = "0.2, 0.3"
 TensorAlgebra = "0.3.10"
 TensorProducts = "0.1.7"
 julia = "1.10"
diff --git a/ext/KroneckerArraysBlockSparseArraysExt/KroneckerArraysBlockSparseArraysExt.jl b/ext/KroneckerArraysBlockSparseArraysExt/KroneckerArraysBlockSparseArraysExt.jl
index a28d98d..7a3f129 100644
--- a/ext/KroneckerArraysBlockSparseArraysExt/KroneckerArraysBlockSparseArraysExt.jl
+++ b/ext/KroneckerArraysBlockSparseArraysExt/KroneckerArraysBlockSparseArraysExt.jl
@@ -23,6 +23,17 @@ function BlockSparseArrays.blockrange(bs::Vector{<:CartesianProduct})
   return blockrange(map(cartesianrange, bs))
 end
 
+using BlockArrays: BlockArrays, mortar
+using BlockSparseArrays: blockrange
+using KroneckerArrays: CartesianProductUnitRange
+# Makes sure that `mortar` results in a `BlockVector` with the correct
+# axes, otherwise the axes would not preserve the Kronecker structure.
+# This is helpful when indexing `BlockUnitRange`, for example:
+# https://github.com/JuliaArrays/BlockArrays.jl/blob/v1.7.1/src/blockaxis.jl#L540-L547
+function BlockArrays.mortar(blocks::AbstractVector{<:CartesianProductUnitRange})
+  return mortar(blocks, (blockrange(map(Base.axes1, blocks)),))
+end
+
 using BlockArrays: AbstractBlockedUnitRange
 using BlockSparseArrays: Block, ZeroBlocks, eachblockaxis, mortar_axis
 using KroneckerArrays: KroneckerArrays, KroneckerArray, ⊗, arg1, arg2, _similar
diff --git a/test/Project.toml b/test/Project.toml
index f37978a..2ea9eec 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -23,7 +23,7 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 Adapt = "4"
 Aqua = "0.8"
 BlockArrays = "1.6"
-BlockSparseArrays = "0.9"
+BlockSparseArrays = "0.9, 0.10"
 DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3.7"
 FillArrays = "1"
@@ -31,7 +31,7 @@ GPUArraysCore = "0.2"
 JLArrays = "0.2"
 KroneckerArrays = "0.1"
 LinearAlgebra = "1.10"
-MatrixAlgebraKit = "0.2"
+MatrixAlgebraKit = "0.2, 0.3"
 SafeTestsets = "0.1"
 StableRNGs = "1.0"
 Suppressor = "0.2"
diff --git a/test/test_blocksparsearrays.jl b/test/test_blocksparsearrays.jl
index db492db..bc27327 100644
--- a/test/test_blocksparsearrays.jl
+++ b/test/test_blocksparsearrays.jl
@@ -32,6 +32,13 @@ arrayts = (Array, JLArray)
   @test blockisequal(arg1(r), blockedrange([2, 3]))
   @test blockisequal(arg2(r), blockedrange([3, 4]))
 
+  r = blockrange([2 × 3, 3 × 4])
+  r′ = r[Block.([2, 1])]
+  @test r′[Block(1)] ≡ cartesianrange(3 × 4, 7:18)
+  @test r′[Block(2)] ≡ cartesianrange(2 × 3, 1:6)
+  @test eachblockaxis(r′)[1] ≡ cartesianrange(3, 4)
+  @test eachblockaxis(r′)[2] ≡ cartesianrange(2, 3)
+
   dev = adapt(arrayt)
   r = blockrange([2 × 2, 3 × 3])
   d = Dict(
@@ -137,13 +144,8 @@ arrayts = (Array, JLArray)
     @test_broken inv(a)
   end
 
-  if arrayt === Array
-    u, s, v = svd_compact(a)
-    @test Array(u * s * v) ≈ Array(a)
-  else
-    # Broken on GPU.
-    @test_broken svd_compact(a)
-  end
+  u, s, v = svd_compact(a)
+  @test Array(u * s * v) ≈ Array(a)
 
   b = a[Block.(1:2), Block(2)]
   @test b[Block(1)] == a[Block(1, 2)]
@@ -236,30 +238,66 @@ end
   @test_broken copy(b)
   @test_broken b[Block(1, 2)]
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   b = @constinferred a * a
   @test typeof(b) === typeof(a)
   @test Array(b) ≈ Array(a) * Array(a)
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   # Type inference is broken for this operation.
   # b = @constinferred a + a
   b = a + a
   @test typeof(b) === typeof(a)
   @test Array(b) ≈ Array(a) + Array(a)
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   # Type inference is broken for this operation.
   # b = @constinferred 3a
   b = 3a
   @test typeof(b) === typeof(a)
   @test Array(b) ≈ 3Array(a)
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   # Type inference is broken for this operation.
   # b = @constinferred a / 3
   b = a / 3
   @test typeof(b) === typeof(a)
   @test Array(b) ≈ Array(a) / 3
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   @test @constinferred(norm(a)) ≈ norm(Array(a))
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   if arrayt === Array
     b = @constinferred exp(a)
     @test Array(b) ≈ exp(Array(a))
@@ -267,21 +305,26 @@ end
     @test_broken exp(a)
   end
 
-  ## if VERSION < v"1.11-" && elt <: Complex
-  ##   # Broken because of type stability issue in Julia v1.10.
-  ##   @test_broken svd_compact(a)
-  if arrayt === Array
-    u, s, v = svd_compact(a)
-    @test u * s * v ≈ a
-    @test blocktype(u) >: blocktype(u)
-    @test eltype(u) === eltype(a)
-    @test blocktype(v) >: blocktype(a)
-    @test eltype(v) === eltype(a)
-    @test eltype(s) === real(eltype(a))
-  else
-    @test_broken svd_compact(a)
-  end
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
+  u, s, v = svd_compact(a)
+  @test u * s * v ≈ a
+  @test blocktype(u) >: blocktype(u)
+  @test eltype(u) === eltype(a)
+  @test blocktype(v) >: blocktype(a)
+  @test eltype(v) === eltype(a)
+  @test eltype(s) === real(eltype(a))
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   if arrayt === Array
     @test Array(inv(a)) ≈ inv(Array(a))
   else
@@ -289,6 +332,12 @@ end
     @test_broken inv(a)
   end
 
+  r = blockrange([2 × 2, 3 × 3])
+  d = Dict(
+    Block(1, 1) => dev(Eye{elt}(2, 2) ⊗ randn(elt, 2, 2)),
+    Block(2, 2) => dev(Eye{elt}(3, 3) ⊗ randn(elt, 3, 3)),
+  )
+  a = dev(blocksparse(d, (r, r)))
   # Broken operations
   b = a[Block.(1:2), Block(2)]
   @test b[Block(1)] == a[Block(1, 2)]