From 478151007de1cb54b0079157b1b2ce565ee2713a Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 16:37:32 -0400
Subject: [PATCH 01/18] rename file to `utility`

---
 src/BlockSparseArrays.jl                             | 2 +-
 src/factorizations/{tensorproducts.jl => utility.jl} | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename src/factorizations/{tensorproducts.jl => utility.jl} (100%)

diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index 4229343..8949bad 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -45,7 +45,7 @@ include("blocksparsearray/blockdiagonalarray.jl")
 include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl")
 
 # factorizations
-include("factorizations/tensorproducts.jl")
+include("factorizations/utility.jl")
 include("factorizations/svd.jl")
 include("factorizations/truncation.jl")
 include("factorizations/qr.jl")
diff --git a/src/factorizations/tensorproducts.jl b/src/factorizations/utility.jl
similarity index 100%
rename from src/factorizations/tensorproducts.jl
rename to src/factorizations/utility.jl

From 8c4278618348390819feaffa8353c28ce6e52497 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 16:37:42 -0400
Subject: [PATCH 02/18] Add blockdiagonalization utility functions

---
 src/factorizations/utility.jl | 33 +++++++++++++++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index bbf02ae..f703f46 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -17,3 +17,36 @@ function supremum(r1::AbstractUnitRange, r2::AbstractUnitRange)
     return r2
   end
 end
+
+function blockdiagonalize(A::AbstractBlockSparseMatrix)
+  # sort in order to avoid depending on internal details such as dictionary order
+  bIs = sort!(collect(eachblockstoredindex(A)); by=Int ∘ last ∘ Tuple)
+
+  # obtain permutation for the blocks that are present
+  rowperm = map(first ∘ Tuple, bIs)
+  colperm = map(last ∘ Tuple, bIs)
+
+  # These checks might be expensive but are convenient for now
+  allunique(rowperm) || throw(ArgumentError("input contains more than one block per row"))
+  allunique(colperm) ||
+    throw(ArgumentError("input contains more than one block per column"))
+
+  # post-process empty rows and columns: this pairs them up in order of occurance,
+  # putting empty rows and columns due to rectangularity at the end
+  emptyrows = setdiff(Block.(1:blocksize(A, 1)), rowperm)
+  append!(rowperm, emptyrows)
+  emptycols = setdiff(Block.(1:blocksize(A, 2)), colperm)
+  append!(colperm, emptycols)
+
+  return A[rowperm, colperm], rowperm, colperm
+end
+
+function isblockdiagonal(A::AbstractBlockSparseMatrix)
+  for bI in eachblockstoredindex(A)
+    row, col = Tuple(bI)
+    row == col || return false
+  end
+  # don't need to check for rows and cols appearing only once
+  # this is guaranteed as long as eachblockstoredindex is unique, which we assume
+  return true
+end

From a3af231beadac5f70928b5b1f13022d88bcdde82 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 16:43:52 -0400
Subject: [PATCH 03/18] move BlockPermutedDiagonalAlgorithm and update
 docstring

---
 src/factorizations/svd.jl     | 17 -----------------
 src/factorizations/utility.jl | 22 ++++++++++++++++++++++
 2 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 8c616d4..9b2c25d 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -3,23 +3,6 @@ using MatrixAlgebraKit:
   MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!
 using TypeParameterAccessors: realtype
 
-"""
-    BlockPermutedDiagonalAlgorithm(A::MatrixAlgebraKit.AbstractAlgorithm)
-  
-A wrapper for `MatrixAlgebraKit.AbstractAlgorithm` that implements the wrapped algorithm on
-a block-by-block basis, which is possible if the input matrix is a block-diagonal matrix or
-a block permuted block-diagonal matrix.
-"""
-struct BlockPermutedDiagonalAlgorithm{F} <: MatrixAlgebraKit.AbstractAlgorithm
-  falg::F
-end
-function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, a::AbstractMatrix)
-  return block_algorithm(alg, typeof(a))
-end
-function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, A::Type{<:AbstractMatrix})
-  return alg.falg(A)
-end
-
 function MatrixAlgebraKit.default_svd_algorithm(
   ::Type{<:AbstractBlockSparseMatrix}; kwargs...
 )
diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index f703f46..75ac825 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -1,3 +1,5 @@
+using MatrixAlgebraKit: MatrixAlgebraKit
+
 function infimum(r1::AbstractUnitRange, r2::AbstractUnitRange)
   (isone(first(r1)) && isone(first(r2))) ||
     throw(ArgumentError("infimum only defined for ranges starting at 1"))
@@ -50,3 +52,23 @@ function isblockdiagonal(A::AbstractBlockSparseMatrix)
   # this is guaranteed as long as eachblockstoredindex is unique, which we assume
   return true
 end
+
+"""
+    BlockPermutedDiagonalAlgorithm([f]) <: MatrixAlgebraKit.AbstractAlgorithm
+
+Type for handling algorithms on a block-by-block basis, which is possible for
+block-diagonal or block-permuted-diagonal input matrices.
+Additionally this wrapper may take a function that, given the input type of the blocks,
+returns the actual algorithm that will be used. This can be leveraged to allow for different
+algorithms for each block.
+"""
+struct BlockPermutedDiagonalAlgorithm{F} <: MatrixAlgebraKit.AbstractAlgorithm
+  falg::F
+end
+
+function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, a::AbstractMatrix)
+  return block_algorithm(alg, typeof(a))
+end
+function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, A::Type{<:AbstractMatrix})
+  return alg.falg(A)
+end

From bae92dab499a19f03764b180e1f0de771abb4681 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 16:46:33 -0400
Subject: [PATCH 04/18] Split blockdiagonal and blockpermuteddiagonal
 algorithms

---
 src/factorizations/utility.jl | 34 +++++++++++++++++++++++++++++-----
 1 file changed, 29 insertions(+), 5 deletions(-)

diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index 75ac825..8bcbd7b 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -53,22 +53,46 @@ function isblockdiagonal(A::AbstractBlockSparseMatrix)
   return true
 end
 
+"""
+    BlockDiagonalAlgorithm([f]) <: MatrixAlgebraKit.AbstractAlgorithm
+
+Type for handling algorithms on a block-by-block basis, which is possible for
+block-diagonal input matrices.
+Additionally this algorithm may take a function that, given the individual blocks, returns
+the algorithm that will be used. This can be leveraged to allow for different algorithms for
+each block.
+"""
+struct BlockDiagonalAlgorithm{F} <: MatrixAlgebraKit.AbstractAlgorithm
+  falg::F
+end
+
+function block_algorithm(alg::BlockDiagonalAlgorithm, a::AbstractMatrix)
+  return block_algorithm(alg, typeof(a))
+end
+function block_algorithm(alg::BlockDiagonalAlgorithm, A::Type{<:AbstractMatrix})
+  return alg.falg(A)
+end
+
 """
     BlockPermutedDiagonalAlgorithm([f]) <: MatrixAlgebraKit.AbstractAlgorithm
 
 Type for handling algorithms on a block-by-block basis, which is possible for
-block-diagonal or block-permuted-diagonal input matrices.
-Additionally this wrapper may take a function that, given the input type of the blocks,
-returns the actual algorithm that will be used. This can be leveraged to allow for different
-algorithms for each block.
+block-diagonal or block-permuted-diagonal input matrices. The algorithms proceed by first
+permuting to a block-diagonal form, and then carrying out the algorithm.
+Additionally this algorithm may take a function that, given the individual blocks, returns
+the algorithm that will be used. This can be leveraged to allow for different algorithms for
+each block.
 """
 struct BlockPermutedDiagonalAlgorithm{F} <: MatrixAlgebraKit.AbstractAlgorithm
   falg::F
 end
-
 function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, a::AbstractMatrix)
   return block_algorithm(alg, typeof(a))
 end
 function block_algorithm(alg::BlockPermutedDiagonalAlgorithm, A::Type{<:AbstractMatrix})
   return alg.falg(A)
 end
+
+function BlockDiagonalAlgorithm(alg::BlockPermutedDiagonalAlgorithm)
+  return BlockDiagonalAlgorithm(alg.falg)
+end

From 91bb6ae68e525f5297bfc47803201f79633ab03f Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 18:10:57 -0400
Subject: [PATCH 05/18] fixup! Add blockdiagonalization utility functions

---
 src/factorizations/utility.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index 8bcbd7b..ee19175 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -53,6 +53,11 @@ function isblockdiagonal(A::AbstractBlockSparseMatrix)
   return true
 end
 
+function isblockpermuteddiagonal(A::AbstractBlockSparseMatrix)
+  return allunique(first.(eachblockstoredindex(A))) &&
+         allunique(last.(eachblockstoredindex(A)))
+end
+
 """
     BlockDiagonalAlgorithm([f]) <: MatrixAlgebraKit.AbstractAlgorithm
 

From 08aa76f092848a91ab54b3c21ebd1b618b8747d1 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 18:11:58 -0400
Subject: [PATCH 06/18] refactor svd to work with new algorithm types

---
 src/factorizations/svd.jl | 251 +++++++++++++++++---------------------
 1 file changed, 111 insertions(+), 140 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 9b2c25d..723fd15 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -29,59 +29,34 @@ function similar_output(
 end
 
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(svd_compact!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
+)
+  return nothing
+end
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-  bmn = min(bm, bn)
-
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
-  u_axes = similar(brows, bmn)
-  v_axes = similar(brows, bmn)
+  # using the property that zip stops as soon as one of the iterators is exhausted
+  s_axes = map(splat(infimum), zip(brows, bcols))
+  s_axis = mortar_axis(s_axes)
+  S_axes = (s_axis, s_axis)
+  U, S, Vᴴ = similar_output(svd_compact!, A, S_axes, alg)
 
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
   for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    u_axes[col] = infimum(brows[row], bcols[col])
-    v_axes[col] = infimum(bcols[col], brows[row])
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    u_axes[col] = infimum(brows[row], bcols[col])
-    v_axes[col] = infimum(bcols[col], brows[row])
-  end
-
-  u_axis = mortar_axis(u_axes)
-  v_axis = mortar_axis(v_axes)
-  S_axes = (u_axis, v_axis)
-  U, S, Vt = similar_output(svd_compact!, A, S_axes, alg)
-
-  # allocate output
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
     block = @view!(A[bI])
     block_alg = block_algorithm(alg, block)
-    U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output(
+    I = first(Tuple(bI)) # == last(Tuple(bI))
+    U[I, I], S[I, I], Vᴴ[I, I] = MatrixAlgebraKit.initialize_output(
       svd_compact!, block, block_alg
     )
   end
 
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(U[Block(row, col)])
-    @view!(Vt[Block(col, col)])
-  end
-
-  return U, S, Vt
+  return U, S, Vᴴ
 end
 
+
 function similar_output(
   ::typeof(svd_full!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm
 )
@@ -93,65 +68,39 @@ function similar_output(
 end
 
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(svd_full!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(svd_full!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-
-  brows = eachblockaxis(axes(A, 1))
-  u_axes = similar(brows)
-
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    u_axes[col] = brows[row]
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    u_axes[col] = brows[row]
-  end
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    u_axes[bn + i] = brows[emptyrows[k]]
-  end
+  return nothing
+end
 
-  u_axis = mortar_axis(u_axes)
-  S_axes = (u_axis, axes(A, 2))
-  U, S, Vt = similar_output(svd_full!, A, S_axes, alg)
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(svd_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
+)
+  U, S, Vᴴ = similar_output(svd_full!, A, axes(A), alg)
 
-  # allocate output
   for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
     block = @view!(A[bI])
     block_alg = block_algorithm(alg, block)
-    U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output(
+    I = first(Tuple(bI)) # == last(Tuple(bI))
+    U[I, I], S[I, I], Vᴴ[I, I] = MatrixAlgebraKit.initialize_output(
       svd_full!, block, block_alg
     )
   end
 
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(U[Block(row, col)])
-    @view!(Vt[Block(col, col)])
-  end
-  # also handle extra rows/cols
-  for i in (length(emptyrows) + 1):length(emptycols)
-    @view!(Vt[Block(emptycols[i], emptycols[i])])
-  end
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    @view!(U[Block(emptyrows[k], bn + i)])
-  end
+  return U, S, Vᴴ
+end
 
-  return U, S, Vt
+function MatrixAlgebraKit.check_input(
+  ::typeof(svd_compact!),
+  A::AbstractBlockSparseMatrix,
+  USVᴴ,
+  ::BlockPermutedDiagonalAlgorithm,
+)
+  @assert isblockpermuteddiagonal(A)
 end
 
 function MatrixAlgebraKit.check_input(
-  ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ)
+  ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ), ::BlockDiagonalAlgorithm
 )
   @assert isa(U, AbstractBlockSparseMatrix) &&
     isa(S, AbstractBlockSparseMatrix) &&
@@ -160,11 +109,19 @@ function MatrixAlgebraKit.check_input(
   @assert real(eltype(A)) == eltype(S)
   @assert axes(A, 1) == axes(U, 1) && axes(A, 2) == axes(Vᴴ, 2)
   @assert axes(S, 1) == axes(S, 2)
+  @assert isblockdiagonal(A)
   return nothing
 end
 
 function MatrixAlgebraKit.check_input(
-  ::typeof(svd_full!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ)
+  ::typeof(svd_full!), A::AbstractBlockSparseMatrix, USVᴴ, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+
+function MatrixAlgebraKit.check_input(
+  ::typeof(svd_full!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ), ::BlockDiagonalAlgorithm
 )
   @assert isa(U, AbstractBlockSparseMatrix) &&
     isa(S, AbstractBlockSparseMatrix) &&
@@ -173,78 +130,92 @@ function MatrixAlgebraKit.check_input(
   @assert real(eltype(A)) == eltype(S)
   @assert axes(A, 1) == axes(U, 1) && axes(A, 2) == axes(Vᴴ, 1) == axes(Vᴴ, 2)
   @assert axes(S, 2) == axes(A, 2)
+  @assert isblockdiagonal(A)
   return nothing
 end
 
 function MatrixAlgebraKit.svd_compact!(
-  A::AbstractBlockSparseMatrix, (U, S, Vᴴ), alg::BlockPermutedDiagonalAlgorithm
+  A::AbstractBlockSparseMatrix, USVᴴ, alg::BlockPermutedDiagonalAlgorithm
 )
-  check_input(svd_compact!, A, (U, S, Vᴴ))
+  check_input(svd_compact!, A, USVᴴ, alg)
 
-  # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    usvᴴ = (@view!(U[brow, bcol]), @view!(S[bcol, bcol]), @view!(Vᴴ[bcol, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    usvᴴ′ = svd_compact!(block, usvᴴ, block_alg)
-    @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place"
-  end
+  Ad, rowperm, colperm = blockdiagonalize(A)
+  Ud, S, Vᴴd = svd_compact!(Ad, BlockDiagonalAlgorithm(alg))
+
+  inv_rowperm = Block.(invperm(Int.(rowperm)))
+  U = Ud[inv_rowperm, :]
+
+  inv_colperm = Block.(invperm(Int.(colperm)))
+  Vᴴ = Vᴴd[:, inv_colperm]
+
+  return U, S, Vᴴ
+end
 
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # U[Block(row, col)] = LinearAlgebra.I
-  # Vᴴ[Block(col, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I)
-    copyto!(@view!(Vᴴ[Block(col, col)]), LinearAlgebra.I)
+function MatrixAlgebraKit.svd_compact!(
+  A::AbstractBlockSparseMatrix, (U, S, Vᴴ), alg::BlockDiagonalAlgorithm
+)
+  check_input(svd_compact!, A, (U, S, Vᴴ), alg)
+
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      usvᴴ = (@view!(U[bI]), @view!(S[bI]), @view!(Vᴴ[bI]))
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      usvᴴ′ = svd_compact!(block, usvᴴ, block_alg)
+      @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place"
+    else
+      copyto!(@view!(U[bI]), LinearAlgebra.I)
+      copyto!(@view!(Vᴴ[bI]), LinearAlgebra.I)
+    end
   end
 
-  return (U, S, Vᴴ)
+  return U, S, Vᴴ
 end
 
 function MatrixAlgebraKit.svd_full!(
-  A::AbstractBlockSparseMatrix, (U, S, Vᴴ), alg::BlockPermutedDiagonalAlgorithm
+  A::AbstractBlockSparseMatrix, USVᴴ, alg::BlockPermutedDiagonalAlgorithm
 )
-  check_input(svd_full!, A, (U, S, Vᴴ))
+  check_input(svd_full!, A, USVᴴ, alg)
 
-  # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    usvᴴ = (@view!(U[brow, bcol]), @view!(S[bcol, bcol]), @view!(Vᴴ[bcol, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    usvᴴ′ = svd_full!(block, usvᴴ, block_alg)
-    @assert usvᴴ === usvᴴ′ "svd_full! might not be in-place"
-  end
+  Ad, rowperm, colperm = blockdiagonalize(A)
+  Ud, S, Vᴴd = svd_full!(Ad, BlockDiagonalAlgorithm(alg))
+
+  inv_rowperm = Block.(invperm(Int.(rowperm)))
+  U = Ud[inv_rowperm, :]
 
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # U[Block(row, col)] = LinearAlgebra.I
-  # Vt[Block(col, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I)
-    copyto!(@view!(Vᴴ[Block(col, col)]), LinearAlgebra.I)
+  inv_colperm = Block.(invperm(Int.(colperm)))
+  Vᴴ = Vᴴd[:, inv_colperm]
+
+  return U, S, Vᴴ
+end
+
+function MatrixAlgebraKit.svd_full!(
+  A::AbstractBlockSparseMatrix, (U, S, Vᴴ), alg::BlockDiagonalAlgorithm
+)
+  check_input(svd_full!, A, (U, S, Vᴴ), alg)
+
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      usvᴴ = (@view!(U[bI]), @view!(S[bI]), @view!(Vᴴ[bI]))
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      usvᴴ′ = svd_full!(block, usvᴴ, block_alg)
+      @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place"
+    else
+      copyto!(@view!(U[bI]), LinearAlgebra.I)
+      copyto!(@view!(Vᴴ[bI]), LinearAlgebra.I)
+    end
   end
 
-  # also handle extra rows/cols
-  for i in (length(emptyrows) + 1):length(emptycols)
-    copyto!(@view!(Vᴴ[Block(emptycols[i], emptycols[i])]), LinearAlgebra.I)
+  # Complete the unitaries for rectangular inputs
+  for I in blocksize(A, 2)+1:blocksize(A, 1)
+    copyto!(@view!(U[Block(I, I)]), LinearAlgebra.I)
   end
-  bn = blocksize(A, 2)
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    copyto!(@view!(U[Block(emptyrows[k], bn + i)]), LinearAlgebra.I)
+  for I in blocksize(A, 1)+1:blocksize(A, 2)
+    copyto!(@view!(Vᴴ[Block(I, I)]), LinearAlgebra.I)
   end
 
-  return (U, S, Vᴴ)
+  return U, S, Vᴴ
 end

From d7911bcddb934050edc998bc48917b795c6b4959 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 18:12:10 -0400
Subject: [PATCH 07/18] slight refactor of tests

---
 test/test_factorizations.jl | 30 +++++++++++++++---------------
 1 file changed, 15 insertions(+), 15 deletions(-)

diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index 651a8f6..cb60b26 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -119,15 +119,15 @@ end
 
 function test_svd(a, (U, S, Vᴴ); full=false)
   # Check that the SVD is correct
-  (U * S * Vᴴ ≈ a) || return false
-  (U' * U ≈ LinearAlgebra.I) || return false
-  (Vᴴ * Vᴴ' ≈ LinearAlgebra.I) || return false
-  full || return true
+  @test (U * S * Vᴴ ≈ a)
+  @test (U' * U ≈ LinearAlgebra.I)
+  @test (Vᴴ * Vᴴ' ≈ LinearAlgebra.I)
+  full || return nothing
 
   # Check factors are unitary
-  (U * U' ≈ LinearAlgebra.I) || return false
-  (Vᴴ' * Vᴴ ≈ LinearAlgebra.I) || return false
-  return true
+  @test (U * U' ≈ LinearAlgebra.I)
+  @test (Vᴴ' * Vᴴ ≈ LinearAlgebra.I)
+  return nothing
 end
 
 blockszs = (
@@ -143,7 +143,7 @@ test_params = Iterators.product(blockszs, eltypes)
 
   # test empty matrix
   usv_empty = svd_compact(a)
-  @test test_svd(a, usv_empty)
+  test_svd(a, usv_empty)
 
   # test blockdiagonal
   rng = StableRNG(123)
@@ -152,13 +152,13 @@ test_params = Iterators.product(blockszs, eltypes)
     a[Block(I.I...)] = rand(rng, T, size(blocks(a)[i]))
   end
   usv = svd_compact(a)
-  @test test_svd(a, usv)
+  test_svd(a, usv)
 
   rng = StableRNG(123)
   perm = Random.randperm(rng, length(m))
   b = a[Block.(perm), Block.(1:length(n))]
   usv = svd_compact(b)
-  @test test_svd(b, usv)
+  test_svd(b, usv)
 
   # test permuted blockdiagonal with missing row/col
   rng = StableRNG(123)
@@ -166,7 +166,7 @@ test_params = Iterators.product(blockszs, eltypes)
   c = copy(b)
   delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed))))
   usv = svd_compact(c)
-  @test test_svd(c, usv)
+  test_svd(c, usv)
 end
 
 # svd_full!
@@ -176,7 +176,7 @@ end
 
   # test empty matrix
   usv_empty = svd_full(a)
-  @test test_svd(a, usv_empty; full=true)
+  test_svd(a, usv_empty; full=true)
 
   # test blockdiagonal
   rng = StableRNG(123)
@@ -185,13 +185,13 @@ end
     a[Block(I.I...)] = rand(rng, T, size(blocks(a)[i]))
   end
   usv = svd_full(a)
-  @test test_svd(a, usv; full=true)
+  test_svd(a, usv; full=true)
 
   rng = StableRNG(123)
   perm = Random.randperm(rng, length(m))
   b = a[Block.(perm), Block.(1:length(n))]
   usv = svd_full(b)
-  @test test_svd(b, usv; full=true)
+  test_svd(b, usv; full=true)
 
   # test permuted blockdiagonal with missing row/col
   rng = StableRNG(123)
@@ -199,7 +199,7 @@ end
   c = copy(b)
   delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed))))
   usv = svd_full(c)
-  @test test_svd(c, usv; full=true)
+  test_svd(c, usv; full=true)
 end
 
 # svd_trunc!

From 535d78cf88d8d8a1f16e666fd90e633783156145 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 4 Aug 2025 18:12:14 -0400
Subject: [PATCH 08/18] add test of issues

---
 test/test_issues.jl | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)
 create mode 100644 test/test_issues.jl

diff --git a/test/test_issues.jl b/test/test_issues.jl
new file mode 100644
index 0000000..aa356ea
--- /dev/null
+++ b/test/test_issues.jl
@@ -0,0 +1,26 @@
+using BlockArrays
+using BlockSparseArrays
+using BlockSparseArrays: blocksparse
+using MatrixAlgebraKit
+using LinearAlgebra: LinearAlgebra
+using Test: @test, @testset
+
+@testset "Issue 162" begin
+  ax = (blockedrange([2, 2]), blockedrange([2, 2, 2]))
+  bs = Dict(Block(1, 1) => randn(2, 2), Block(2, 3) => randn(2, 2))
+  a = blocksparse(bs, ax)
+  U, S, Vᴴ = svd_compact(a)
+
+  @test U * S * Vᴴ ≈ a
+  @test U' * U ≈ LinearAlgebra.I
+  @test Vᴴ * Vᴴ' ≈ LinearAlgebra.I
+
+  U, S, Vᴴ = svd_full(a);
+
+  @test U * S * Vᴴ ≈ a
+  @test U' * U ≈ LinearAlgebra.I
+  @test U * U' ≈ LinearAlgebra.I
+  @test Vᴴ * Vᴴ' ≈ LinearAlgebra.I
+  @test Vᴴ' * Vᴴ ≈ LinearAlgebra.I
+end
+

From 54570950f9c0a9184c52d81f44ab9937d3fb1f20 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 4 Aug 2025 22:19:39 +0000
Subject: [PATCH 09/18] Formatter

---
 src/factorizations/svd.jl | 5 ++---
 test/test_issues.jl       | 1 -
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 723fd15..578c518 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -56,7 +56,6 @@ function MatrixAlgebraKit.initialize_output(
   return U, S, Vᴴ
 end
 
-
 function similar_output(
   ::typeof(svd_full!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm
 )
@@ -210,10 +209,10 @@ function MatrixAlgebraKit.svd_full!(
   end
 
   # Complete the unitaries for rectangular inputs
-  for I in blocksize(A, 2)+1:blocksize(A, 1)
+  for I in (blocksize(A, 2) + 1):blocksize(A, 1)
     copyto!(@view!(U[Block(I, I)]), LinearAlgebra.I)
   end
-  for I in blocksize(A, 1)+1:blocksize(A, 2)
+  for I in (blocksize(A, 1) + 1):blocksize(A, 2)
     copyto!(@view!(Vᴴ[Block(I, I)]), LinearAlgebra.I)
   end
 
diff --git a/test/test_issues.jl b/test/test_issues.jl
index aa356ea..0ce92ca 100644
--- a/test/test_issues.jl
+++ b/test/test_issues.jl
@@ -23,4 +23,3 @@ using Test: @test, @testset
   @test Vᴴ * Vᴴ' ≈ LinearAlgebra.I
   @test Vᴴ' * Vᴴ ≈ LinearAlgebra.I
 end
-

From 4c632e83136c317445250e24446c398c55ac5d82 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 12 Aug 2025 11:36:02 +0200
Subject: [PATCH 10/18] change transform into functions

---
 src/factorizations/svd.jl     | 20 ++++++--------------
 src/factorizations/utility.jl |  8 +++++++-
 2 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 578c518..fd93a5c 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -138,14 +138,10 @@ function MatrixAlgebraKit.svd_compact!(
 )
   check_input(svd_compact!, A, USVᴴ, alg)
 
-  Ad, rowperm, colperm = blockdiagonalize(A)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
   Ud, S, Vᴴd = svd_compact!(Ad, BlockDiagonalAlgorithm(alg))
-
-  inv_rowperm = Block.(invperm(Int.(rowperm)))
-  U = Ud[inv_rowperm, :]
-
-  inv_colperm = Block.(invperm(Int.(colperm)))
-  Vᴴ = Vᴴd[:, inv_colperm]
+  U = transform_rows(Ud)
+  Vᴴ = transform_cols(Vᴴd)
 
   return U, S, Vᴴ
 end
@@ -177,14 +173,10 @@ function MatrixAlgebraKit.svd_full!(
 )
   check_input(svd_full!, A, USVᴴ, alg)
 
-  Ad, rowperm, colperm = blockdiagonalize(A)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
   Ud, S, Vᴴd = svd_full!(Ad, BlockDiagonalAlgorithm(alg))
-
-  inv_rowperm = Block.(invperm(Int.(rowperm)))
-  U = Ud[inv_rowperm, :]
-
-  inv_colperm = Block.(invperm(Int.(colperm)))
-  Vᴴ = Vᴴd[:, inv_colperm]
+  U = transform_rows(Ud)
+  Vᴴ = transform_cols(Vᴴd)
 
   return U, S, Vᴴ
 end
diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index ee19175..f56b44b 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -40,7 +40,13 @@ function blockdiagonalize(A::AbstractBlockSparseMatrix)
   emptycols = setdiff(Block.(1:blocksize(A, 2)), colperm)
   append!(colperm, emptycols)
 
-  return A[rowperm, colperm], rowperm, colperm
+  invrowperm = Block.(invperm(Int.(rowperm)))
+  transform_rows(A) = A[invrowperm, :]
+
+  invcolperm = Block.(invperm(Int.(colperm)))
+  transform_cols(A) = A[:, invcolperm]
+
+  return A[rowperm, colperm], transform_rows, transform_cols
 end
 
 function isblockdiagonal(A::AbstractBlockSparseMatrix)

From 2fa94761619de9a44fb47ac9eb998eb008df802f Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 13 Aug 2025 11:16:26 +0200
Subject: [PATCH 11/18] simplify initializations

---
 src/factorizations/svd.jl | 73 ++++++++++++---------------------------
 1 file changed, 22 insertions(+), 51 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index fd93a5c..6c5f6b2 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -11,21 +11,11 @@ function MatrixAlgebraKit.default_svd_algorithm(
   end
 end
 
-function output_type(::typeof(svd_compact!), A::Type{<:AbstractMatrix{T}}) where {T}
-  USVᴴ = Base.promote_op(svd_compact!, A)
-  !isconcretetype(USVᴴ) &&
-    return Tuple{AbstractMatrix{T},AbstractMatrix{realtype(T)},AbstractMatrix{T}}
-  return USVᴴ
-end
-
-function similar_output(
-  ::typeof(svd_compact!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm
-)
-  BU, BS, BVᴴ = fieldtypes(output_type(svd_compact!, blocktype(A)))
-  U = similar(A, BlockType(BU), (axes(A, 1), S_axes[1]))
-  S = similar(A, BlockType(BS), S_axes)
-  Vᴴ = similar(A, BlockType(BVᴴ), (S_axes[2], axes(A, 2)))
-  return U, S, Vᴴ
+function output_type(
+  f::Union{typeof(svd_compact!),typeof(svd_full!)}, A::Type{<:AbstractMatrix{T}}
+) where {T}
+  USVᴴ = Base.promote_op(f, A)
+  return isconcretetype(USVᴴ) ? USVᴴ : Tuple{AbstractMatrix{T},AbstractMatrix{realtype(T)},AbstractMatrix{T}}
 end
 
 function MatrixAlgebraKit.initialize_output(
@@ -42,28 +32,13 @@ function MatrixAlgebraKit.initialize_output(
   s_axes = map(splat(infimum), zip(brows, bcols))
   s_axis = mortar_axis(s_axes)
   S_axes = (s_axis, s_axis)
-  U, S, Vᴴ = similar_output(svd_compact!, A, S_axes, alg)
-
-  for bI in eachblockstoredindex(A)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    I = first(Tuple(bI)) # == last(Tuple(bI))
-    U[I, I], S[I, I], Vᴴ[I, I] = MatrixAlgebraKit.initialize_output(
-      svd_compact!, block, block_alg
-    )
-  end
 
-  return U, S, Vᴴ
-end
+  BU, BS, BVᴴ = fieldtypes(output_type(svd_compact!, blocktype(A)))
+  U = similar(A, BlockType(BU), (axes(A, 1), S_axes[1]))
+  S = similar(A, BlockType(BS), S_axes)
+  Vᴴ = similar(A, BlockType(BVᴴ), (S_axes[2], axes(A, 2)))
 
-function similar_output(
-  ::typeof(svd_full!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm
-)
-  U = similar(A, axes(A, 1), S_axes[1])
-  T = real(eltype(A))
-  S = similar(A, T, S_axes)
-  Vt = similar(A, S_axes[2], axes(A, 2))
-  return U, S, Vt
+  return U, S, Vᴴ
 end
 
 function MatrixAlgebraKit.initialize_output(
@@ -75,16 +50,10 @@ end
 function MatrixAlgebraKit.initialize_output(
   ::typeof(svd_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
 )
-  U, S, Vᴴ = similar_output(svd_full!, A, axes(A), alg)
-
-  for bI in eachblockstoredindex(A)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    I = first(Tuple(bI)) # == last(Tuple(bI))
-    U[I, I], S[I, I], Vᴴ[I, I] = MatrixAlgebraKit.initialize_output(
-      svd_full!, block, block_alg
-    )
-  end
+  BU, BS, BVᴴ = fieldtypes(output_type(svd_full!, blocktype(A)))
+  U = similar(A, BlockType(BU), (axes(A, 1), axes(A, 1)))
+  S = similar(A, BlockType(BS), axes(A))
+  Vᴴ = similar(A, BlockType(BVᴴ), (axes(A, 2), axes(A, 2)))
 
   return U, S, Vᴴ
 end
@@ -154,11 +123,12 @@ function MatrixAlgebraKit.svd_compact!(
   for I in 1:min(blocksize(A)...)
     bI = Block(I, I)
     if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
-      usvᴴ = (@view!(U[bI]), @view!(S[bI]), @view!(Vᴴ[bI]))
       block = @view!(A[bI])
       block_alg = block_algorithm(alg, block)
-      usvᴴ′ = svd_compact!(block, usvᴴ, block_alg)
-      @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place"
+      bU, bS, bVᴴ = svd_compact!(block, block_alg)
+      U[bI] = bU
+      S[bI] = bS
+      Vᴴ[bI] = bVᴴ
     else
       copyto!(@view!(U[bI]), LinearAlgebra.I)
       copyto!(@view!(Vᴴ[bI]), LinearAlgebra.I)
@@ -189,11 +159,12 @@ function MatrixAlgebraKit.svd_full!(
   for I in 1:min(blocksize(A)...)
     bI = Block(I, I)
     if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
-      usvᴴ = (@view!(U[bI]), @view!(S[bI]), @view!(Vᴴ[bI]))
       block = @view!(A[bI])
       block_alg = block_algorithm(alg, block)
-      usvᴴ′ = svd_full!(block, usvᴴ, block_alg)
-      @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place"
+      bU, bS, bVᴴ = svd_full!(block, block_alg)
+      U[bI] = bU
+      S[bI] = bS
+      Vᴴ[bI] = bVᴴ
     else
       copyto!(@view!(U[bI]), LinearAlgebra.I)
       copyto!(@view!(Vᴴ[bI]), LinearAlgebra.I)

From cecf9673954c44e6806684201bbfa42e6991e830 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 13 Aug 2025 11:54:54 +0200
Subject: [PATCH 12/18] Refactor QR

---
 src/factorizations/qr.jl | 247 ++++++++++++++-------------------------
 1 file changed, 89 insertions(+), 158 deletions(-)

diff --git a/src/factorizations/qr.jl b/src/factorizations/qr.jl
index 10c40d8..99804aa 100644
--- a/src/factorizations/qr.jl
+++ b/src/factorizations/qr.jl
@@ -1,5 +1,4 @@
-using MatrixAlgebraKit:
-  MatrixAlgebraKit, default_qr_algorithm, lq_compact!, lq_full!, qr_compact!, qr_full!
+using MatrixAlgebraKit: MatrixAlgebraKit, default_qr_algorithm, qr_compact!, qr_full!
 
 function MatrixAlgebraKit.default_qr_algorithm(
   ::Type{<:AbstractBlockSparseMatrix}; kwargs...
@@ -9,211 +8,143 @@ function MatrixAlgebraKit.default_qr_algorithm(
   end
 end
 
-function similar_output(
-  ::typeof(qr_compact!), A, R_axis, alg::MatrixAlgebraKit.AbstractAlgorithm
-)
-  Q = similar(A, axes(A, 1), R_axis)
-  R = similar(A, R_axis, axes(A, 2))
-  return Q, R
+function output_type(
+  f::Union{typeof(qr_compact!),typeof(qr_full!)}, A::Type{<:AbstractMatrix{T}}
+) where {T}
+  QR = Base.promote_op(f, A)
+  return isconcretetype(QR) ? QR : Tuple{AbstractMatrix{T},AbstractMatrix{T}}
 end
 
-function similar_output(
-  ::typeof(qr_full!), A, R_axis, alg::MatrixAlgebraKit.AbstractAlgorithm
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(qr_compact!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
-  Q = similar(A, axes(A, 1), R_axis)
-  R = similar(A, R_axis, axes(A, 2))
-  return Q, R
+  return nothing
 end
-
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(qr_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(qr_compact!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-  bmn = min(bm, bn)
-
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
-  r_axes = similar(brows, bmn)
-
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    len = minimum(length, (brows[row], bcols[col]))
-    r_axes[col] = brows[row][Base.OneTo(len)]
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    len = minimum(length, (brows[row], bcols[col]))
-    r_axes[col] = brows[row][Base.OneTo(len)]
-  end
-
+  # using the property that zip stops as soon as one of the iterators is exhausted
+  r_axes = map(splat(infimum), zip(brows, bcols))
   r_axis = mortar_axis(r_axes)
-  Q, R = similar_output(qr_compact!, A, r_axis, alg)
-
-  # allocate output
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    Q[brow, bcol], R[bcol, bcol] = MatrixAlgebraKit.initialize_output(
-      qr_compact!, block, block_alg
-    )
-  end
 
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(Q[Block(row, col)])
-  end
+  BQ, BR = fieldtypes(output_type(qr_compact!, blocktype(A)))
+  Q = similar(A, BlockType(BQ), (axes(A, 1), r_axis))
+  R = similar(A, BlockType(BR), (r_axis, axes(A, 2)))
 
   return Q, R
 end
 
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(qr_full!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(qr_full!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-
-  brows = eachblockaxis(axes(A, 1))
-  r_axes = copy(brows)
-
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    r_axes[col] = brows[row]
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    r_axes[col] = brows[row]
-  end
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    r_axes[bn + i] = brows[emptyrows[k]]
-  end
-
-  r_axis = mortar_axis(r_axes)
-  Q, R = similar_output(qr_full!, A, r_axis, alg)
-
-  # allocate output
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    Q[brow, bcol], R[bcol, bcol] = MatrixAlgebraKit.initialize_output(
-      qr_full!, block, block_alg
-    )
-  end
-
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(Q[Block(row, col)])
-  end
-  # also handle extra rows/cols
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    @view!(Q[Block(emptyrows[k], bn + i)])
-  end
-
+  return nothing
+end
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(qr_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
+)
+  BQ, BR = fieldtypes(output_type(qr_full!, blocktype(A)))
+  Q = similar(A, BlockType(BQ), (axes(A, 1), axes(A, 1)))
+  R = similar(A, BlockType(BR), (axes(A, 1), axes(A, 2)))
   return Q, R
 end
 
 function MatrixAlgebraKit.check_input(
-  ::typeof(qr_compact!), A::AbstractBlockSparseMatrix, QR
+  ::typeof(qr_compact!), A::AbstractBlockSparseMatrix, QR, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(qr_compact!), A::AbstractBlockSparseMatrix, (Q, R), ::BlockDiagonalAlgorithm
 )
-  Q, R = QR
   @assert isa(Q, AbstractBlockSparseMatrix) && isa(R, AbstractBlockSparseMatrix)
   @assert eltype(A) == eltype(Q) == eltype(R)
   @assert axes(A, 1) == axes(Q, 1) && axes(A, 2) == axes(R, 2)
   @assert axes(Q, 2) == axes(R, 1)
-
+  @assert isblockdiagonal(A)
   return nothing
 end
 
-function MatrixAlgebraKit.check_input(::typeof(qr_full!), A::AbstractBlockSparseMatrix, QR)
-  Q, R = QR
+function MatrixAlgebraKit.check_input(
+  ::typeof(qr_full!), A::AbstractBlockSparseMatrix, QR, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(qr_full!), A::AbstractBlockSparseMatrix, (Q, R), ::BlockDiagonalAlgorithm
+)
   @assert isa(Q, AbstractBlockSparseMatrix) && isa(R, AbstractBlockSparseMatrix)
   @assert eltype(A) == eltype(Q) == eltype(R)
   @assert axes(A, 1) == axes(Q, 1) && axes(A, 2) == axes(R, 2)
   @assert axes(Q, 2) == axes(R, 1)
-
+  @assert isblockdiagonal(A)
   return nothing
 end
 
 function MatrixAlgebraKit.qr_compact!(
   A::AbstractBlockSparseMatrix, QR, alg::BlockPermutedDiagonalAlgorithm
 )
-  MatrixAlgebraKit.check_input(qr_compact!, A, QR)
-  Q, R = QR
+  check_input(qr_compact!, A, QR, alg)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Qd, Rd = qr_compact!(Ad, BlockDiagonalAlgorithm(alg))
+  Q = transform_rows(Qd)
+  R = transform_cols(Rd)
+  return Q, R
+end
 
-  # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    qr = (@view!(Q[brow, bcol]), @view!(R[bcol, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    qr′ = qr_compact!(block, qr, block_alg)
-    @assert qr === qr′ "qr_compact! might not be in-place"
-  end
+function MatrixAlgebraKit.qr_compact!(
+  A::AbstractBlockSparseMatrix, (Q, R), alg::BlockDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(qr_compact!, A, (Q, R), alg)
 
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # Q[Block(row, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(Q[Block(row, col)]), LinearAlgebra.I)
+  # do decomposition on each block
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      bQ, bR = qr_compact!(block, block_alg)
+      Q[bI] = bQ
+      R[bI] = bR
+    else
+      copyto!(@view!(Q[bI]), LinearAlgebra.I)
+    end
   end
 
-  return QR
+  return Q, R
 end
 
 function MatrixAlgebraKit.qr_full!(
   A::AbstractBlockSparseMatrix, QR, alg::BlockPermutedDiagonalAlgorithm
 )
-  MatrixAlgebraKit.check_input(qr_full!, A, QR)
-  Q, R = QR
-
-  # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    qr = (@view!(Q[brow, bcol]), @view!(R[bcol, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    qr′ = qr_full!(block, qr, block_alg)
-    @assert qr === qr′ "qr_full! might not be in-place"
-  end
+  check_input(qr_full!, A, QR, alg)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Qd, Rd = qr_full!(Ad, BlockDiagonalAlgorithm(alg))
+  Q = transform_rows(Qd)
+  R = transform_cols(Rd)
+  return Q, R
+end
 
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # Q[Block(row, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(Q[Block(row, col)]), LinearAlgebra.I)
-  end
+function MatrixAlgebraKit.qr_full!(
+  A::AbstractBlockSparseMatrix, (Q, R), alg::BlockDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(qr_full!, A, (Q, R), alg)
 
-  # also handle extra rows/cols
-  bn = blocksize(A, 2)
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    copyto!(@view!(Q[Block(emptyrows[k], bn + i)]), LinearAlgebra.I)
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      bQ, bR = qr_full!(block, block_alg)
+      Q[bI] = bQ
+      R[bI] = bR
+    else
+      copyto!(@view!(Q[bI]), LinearAlgebra.I)
+    end
   end
 
-  return QR
+  return Q, R
 end

From d926f52a1ed8d5b276e2379147589ef15ef3d757 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 12:39:31 +0200
Subject: [PATCH 13/18] Refactor LQ

---
 src/factorizations/lq.jl | 245 ++++++++++++++-------------------------
 1 file changed, 87 insertions(+), 158 deletions(-)

diff --git a/src/factorizations/lq.jl b/src/factorizations/lq.jl
index 5dc1301..306849a 100644
--- a/src/factorizations/lq.jl
+++ b/src/factorizations/lq.jl
@@ -8,211 +8,140 @@ function MatrixAlgebraKit.default_lq_algorithm(
   end
 end
 
-function similar_output(
-  ::typeof(lq_compact!), A, L_axis, alg::MatrixAlgebraKit.AbstractAlgorithm
-)
-  L = similar(A, axes(A, 1), L_axis)
-  Q = similar(A, L_axis, axes(A, 2))
-  return L, Q
+function output_type(
+  f::Union{typeof(lq_compact!),typeof(lq_full!)}, A::Type{<:AbstractMatrix{T}}
+) where {T}
+  LQ = Base.promote_op(f, A)
+  return isconcretetype(LQ) ? LQ : Tuple{AbstractMatrix{T},AbstractMatrix{T}}
 end
 
-function similar_output(
-  ::typeof(lq_full!), A, L_axis, alg::MatrixAlgebraKit.AbstractAlgorithm
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(lq_compact!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
-  L = similar(A, axes(A, 1), L_axis)
-  Q = similar(A, L_axis, axes(A, 2))
-  return L, Q
+  return nothing
 end
-
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-  bmn = min(bm, bn)
-
   brows = eachblockaxis(axes(A, 1))
   bcols = eachblockaxis(axes(A, 2))
-  l_axes = similar(brows, bmn)
-
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    len = minimum(length, (brows[row], bcols[col]))
-    l_axes[row] = bcols[col][Base.OneTo(len)]
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    len = minimum(length, (brows[row], bcols[col]))
-    l_axes[row] = bcols[col][Base.OneTo(len)]
-  end
-
+  # using the property that zip stops as soon as one of the iterators is exhausted
+  l_axes = map(splat(infimum), zip(brows, bcols))
   l_axis = mortar_axis(l_axes)
-  L, Q = similar_output(lq_compact!, A, l_axis, alg)
 
-  # allocate output
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    L[brow, brow], Q[brow, bcol] = MatrixAlgebraKit.initialize_output(
-      lq_compact!, block, block_alg
-    )
-  end
-
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(Q[Block(row, col)])
-  end
+  BL, BQ = fieldtypes(output_type(lq_compact!, blocktype(A)))
+  L = similar(A, BlockType(BL), (axes(A, 1), l_axis))
+  Q = similar(A, BlockType(BQ), (l_axis, axes(A, 2)))
 
   return L, Q
 end
 
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(lq_full!), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+  ::typeof(lq_full!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
-  bm, bn = blocksize(A)
-
-  bcols = eachblockaxis(axes(A, 2))
-  l_axes = copy(bcols)
-
-  # fill in values for blocks that are present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  for bI in eachblockstoredindex(A)
-    row, col = Int.(Tuple(bI))
-    l_axes[row] = bcols[col]
-  end
-
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  # this is a convention, which at least gives the expected results for blockdiagonal
-  emptyrows = setdiff(1:bm, browIs)
-  emptycols = setdiff(1:bn, bcolIs)
-  for (row, col) in zip(emptyrows, emptycols)
-    l_axes[row] = bcols[col]
-  end
-  for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows))
-    l_axes[bn + i] = bcols[emptycols[k]]
-  end
-
-  l_axis = mortar_axis(l_axes)
-  L, Q = similar_output(lq_full!, A, l_axis, alg)
-
-  # allocate output
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    L[brow, brow], Q[brow, bcol] = MatrixAlgebraKit.initialize_output(
-      lq_full!, block, block_alg
-    )
-  end
-
-  # allocate output for blocks that aren't present -- do we also fill identities here?
-  for (row, col) in zip(emptyrows, emptycols)
-    @view!(Q[Block(row, col)])
-  end
-  # also handle extra rows/cols
-  for (i, k) in enumerate((length(emptyrows) + 1):length(emptycols))
-    @view!(Q[Block(bm + i, emptycols[k])])
-  end
-
+  return nothing
+end
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(lq_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
+)
+  BL, BQ = fieldtypes(output_type(lq_full!, blocktype(A)))
+  L = similar(A, BlockType(BL), (axes(A, 1), axes(A, 2)))
+  Q = similar(A, BlockType(BQ), (axes(A, 2), axes(A, 2)))
   return L, Q
 end
 
 function MatrixAlgebraKit.check_input(
-  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix, LQ
+  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix, LQ, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix, (L, Q), ::BlockDiagonalAlgorithm
 )
-  L, Q = LQ
   @assert isa(L, AbstractBlockSparseMatrix) && isa(Q, AbstractBlockSparseMatrix)
   @assert eltype(A) == eltype(L) == eltype(Q)
   @assert axes(A, 1) == axes(L, 1) && axes(A, 2) == axes(Q, 2)
   @assert axes(L, 2) == axes(Q, 1)
-
+  @assert isblockdiagonal(A)
   return nothing
 end
-
-function MatrixAlgebraKit.check_input(::typeof(lq_full!), A::AbstractBlockSparseMatrix, LQ)
-  L, Q = LQ
+function MatrixAlgebraKit.check_input(
+  ::typeof(lq_full!), A::AbstractBlockSparseMatrix, LQ, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(lq_full!), A::AbstractBlockSparseMatrix, (L, Q), ::BlockDiagonalAlgorithm
+)
   @assert isa(L, AbstractBlockSparseMatrix) && isa(Q, AbstractBlockSparseMatrix)
   @assert eltype(A) == eltype(L) == eltype(Q)
   @assert axes(A, 1) == axes(L, 1) && axes(A, 2) == axes(Q, 2)
   @assert axes(L, 2) == axes(Q, 1)
-
+  @assert isblockdiagonal(A)
   return nothing
 end
 
 function MatrixAlgebraKit.lq_compact!(
   A::AbstractBlockSparseMatrix, LQ, alg::BlockPermutedDiagonalAlgorithm
 )
-  MatrixAlgebraKit.check_input(lq_compact!, A, LQ)
-  L, Q = LQ
+  MatrixAlgebraKit.check_input(lq_compact!, A, LQ, alg)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ld, Qd = lq_compact!(Ad, BlockDiagonalAlgorithm(alg))
+  L = transform_rows(Ld)
+  Q = transform_cols(Qd)
+  return L, Q
+end
+function MatrixAlgebraKit.lq_compact!(
+  A::AbstractBlockSparseMatrix, (L, Q), alg::BlockDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(lq_compact!, A, (L, Q), alg)
 
   # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    lq = (@view!(L[brow, brow]), @view!(Q[brow, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    lq′ = lq_compact!(block, lq, block_alg)
-    @assert lq === lq′ "lq_compact! might not be in-place"
-  end
-
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # Q[Block(row, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(Q[Block(row, col)]), LinearAlgebra.I)
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      bL, bQ = lq_compact!(block, block_alg)
+      L[bI] = bL
+      Q[bI] = bQ
+    else
+      copyto!(@view!(Q[bI]), LinearAlgebra.I)
+    end
   end
 
-  return LQ
+  return L, Q
 end
 
 function MatrixAlgebraKit.lq_full!(
   A::AbstractBlockSparseMatrix, LQ, alg::BlockPermutedDiagonalAlgorithm
 )
-  MatrixAlgebraKit.check_input(lq_full!, A, LQ)
-  L, Q = LQ
-
-  # do decomposition on each block
-  for bI in eachblockstoredindex(A)
-    brow, bcol = Tuple(bI)
-    lq = (@view!(L[brow, brow]), @view!(Q[brow, bcol]))
-    block = @view!(A[bI])
-    block_alg = block_algorithm(alg, block)
-    lq′ = lq_full!(block, lq, block_alg)
-    @assert lq === lq′ "lq_full! might not be in-place"
-  end
-
-  # fill in identities for blocks that aren't present
-  bIs = collect(eachblockstoredindex(A))
-  browIs = Int.(first.(Tuple.(bIs)))
-  bcolIs = Int.(last.(Tuple.(bIs)))
-  emptyrows = setdiff(1:blocksize(A, 1), browIs)
-  emptycols = setdiff(1:blocksize(A, 2), bcolIs)
-  # needs copyto! instead because size(::LinearAlgebra.I) doesn't work
-  # Q[Block(row, col)] = LinearAlgebra.I
-  for (row, col) in zip(emptyrows, emptycols)
-    copyto!(@view!(Q[Block(row, col)]), LinearAlgebra.I)
-  end
+  MatrixAlgebraKit.check_input(lq_full!, A, LQ, alg)
+  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ld, Qd = lq_full!(Ad, BlockDiagonalAlgorithm(alg))
+  L = transform_rows(Ld)
+  Q = transform_cols(Qd)
+  return L, Q
+end
+function MatrixAlgebraKit.lq_full!(
+  A::AbstractBlockSparseMatrix, (L, Q), alg::BlockDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(lq_full!, A, (L, Q), alg)
 
-  # also handle extra rows/cols
-  bm = blocksize(A, 1)
-  for (i, k) in enumerate((length(emptyrows) + 1):length(emptycols))
-    copyto!(@view!(Q[Block(bm + i, emptycols[k])]), LinearAlgebra.I)
+  for I in 1:min(blocksize(A)...)
+    bI = Block(I, I)
+    if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+      block = @view!(A[bI])
+      block_alg = block_algorithm(alg, block)
+      bL, bQ = lq_full!(block, block_alg)
+      L[bI] = bL
+      Q[bI] = bQ
+    else
+      copyto!(@view!(Q[bI]), LinearAlgebra.I)
+    end
   end
 
-  return LQ
+  return L, Q
 end

From 4260a998260c0454eec7f004fce4fa0846f72856 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 15:05:33 +0200
Subject: [PATCH 14/18] Refactor eig

---
 src/factorizations/eig.jl | 89 ++++++++++++++++++++++++++++-----------
 1 file changed, 64 insertions(+), 25 deletions(-)

diff --git a/src/factorizations/eig.jl b/src/factorizations/eig.jl
index 502ac10..edfa466 100644
--- a/src/factorizations/eig.jl
+++ b/src/factorizations/eig.jl
@@ -4,7 +4,6 @@ using LinearAlgebra: LinearAlgebra, Diagonal
 using MatrixAlgebraKit:
   MatrixAlgebraKit,
   TruncationStrategy,
-  check_input,
   default_eig_algorithm,
   default_eigh_algorithm,
   diagview,
@@ -26,41 +25,67 @@ for f in [:default_eig_algorithm, :default_eigh_algorithm]
   end
 end
 
+function output_type(::typeof(eig_full!), A::Type{<:AbstractMatrix{T}}) where {T}
+  DV = Base.promote_op(eig_full!, A)
+  return if isconcretetype(DV)
+    DV
+  else
+    Tuple{AbstractMatrix{complex(T)},AbstractMatrix{complex(T)}}
+  end
+end
+function output_type(::typeof(eigh_full!), A::Type{<:AbstractMatrix{T}}) where {T}
+  DV = Base.promote_op(eigh_full!, A)
+  return isconcretetype(DV) ? DV : Tuple{AbstractMatrix{real(T)},AbstractMatrix{T}}
+end
+
+function MatrixAlgebraKit.initialize_output(
+  ::Union{typeof(eig_full!),typeof(eigh_full!)},
+  ::AbstractBlockSparseMatrix,
+  ::BlockPermutedDiagonalAlgorithm,
+)
+  return nothing
+end
+
+function MatrixAlgebraKit.check_input(
+  ::Union{typeof(eig_full!),typeof(eigh_full!)},
+  A::AbstractBlockSparseMatrix,
+  DV,
+  ::BlockPermutedDiagonalAlgorithm,
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
 function MatrixAlgebraKit.check_input(
-  ::typeof(eig_full!), A::AbstractBlockSparseMatrix, (D, V)
+  ::typeof(eig_full!), A::AbstractBlockSparseMatrix, (D, V), ::BlockDiagonalAlgorithm
 )
   @assert isa(D, AbstractBlockSparseMatrix) && isa(V, AbstractBlockSparseMatrix)
   @assert eltype(V) === eltype(D) === complex(eltype(A))
   @assert axes(A, 1) == axes(A, 2)
   @assert axes(A) == axes(D) == axes(V)
+  @assert isblockdiagonal(A)
   return nothing
 end
 function MatrixAlgebraKit.check_input(
-  ::typeof(eigh_full!), A::AbstractBlockSparseMatrix, (D, V)
+  ::typeof(eigh_full!), A::AbstractBlockSparseMatrix, (D, V), ::BlockDiagonalAlgorithm
 )
   @assert isa(D, AbstractBlockSparseMatrix) && isa(V, AbstractBlockSparseMatrix)
   @assert eltype(V) === eltype(A)
   @assert eltype(D) === real(eltype(A))
   @assert axes(A, 1) == axes(A, 2)
   @assert axes(A) == axes(D) == axes(V)
+  @assert isblockdiagonal(A)
   return nothing
 end
 
-function output_type(f::typeof(eig_full!), A::Type{<:AbstractMatrix{T}}) where {T}
-  DV = Base.promote_op(f, A)
-  !isconcretetype(DV) && return Tuple{AbstractMatrix{complex(T)},AbstractMatrix{complex(T)}}
-  return DV
-end
-function output_type(f::typeof(eigh_full!), A::Type{<:AbstractMatrix{T}}) where {T}
-  DV = Base.promote_op(f, A)
-  !isconcretetype(DV) && return Tuple{AbstractMatrix{real(T)},AbstractMatrix{T}}
-  return DV
-end
-
 for f in [:eig_full!, :eigh_full!]
   @eval begin
     function MatrixAlgebraKit.initialize_output(
       ::typeof($f), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+    )
+      return nothing
+    end
+    function MatrixAlgebraKit.initialize_output(
+      ::typeof($f), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
     )
       Td, Tv = fieldtypes(output_type($f, blocktype(A)))
       D = similar(A, BlockType(Td))
@@ -68,18 +93,32 @@ for f in [:eig_full!, :eigh_full!]
       return (D, V)
     end
     function MatrixAlgebraKit.$f(
-      A::AbstractBlockSparseMatrix, (D, V), alg::BlockPermutedDiagonalAlgorithm
+      A::AbstractBlockSparseMatrix, DV, alg::BlockPermutedDiagonalAlgorithm
     )
-      check_input($f, A, (D, V))
-      for I in eachstoredblockdiagindex(A)
-        block = @view!(A[I])
-        block_alg = block_algorithm(alg, block)
-        D[I], V[I] = $f(block, block_alg)
-      end
-      for I in eachunstoredblockdiagindex(A)
-        # TODO: Support setting `LinearAlgebra.I` directly, and/or
-        # using `FillArrays.Eye`.
-        V[I] = LinearAlgebra.I(size(@view(V[I]), 1))
+      MatrixAlgebraKit.check_input($f, A, DV, alg)
+      Ad, transform_rows, transform_cols = blockdiagonalize(A)
+      Dd, Vd = $f(Ad, BlockDiagonalAlgorithm(alg))
+      D = transform_rows(Dd)
+      V = transform_cols(Vd)
+      return D, V
+    end
+    function MatrixAlgebraKit.$f(
+      A::AbstractBlockSparseMatrix, (D, V), alg::BlockDiagonalAlgorithm
+    )
+      MatrixAlgebraKit.check_input($f, A, (D, V), alg)
+
+      # do decomposition on each block
+      for I in 1:min(blocksize(A)...)
+        bI = Block(I, I)
+        if isstored(blocks(A), CartesianIndex(I, I)) # TODO: isblockstored
+          block = @view!(A[bI])
+          block_alg = block_algorithm(alg, block)
+          bD, bV = $f(block, block_alg)
+          D[bI] = bD
+          V[bI] = bV
+        else
+          copyto!(@view!(V[bI]), LinearAlgebra.I)
+        end
       end
       return (D, V)
     end

From b175dfc1b82c43c3ba92002561e50f3864cd6134 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 15:40:08 +0200
Subject: [PATCH 15/18] add eigvals and svdvals, refactor transformers

---
 src/factorizations/eig.jl     | 55 ++++++++++++++++--------
 src/factorizations/lq.jl      | 12 +++---
 src/factorizations/qr.jl      | 12 +++---
 src/factorizations/svd.jl     | 78 ++++++++++++++++++++++++++++++-----
 src/factorizations/utility.jl |  9 ++--
 5 files changed, 122 insertions(+), 44 deletions(-)

diff --git a/src/factorizations/eig.jl b/src/factorizations/eig.jl
index edfa466..ce1288b 100644
--- a/src/factorizations/eig.jl
+++ b/src/factorizations/eig.jl
@@ -1,19 +1,9 @@
 using BlockArrays: blocksizes
 using DiagonalArrays: diagonal
 using LinearAlgebra: LinearAlgebra, Diagonal
-using MatrixAlgebraKit:
-  MatrixAlgebraKit,
-  TruncationStrategy,
-  default_eig_algorithm,
-  default_eigh_algorithm,
-  diagview,
-  eig_full!,
-  eig_trunc!,
-  eig_vals!,
-  eigh_full!,
-  eigh_trunc!,
-  eigh_vals!,
-  findtruncated
+using MatrixAlgebraKit: MatrixAlgebraKit, diagview
+using MatrixAlgebraKit: default_eig_algorithm, eig_full!, eig_vals!
+using MatrixAlgebraKit: default_eigh_algorithm, eigh_full!, eigh_vals!
 
 for f in [:default_eig_algorithm, :default_eigh_algorithm]
   @eval begin
@@ -96,10 +86,10 @@ for f in [:eig_full!, :eigh_full!]
       A::AbstractBlockSparseMatrix, DV, alg::BlockPermutedDiagonalAlgorithm
     )
       MatrixAlgebraKit.check_input($f, A, DV, alg)
-      Ad, transform_rows, transform_cols = blockdiagonalize(A)
+      Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
       Dd, Vd = $f(Ad, BlockDiagonalAlgorithm(alg))
-      D = transform_rows(Dd)
-      V = transform_cols(Vd)
+      D = transform_rows(Dd, invrowperm)
+      V = transform_cols(Vd, invcolperm)
       return D, V
     end
     function MatrixAlgebraKit.$f(
@@ -140,16 +130,47 @@ for f in [:eig_vals!, :eigh_vals!]
   @eval begin
     function MatrixAlgebraKit.initialize_output(
       ::typeof($f), A::AbstractBlockSparseMatrix, alg::BlockPermutedDiagonalAlgorithm
+    )
+      return nothing
+    end
+    function MatrixAlgebraKit.initialize_output(
+      ::typeof($f), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
     )
       T = output_type($f, blocktype(A))
       return similar(A, BlockType(T), axes(A, 1))
     end
+    function MatrixAlgebraKit.check_input(
+      ::typeof($f), A::AbstractBlockSparseMatrix, D, ::BlockPermutedDiagonalAlgorithm
+    )
+      @assert isblockpermuteddiagonal(A)
+      return nothing
+    end
+    function MatrixAlgebraKit.check_input(
+      ::typeof($f), A::AbstractBlockSparseMatrix, D, ::BlockDiagonalAlgorithm
+    )
+      @assert isa(D, AbstractBlockSparseVector)
+      @assert eltype(D) === $(f == :eig_vals! ? complex : real)(eltype(A))
+      @assert axes(A, 1) == axes(A, 2)
+      @assert (axes(A, 1),) == axes(D)
+      @assert isblockdiagonal(A)
+      return nothing
+    end
+
     function MatrixAlgebraKit.$f(
       A::AbstractBlockSparseMatrix, D, alg::BlockPermutedDiagonalAlgorithm
     )
+      MatrixAlgebraKit.check_input($f, A, D, alg)
+      Ad, (invrowperm, _) = blockdiagonalize(A)
+      Dd = $f(Ad, BlockDiagonalAlgorithm(alg))
+      return transform_rows(Dd, invrowperm)
+    end
+    function MatrixAlgebraKit.$f(
+      A::AbstractBlockSparseMatrix, D, alg::BlockDiagonalAlgorithm
+    )
+      MatrixAlgebraKit.check_input($f, A, D, alg)
       for I in eachblockstoredindex(A)
         block = @view!(A[I])
-        D[I] = $f(block, block_algorithm(alg, block))
+        D[Tuple(I)[1]] = $f(block, block_algorithm(alg, block))
       end
       return D
     end
diff --git a/src/factorizations/lq.jl b/src/factorizations/lq.jl
index 306849a..1f5bed4 100644
--- a/src/factorizations/lq.jl
+++ b/src/factorizations/lq.jl
@@ -87,10 +87,10 @@ function MatrixAlgebraKit.lq_compact!(
   A::AbstractBlockSparseMatrix, LQ, alg::BlockPermutedDiagonalAlgorithm
 )
   MatrixAlgebraKit.check_input(lq_compact!, A, LQ, alg)
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Ld, Qd = lq_compact!(Ad, BlockDiagonalAlgorithm(alg))
-  L = transform_rows(Ld)
-  Q = transform_cols(Qd)
+  L = transform_rows(Ld, invrowperm)
+  Q = transform_cols(Qd, invcolperm)
   return L, Q
 end
 function MatrixAlgebraKit.lq_compact!(
@@ -119,10 +119,10 @@ function MatrixAlgebraKit.lq_full!(
   A::AbstractBlockSparseMatrix, LQ, alg::BlockPermutedDiagonalAlgorithm
 )
   MatrixAlgebraKit.check_input(lq_full!, A, LQ, alg)
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Ld, Qd = lq_full!(Ad, BlockDiagonalAlgorithm(alg))
-  L = transform_rows(Ld)
-  Q = transform_cols(Qd)
+  L = transform_rows(Ld, invrowperm)
+  Q = transform_cols(Qd, invcolperm)
   return L, Q
 end
 function MatrixAlgebraKit.lq_full!(
diff --git a/src/factorizations/qr.jl b/src/factorizations/qr.jl
index 99804aa..a564c1e 100644
--- a/src/factorizations/qr.jl
+++ b/src/factorizations/qr.jl
@@ -88,10 +88,10 @@ function MatrixAlgebraKit.qr_compact!(
   A::AbstractBlockSparseMatrix, QR, alg::BlockPermutedDiagonalAlgorithm
 )
   check_input(qr_compact!, A, QR, alg)
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Qd, Rd = qr_compact!(Ad, BlockDiagonalAlgorithm(alg))
-  Q = transform_rows(Qd)
-  R = transform_cols(Rd)
+  Q = transform_rows(Qd, invrowperm)
+  R = transform_cols(Rd, invcolperm)
   return Q, R
 end
 
@@ -121,10 +121,10 @@ function MatrixAlgebraKit.qr_full!(
   A::AbstractBlockSparseMatrix, QR, alg::BlockPermutedDiagonalAlgorithm
 )
   check_input(qr_full!, A, QR, alg)
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Qd, Rd = qr_full!(Ad, BlockDiagonalAlgorithm(alg))
-  Q = transform_rows(Qd)
-  R = transform_cols(Rd)
+  Q = transform_rows(Qd, invrowperm)
+  R = transform_cols(Rd, invcolperm)
   return Q, R
 end
 
diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 6c5f6b2..cb8ac14 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -1,6 +1,6 @@
 using DiagonalArrays: diagonaltype
 using MatrixAlgebraKit:
-  MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!
+  MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!, svd_vals!
 using TypeParameterAccessors: realtype
 
 function MatrixAlgebraKit.default_svd_algorithm(
@@ -15,7 +15,15 @@ function output_type(
   f::Union{typeof(svd_compact!),typeof(svd_full!)}, A::Type{<:AbstractMatrix{T}}
 ) where {T}
   USVᴴ = Base.promote_op(f, A)
-  return isconcretetype(USVᴴ) ? USVᴴ : Tuple{AbstractMatrix{T},AbstractMatrix{realtype(T)},AbstractMatrix{T}}
+  return if isconcretetype(USVᴴ)
+    USVᴴ
+  else
+    Tuple{AbstractMatrix{T},AbstractMatrix{realtype(T)},AbstractMatrix{T}}
+  end
+end
+function output_type(::typeof(svd_vals!), A::Type{<:AbstractMatrix{T}}) where {T}
+  S = Base.promote_op(svd_vals!, A)
+  return isconcretetype(S) ? S : AbstractVector{real(T)}
 end
 
 function MatrixAlgebraKit.initialize_output(
@@ -46,7 +54,6 @@ function MatrixAlgebraKit.initialize_output(
 )
   return nothing
 end
-
 function MatrixAlgebraKit.initialize_output(
   ::typeof(svd_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
 )
@@ -58,6 +65,24 @@ function MatrixAlgebraKit.initialize_output(
   return U, S, Vᴴ
 end
 
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(svd_vals!), ::AbstractBlockSparseMatrix, ::BlockDiagonalAlgorithm
+)
+  return nothing
+end
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(svd_vals!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm
+)
+  brows = eachblockaxis(axes(A, 1))
+  bcols = eachblockaxis(axes(A, 2))
+  # using the property that zip stops as soon as one of the iterators is exhausted
+  s_axes = map(splat(infimum), zip(brows, bcols))
+  s_axis = mortar_axis(s_axes)
+
+  BS = output_type(svd_vals!, blocktype(A))
+  return similar(A, BlockType(BS), S_axes)
+end
+
 function MatrixAlgebraKit.check_input(
   ::typeof(svd_compact!),
   A::AbstractBlockSparseMatrix,
@@ -66,7 +91,6 @@ function MatrixAlgebraKit.check_input(
 )
   @assert isblockpermuteddiagonal(A)
 end
-
 function MatrixAlgebraKit.check_input(
   ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ), ::BlockDiagonalAlgorithm
 )
@@ -87,7 +111,6 @@ function MatrixAlgebraKit.check_input(
   @assert isblockpermuteddiagonal(A)
   return nothing
 end
-
 function MatrixAlgebraKit.check_input(
   ::typeof(svd_full!), A::AbstractBlockSparseMatrix, (U, S, Vᴴ), ::BlockDiagonalAlgorithm
 )
@@ -102,15 +125,30 @@ function MatrixAlgebraKit.check_input(
   return nothing
 end
 
+function MatrixAlgebraKit.check_input(
+  ::typeof(svd_vals!), A::AbstractBlockSparseMatrix, S, ::BlockPermutedDiagonalAlgorithm
+)
+  @assert isblockpermuteddiagonal(A)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(svd_vals!), A::AbstractBlockSparseMatrix, S, ::BlockDiagonalAlgorithm
+)
+  @assert isa(S, AbstractBlockSparseVector)
+  @assert real(eltype(A)) == eltype(S)
+  @assert isblockdiagonal(A)
+  return nothing
+end
+
 function MatrixAlgebraKit.svd_compact!(
   A::AbstractBlockSparseMatrix, USVᴴ, alg::BlockPermutedDiagonalAlgorithm
 )
   check_input(svd_compact!, A, USVᴴ, alg)
 
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Ud, S, Vᴴd = svd_compact!(Ad, BlockDiagonalAlgorithm(alg))
-  U = transform_rows(Ud)
-  Vᴴ = transform_cols(Vᴴd)
+  U = transform_rows(Ud, invrowperm)
+  Vᴴ = transform_cols(Vᴴd, invcolperm)
 
   return U, S, Vᴴ
 end
@@ -143,10 +181,10 @@ function MatrixAlgebraKit.svd_full!(
 )
   check_input(svd_full!, A, USVᴴ, alg)
 
-  Ad, transform_rows, transform_cols = blockdiagonalize(A)
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
   Ud, S, Vᴴd = svd_full!(Ad, BlockDiagonalAlgorithm(alg))
-  U = transform_rows(Ud)
-  Vᴴ = transform_cols(Vᴴd)
+  U = transform_rows(Ud, invrowperm)
+  Vᴴ = transform_cols(Vᴴd, invcolperm)
 
   return U, S, Vᴴ
 end
@@ -181,3 +219,21 @@ function MatrixAlgebraKit.svd_full!(
 
   return U, S, Vᴴ
 end
+
+function MatrixAlgebraKit.svd_vals!(
+  A::AbstractBlockSparseMatrix, S, alg::BlockPermutedDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(svd_vals!, A, S, alg)
+  Ad, _ = blockdiagonalize(A)
+  return svd_vals!(Ad, BlockDiagonalAlgorithm(alg))
+end
+function MatrixAlgebraKit.svd_vals!(
+  A::AbstractBlockSparseMatrix, S, alg::BlockDiagonalAlgorithm
+)
+  MatrixAlgebraKit.check_input(svd_vals!, A, S, alg)
+  for I in eachblockstoredindex(A)
+    block = @view!(A[I])
+    S[Tuple(I)[1]] = $f(block, block_algorithm(alg, block))
+  end
+  return S
+end
diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl
index f56b44b..b269c0f 100644
--- a/src/factorizations/utility.jl
+++ b/src/factorizations/utility.jl
@@ -20,6 +20,10 @@ function supremum(r1::AbstractUnitRange, r2::AbstractUnitRange)
   end
 end
 
+transform_rows(A::AbstractMatrix, invrowperm) = A[invrowperm, :]
+transform_rows(A::AbstractVector, invrowperm) = A[invrowperm]
+transform_cols(A::AbstractMatrix, invcolperm) = A[:, invcolperm]
+
 function blockdiagonalize(A::AbstractBlockSparseMatrix)
   # sort in order to avoid depending on internal details such as dictionary order
   bIs = sort!(collect(eachblockstoredindex(A)); by=Int ∘ last ∘ Tuple)
@@ -41,12 +45,9 @@ function blockdiagonalize(A::AbstractBlockSparseMatrix)
   append!(colperm, emptycols)
 
   invrowperm = Block.(invperm(Int.(rowperm)))
-  transform_rows(A) = A[invrowperm, :]
-
   invcolperm = Block.(invperm(Int.(colperm)))
-  transform_cols(A) = A[:, invcolperm]
 
-  return A[rowperm, colperm], transform_rows, transform_cols
+  return A[rowperm, colperm], (invrowperm, invcolperm)
 end
 
 function isblockdiagonal(A::AbstractBlockSparseMatrix)

From a992d38be6fecc362a28e3ca630eeb5808e59722 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 16:00:56 +0200
Subject: [PATCH 16/18] Also fix factorizations

---
 src/factorizations/truncation.jl | 58 +++++++++++++++++++++-----------
 1 file changed, 38 insertions(+), 20 deletions(-)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index d173f88..a28995b 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -1,4 +1,6 @@
 using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  TruncatedAlgorithm,
   TruncationStrategy,
   diagview,
   eig_trunc!,
@@ -8,34 +10,50 @@ using MatrixAlgebraKit:
   truncate!
 
 """
-    BlockPermutedDiagonalTruncationStrategy(strategy::TruncationStrategy)
+    BlockDiagonalTruncationStrategy(strategy::TruncationStrategy)
 
 A wrapper for `TruncationStrategy` that implements the wrapped strategy on a block-by-block
-basis, which is possible if the input matrix is a block-diagonal matrix or a block permuted
-block-diagonal matrix.
+basis, which is possible if the input matrix is a block-diagonal matrix.
 """
-struct BlockPermutedDiagonalTruncationStrategy{T<:TruncationStrategy} <: TruncationStrategy
+struct BlockDiagonalTruncationStrategy{T<:TruncationStrategy} <: TruncationStrategy
   strategy::T
 end
 
-function MatrixAlgebraKit.truncate!(
-  ::typeof(svd_trunc!),
-  (U, S, Vᴴ)::NTuple{3,AbstractBlockSparseMatrix},
-  strategy::TruncationStrategy,
+function BlockDiagonalTruncationStrategy(alg::BlockPermutedDiagonalAlgorithm)
+  return BlockDiagonalTruncationStrategy(alg.strategy)
+end
+
+function MatrixAlgebraKit.svd_trunc!(
+  A::AbstractBlockSparseMatrix,
+  out,
+  alg::TruncatedAlgorithm{<:BlockPermutedDiagonalAlgorithm},
 )
-  # TODO assert blockdiagonal
-  return truncate!(
-    svd_trunc!, (U, S, Vᴴ), BlockPermutedDiagonalTruncationStrategy(strategy)
-  )
+  Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
+  blockalg = BlockDiagonalAlgorithm(alg.alg)
+  blockstrategy = BlockDiagonalTruncationStrategy(alg.trunc)
+  Ud, S, Vᴴd = svd_trunc!(Ad, TruncatedAlgorithm(blockalg, blockstrategy))
+
+  U = transform_rows(Ud, invrowperm)
+  Vᴴ = transform_cols(Vᴴd, invcolperm)
+
+  return U, S, Vᴴ
 end
+
 for f in [:eig_trunc!, :eigh_trunc!]
   @eval begin
-    function MatrixAlgebraKit.truncate!(
-      ::typeof($f),
-      (D, V)::NTuple{2,AbstractBlockSparseMatrix},
-      strategy::TruncationStrategy,
+    function MatrixAlgebraKit.$f(
+      A::AbstractBlockSparseMatrix,
+      out,
+      alg::TruncatedAlgorithm{<:BlockPermutedDiagonalAlgorithm},
     )
-      return truncate!($f, (D, V), BlockPermutedDiagonalTruncationStrategy(strategy))
+      Ad, (invrowperm, invcolperm) = blockdiagonalize(A)
+      blockalg = BlockDiagonalAlgorithm(alg.alg)
+      blockstrategy = BlockDiagonalTruncationStrategy(alg.trunc)
+      Dd, Vd = $f(Ad, TruncatedAlgorithm(blockalg, blockstrategy))
+
+      D = transform_rows(Dd, invrowperm)
+      V = transform_cols(Vd, invcolperm)
+      return D, V
     end
   end
 end
@@ -43,7 +61,7 @@ end
 # cannot use regular slicing here: I want to slice without altering blockstructure
 # solution: use boolean indexing and slice the mask, effectively cheaply inverting the map
 function MatrixAlgebraKit.findtruncated(
-  values::AbstractVector, strategy::BlockPermutedDiagonalTruncationStrategy
+  values::AbstractVector, strategy::BlockDiagonalTruncationStrategy
 )
   ind = findtruncated(Vector(values), strategy.strategy)
   indexmask = falses(length(values))
@@ -66,7 +84,7 @@ end
 function MatrixAlgebraKit.truncate!(
   ::typeof(svd_trunc!),
   (U, S, Vᴴ)::NTuple{3,AbstractBlockSparseMatrix},
-  strategy::BlockPermutedDiagonalTruncationStrategy,
+  strategy::BlockDiagonalTruncationStrategy,
 )
   I = findtruncated(diag(S), strategy)
   return (U[:, I], S[I, I], Vᴴ[I, :])
@@ -76,7 +94,7 @@ for f in [:eig_trunc!, :eigh_trunc!]
     function MatrixAlgebraKit.truncate!(
       ::typeof($f),
       (D, V)::NTuple{2,AbstractBlockSparseMatrix},
-      strategy::BlockPermutedDiagonalTruncationStrategy,
+      strategy::BlockDiagonalTruncationStrategy,
     )
       I = findtruncated(diag(D), strategy)
       return (D[I, I], V[:, I])

From ec0b9733487e9b14207fb6182e4a2d17dc6d5368 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 16:01:45 +0200
Subject: [PATCH 17/18] Bump v0.10

---
 Project.toml          | 2 +-
 docs/Project.toml     | 2 +-
 examples/Project.toml | 2 +-
 test/Project.toml     | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0ca7f6c..1905cdc 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.9.1"
+version = "0.10.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/docs/Project.toml b/docs/Project.toml
index edbaa0b..6f5417f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.9"
+BlockSparseArrays = "0.10"
 Documenter = "1"
 Literate = "2"
diff --git a/examples/Project.toml b/examples/Project.toml
index 7e02e0d..14c5df3 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.9"
+BlockSparseArrays = "0.10"
 Test = "1"
diff --git a/test/Project.toml b/test/Project.toml
index 12fb6f0..c084b9f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -25,7 +25,7 @@ Adapt = "4"
 Aqua = "0.8"
 ArrayLayouts = "1"
 BlockArrays = "1"
-BlockSparseArrays = "0.9"
+BlockSparseArrays = "0.10"
 DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3"
 GPUArraysCore = "0.2"

From ca1350c372974147432e48f9e0e02d2b37078be2 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 15 Aug 2025 16:05:36 +0200
Subject: [PATCH 18/18] small fixes

---
 src/factorizations/svd.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index cb8ac14..c2482f7 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -66,7 +66,7 @@ function MatrixAlgebraKit.initialize_output(
 end
 
 function MatrixAlgebraKit.initialize_output(
-  ::typeof(svd_vals!), ::AbstractBlockSparseMatrix, ::BlockDiagonalAlgorithm
+  ::typeof(svd_vals!), ::AbstractBlockSparseMatrix, ::BlockPermutedDiagonalAlgorithm
 )
   return nothing
 end
@@ -233,7 +233,7 @@ function MatrixAlgebraKit.svd_vals!(
   MatrixAlgebraKit.check_input(svd_vals!, A, S, alg)
   for I in eachblockstoredindex(A)
     block = @view!(A[I])
-    S[Tuple(I)[1]] = $f(block, block_algorithm(alg, block))
+    S[Tuple(I)[1]] = svd_vals!(block, block_algorithm(alg, block))
   end
   return S
 end