From 150b8f5efb6076ad7139701b8a2a83e47be9902d Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 17 Dec 2024 09:53:54 +0100
Subject: [PATCH 01/21] Migrate changes from
 https://github.com/ITensor/ITensors.jl/pull/1572

---
 .../BlockArraysExtensions.jl                  |  25 +++
 src/BlockSparseArrays.jl                      |  12 +
 .../abstractblocksparsematrix.jl              |  21 ++
 src/blocksparsearray/blockdiagonalarray.jl    |  59 +++++
 src/factorizations/svd.jl                     | 206 ++++++++++++++++++
 test/basics/test_svd.jl                       | 164 ++++++++++++++
 6 files changed, 487 insertions(+)
 create mode 100644 src/blocksparsearray/blockdiagonalarray.jl
 create mode 100644 src/factorizations/svd.jl
 create mode 100644 test/basics/test_svd.jl

diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl
index 3208508..e0a3c7a 100644
--- a/src/BlockArraysExtensions/BlockArraysExtensions.jl
+++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl
@@ -598,3 +598,28 @@ macro view!(expr)
   @capture(expr, array_[indices__])
   return :(view!($(esc(array)), $(esc.(indices)...)))
 end
+
+# SVD additions
+# -------------
+using LinearAlgebra: Algorithm
+using BlockArrays: BlockedMatrix
+
+# svd first calls `eigencopy_oftype` to create something that can be in-place SVD'd
+# Here, we hijack this system to determine if there is any structure we can exploit
+# default: SVD is most efficient with BlockedArray
+function eigencopy_oftype(A::AbstractBlockArray, S)
+  return BlockedMatrix{S}(A)
+end
+
+function svd!(A::BlockedMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A))
+  F = svd!(parent(A); full, alg)
+
+  # restore block pattern
+  m = length(F.S)
+  bax1, bax2, bax3 = axes(A, 1), blockedrange([m]), axes(A, 2)
+
+  u = BlockedArray(F.U, (bax1, bax2))
+  s = BlockedVector(F.S, (bax2,))
+  vt = BlockedArray(F.Vt, (bax2, bax3))
+  return SVD(u, s, vt)
+end
diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index 1425f1d..f125e2b 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -7,7 +7,13 @@ export BlockSparseArray,
   eachblockstoredindex,
   eachstoredblock
 
+# factorizations
+include("factorizations/svd.jl")
+
+# possible upstream contributions
 include("BlockArraysExtensions/BlockArraysExtensions.jl")
+
+# interface functions that don't have to specialize
 include("blocksparsearrayinterface/blocksparsearrayinterface.jl")
 include("blocksparsearrayinterface/linearalgebra.jl")
 include("blocksparsearrayinterface/getunstoredblock.jl")
@@ -16,6 +22,8 @@ include("blocksparsearrayinterface/map.jl")
 include("blocksparsearrayinterface/arraylayouts.jl")
 include("blocksparsearrayinterface/views.jl")
 include("blocksparsearrayinterface/cat.jl")
+
+# functions defined for any abstractblocksparsearray
 include("abstractblocksparsearray/abstractblocksparsearray.jl")
 include("abstractblocksparsearray/wrappedabstractblocksparsearray.jl")
 include("abstractblocksparsearray/abstractblocksparsematrix.jl")
@@ -27,8 +35,12 @@ include("abstractblocksparsearray/broadcast.jl")
 include("abstractblocksparsearray/map.jl")
 include("abstractblocksparsearray/linearalgebra.jl")
 include("abstractblocksparsearray/cat.jl")
+
+# functions specifically for BlockSparseArray
 include("blocksparsearray/defaults.jl")
 include("blocksparsearray/blocksparsearray.jl")
+include("blocksparsearray/blockdiagonalarray.jl")
+
 include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl")
 
 end
diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 0c2c578..cb5e374 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -1 +1,22 @@
 const AbstractBlockSparseMatrix{T} = AbstractBlockSparseArray{T,2}
+
+# SVD is implemented by trying to
+# 1. Attempt to find a block-diagonal implementation by permuting
+# 2. Fallback to AbstractBlockArray implementation via BlockedArray
+function svd(
+  A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A)
+)
+  T = LinearAlgebra.eigtype(eltype(A))
+  A′ = try_to_blockdiagonal(A)
+
+  if isnothing(A′)
+    # not block-diagonal, fall back to dense case
+    Adense = eigencopy_oftype(A, T)
+    return svd!(Adense; full, alg)
+  end
+
+  # compute block-by-block and permute back
+  A″, (I, J) = A′
+  F = svd!(eigencopy_oftype(A″, T); full, alg)
+  return SVD(F.U[Block.(I), Block.(J)], F.S, F.Vt)
+end
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
new file mode 100644
index 0000000..b900101
--- /dev/null
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -0,0 +1,59 @@
+# type alias for block-diagonal
+using LinearAlgebra: Diagonal
+
+const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{
+  T,A,Diagonal{A,V},Axes
+}
+
+function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
+  return BlockSparseArray(
+    Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
+  )
+end
+
+# Cast to block-diagonal implementation if permuted-blockdiagonal
+function try_to_blockdiagonal_perm(A)
+  inds = map(x -> Int.(Tuple(x)), vec(collect(block_stored_indices(A))))
+  I = first.(inds)
+  allunique(I) || return nothing
+  J = last.(inds)
+  p = sortperm(J)
+  Jsorted = J[p]
+  allunique(Jsorted) || return nothing
+  return Block.(I[p], Jsorted)
+end
+
+"""
+    try_to_blockdiagonal(A)
+
+Attempt to find a permutation of blocks that makes `A` blockdiagonal. If unsuccesful,
+returns nothing, otherwise returns both the blockdiagonal `B` as well as the permutation `I, J`.
+"""
+function try_to_blockdiagonal(A::AbstractBlockSparseMatrix)
+  perm = try_to_blockdiagonal_perm(A)
+  isnothing(perm) && return perm
+  I = first.(Tuple.(perm))
+  J = last.(Tuple.(perm))
+  diagblocks = map(invperm(I), J) do i, j
+    return A[Block(i, j)]
+  end
+  return BlockDiagonal(diagblocks), perm
+end
+
+# SVD implementation
+function eigencopy_oftype(A::BlockDiagonal, S)
+  diag = map(Base.Fix2(eigencopy_oftype, S), A.blocks.diag)
+  return BlockDiagonal(diag)
+end
+
+function svd(A::BlockDiagonal; kwargs...)
+  return svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...)
+end
+function svd!(A::BlockDiagonal; full::Bool=false, alg::Algorithm=default_svd_alg(A))
+  # TODO: handle full
+  F = map(a -> svd!(a; full, alg), blocks(A).diag)
+  Us = map(Base.Fix2(getproperty, :U), F)
+  Ss = map(Base.Fix2(getproperty, :S), F)
+  Vts = map(Base.Fix2(getproperty, :Vt), F)
+  return SVD(BlockDiagonal(Us), mortar(Ss), BlockDiagonal(Vts))
+end
diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
new file mode 100644
index 0000000..9e46dc6
--- /dev/null
+++ b/src/factorizations/svd.jl
@@ -0,0 +1,206 @@
+using LinearAlgebra:
+  LinearAlgebra, Factorization, Algorithm, default_svd_alg, Adjoint, Transpose
+using BlockArrays: AbstractBlockMatrix, BlockedArray, BlockedMatrix, BlockedVector
+using BlockArrays: BlockLayout
+
+# Singular Value Decomposition:
+# need new type to deal with U and V having possible different types
+# this is basically a carbon copy of the LinearAlgebra implementation.
+# additionally, by default we implement a fallback to the LinearAlgebra implementation
+# in hope to support as many foreign types as possible that chose to extend those methods.
+
+# TODO: add this to MatrixFactorizations
+# TODO: decide where this goes
+# TODO: decide whether or not to restrict types to be blocked.
+"""
+    SVD <: Factorization
+
+Matrix factorization type of the singular value decomposition (SVD) of a matrix `A`.
+This is the return type of [`svd(_)`](@ref), the corresponding matrix factorization function.
+
+If `F::SVD` is the factorization object, `U`, `S`, `V` and `Vt` can be obtained
+via `F.U`, `F.S`, `F.V` and `F.Vt`, such that `A = U * Diagonal(S) * Vt`.
+The singular values in `S` are sorted in descending order.
+
+Iterating the decomposition produces the components `U`, `S`, and `V`.
+
+# Examples
+```jldoctest
+julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
+4×5 Matrix{Float64}:
+ 1.0  0.0  0.0  0.0  2.0
+ 0.0  0.0  3.0  0.0  0.0
+ 0.0  0.0  0.0  0.0  0.0
+ 0.0  2.0  0.0  0.0  0.0
+
+julia> F = svd(A)
+SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
+U factor:
+4×4 Matrix{Float64}:
+ 0.0  1.0   0.0  0.0
+ 1.0  0.0   0.0  0.0
+ 0.0  0.0   0.0  1.0
+ 0.0  0.0  -1.0  0.0
+singular values:
+4-element Vector{Float64}:
+ 3.0
+ 2.23606797749979
+ 2.0
+ 0.0
+Vt factor:
+4×5 Matrix{Float64}:
+ -0.0        0.0  1.0  -0.0  0.0
+  0.447214   0.0  0.0   0.0  0.894427
+  0.0       -1.0  0.0   0.0  0.0
+  0.0        0.0  0.0   1.0  0.0
+
+julia> F.U * Diagonal(F.S) * F.Vt
+4×5 Matrix{Float64}:
+ 1.0  0.0  0.0  0.0  2.0
+ 0.0  0.0  3.0  0.0  0.0
+ 0.0  0.0  0.0  0.0  0.0
+ 0.0  2.0  0.0  0.0  0.0
+
+julia> u, s, v = F; # destructuring via iteration
+
+julia> u == F.U && s == F.S && v == F.V
+true
+```
+"""
+struct SVD{T,Tr,M<:AbstractArray{T},C<:AbstractVector{Tr},N<:AbstractArray{T}} <:
+       Factorization{T}
+  U::M
+  S::C
+  Vt::N
+  function SVD{T,Tr,M,C,N}(
+    U, S, Vt
+  ) where {T,Tr,M<:AbstractArray{T},C<:AbstractVector{Tr},N<:AbstractArray{T}}
+    Base.require_one_based_indexing(U, S, Vt)
+    return new{T,Tr,M,C,N}(U, S, Vt)
+  end
+end
+function SVD(U::AbstractArray{T}, S::AbstractVector{Tr}, Vt::AbstractArray{T}) where {T,Tr}
+  return SVD{T,Tr,typeof(U),typeof(S),typeof(Vt)}(U, S, Vt)
+end
+function SVD{T}(U::AbstractArray, S::AbstractVector{Tr}, Vt::AbstractArray) where {T,Tr}
+  return SVD(
+    convert(AbstractArray{T}, U),
+    convert(AbstractVector{Tr}, S),
+    convert(AbstractArray{T}, Vt),
+  )
+end
+
+function SVD{T}(F::SVD) where {T}
+  return SVD(
+    convert(AbstractMatrix{T}, F.U),
+    convert(AbstractVector{real(T)}, F.S),
+    convert(AbstractMatrix{T}, F.Vt),
+  )
+end
+LinearAlgebra.Factorization{T}(F::SVD) where {T} = SVD{T}(F)
+
+# iteration for destructuring into components
+Base.iterate(S::SVD) = (S.U, Val(:S))
+Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V))
+Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done))
+Base.iterate(::SVD, ::Val{:done}) = nothing
+
+function Base.getproperty(F::SVD, d::Symbol)
+  if d === :V
+    return getfield(F, :Vt)'
+  else
+    return getfield(F, d)
+  end
+end
+
+function Base.propertynames(F::SVD, private::Bool=false)
+  return private ? (:V, fieldnames(typeof(F))...) : (:U, :S, :V, :Vt)
+end
+
+Base.size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim)
+Base.size(A::SVD) = (size(A, 1), size(A, 2))
+
+function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::SVD)
+  summary(io, F)
+  println(io)
+  println(io, "U factor:")
+  show(io, mime, F.U)
+  println(io, "\nsingular values:")
+  show(io, mime, F.S)
+  println(io, "\nVt factor:")
+  return show(io, mime, F.Vt)
+end
+
+Base.adjoint(usv::SVD) = SVD(adjoint(usv.Vt), usv.S, adjoint(usv.U))
+Base.transpose(usv::SVD) = SVD(transpose(usv.Vt), usv.S, transpose(usv.U))
+
+# Conversion
+Base.AbstractMatrix(F::SVD) = (F.U * Diagonal(F.S)) * F.Vt
+Base.AbstractArray(F::SVD) = AbstractMatrix(F)
+Base.Matrix(F::SVD) = Array(AbstractArray(F))
+Base.Array(F::SVD) = Matrix(F)
+SVD(usv::SVD) = usv
+SVD(usv::LinearAlgebra.SVD) = SVD(usv.U, usv.S, usv.Vt)
+
+# functions default to LinearAlgebra
+# ----------------------------------
+"""
+    svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
+
+`svd!` is the same as [`svd`](@ref), but saves space by
+overwriting the input `A`, instead of creating a copy. See documentation of [`svd`](@ref) for details.
+"""
+svd!(A; kwargs...) = SVD(LinearAlgebra.svd!(A; kwargs...))
+
+"""
+    svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
+
+Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.
+
+`U`, `S`, `V` and `Vt` can be obtained from the factorization `F` with `F.U`,
+`F.S`, `F.V` and `F.Vt`, such that `A = U * Diagonal(S) * Vt`.
+The algorithm produces `Vt` and hence `Vt` is more efficient to extract than `V`.
+The singular values in `S` are sorted in descending order.
+
+Iterating the decomposition produces the components `U`, `S`, and `V`.
+
+If `full = false` (default), a "thin" SVD is returned. For an ``M
+\\times N`` matrix `A`, in the full factorization `U` is ``M \\times M``
+and `V` is ``N \\times N``, while in the thin factorization `U` is ``M
+\\times K`` and `V` is ``N \\times K``, where ``K = \\min(M,N)`` is the
+number of singular values.
+
+`alg` specifies which algorithm and LAPACK method to use for SVD:
+- `alg = DivideAndConquer()` (default): Calls `LAPACK.gesdd!`.
+- `alg = QRIteration()`: Calls `LAPACK.gesvd!` (typically slower but more accurate) .
+
+!!! compat "Julia 1.3"
+    The `alg` keyword argument requires Julia 1.3 or later.
+
+# Examples
+```jldoctest
+julia> A = rand(4,3);
+
+julia> F = svd(A); # Store the Factorization Object
+
+julia> A ≈ F.U * Diagonal(F.S) * F.Vt
+true
+
+julia> U, S, V = F; # destructuring via iteration
+
+julia> A ≈ U * Diagonal(S) * V'
+true
+
+julia> Uonly, = svd(A); # Store U only
+
+julia> Uonly == U
+true
+```
+"""
+svd(A; kwargs...) =
+  SVD(svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...))
+
+LinearAlgebra.svdvals(usv::SVD{<:Any,T}) where {T} = (usv.S)::Vector{T}
+
+# Added here to avoid type-piracy
+eigencopy_oftype(A, S) = LinearAlgebra.eigencopy_oftype(A, S)
diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl
new file mode 100644
index 0000000..db9655d
--- /dev/null
+++ b/test/basics/test_svd.jl
@@ -0,0 +1,164 @@
+using Test
+using BlockSparseArrays
+using BlockSparseArrays:
+    BlockSparseArray, svd, tsvd, notrunc, truncbelow, truncdim, BlockDiagonal
+using BlockArrays
+using LinearAlgebra: LinearAlgebra, Diagonal, svdvals
+
+function test_svd(a, usv)
+    U, S, V = usv
+
+    @test U * Diagonal(S) * V' ≈ a
+    @test U' * U ≈ LinearAlgebra.I
+    @test V' * V ≈ LinearAlgebra.I
+end
+
+# regular matrix
+# --------------
+sizes = ((3, 3), (4, 3), (3, 4))
+eltypes = (Float32, Float64, ComplexF64)
+@testset "($m, $n) Matrix{$T}" for ((m, n), T) in Iterators.product(sizes, eltypes)
+    a = rand(m, n)
+    usv = @inferred svd(a)
+    test_svd(a, usv)
+
+    # TODO: type unstable?
+    usv2 = tsvd(a)
+    test_svd(a, usv2)
+
+    usv3 = tsvd(a; trunc=truncdim(2))
+    @test length(usv3.S) == 2
+    @test usv3.U' * usv3.U ≈ LinearAlgebra.I
+    @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
+
+    s = usv3.S[end]
+    usv4 = tsvd(a; trunc=truncbelow(s))
+    @test length(usv4.S) == 2
+    @test usv4.U' * usv4.U ≈ LinearAlgebra.I
+    @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
+end
+
+# block matrix
+# ------------
+blockszs = (([2, 2], [2, 2]), ([2, 2], [2, 3]), ([2, 2, 1], [2, 3]), ([2, 3], [2]))
+@testset "($m, $n) BlockMatrix{$T}" for ((m, n), T) in Iterators.product(blockszs, eltypes)
+    a = mortar([rand(T, i, j) for i in m, j in n])
+    usv = svd(a)
+    test_svd(a, usv)
+    @test usv.U isa BlockedMatrix
+    @test usv.Vt isa BlockedMatrix
+    @test usv.S isa BlockedVector
+
+    usv2 = tsvd(a)
+    test_svd(a, usv2)
+    @test usv.U isa BlockedMatrix
+    @test usv.Vt isa BlockedMatrix
+    @test usv.S isa BlockedVector
+
+    usv3 = tsvd(a; trunc=truncdim(2))
+    @test length(usv3.S) == 2
+    @test usv3.U' * usv3.U ≈ LinearAlgebra.I
+    @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
+    @test usv.U isa BlockedMatrix
+    @test usv.Vt isa BlockedMatrix
+    @test usv.S isa BlockedVector
+
+    s = usv3.S[end]
+    usv4 = tsvd(a; trunc=truncbelow(s))
+    @test length(usv4.S) == 2
+    @test usv4.U' * usv4.U ≈ LinearAlgebra.I
+    @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
+    @test usv.U isa BlockedMatrix
+    @test usv.Vt isa BlockedMatrix
+    @test usv.S isa BlockedVector
+end
+
+# Block-Diagonal matrices
+# -----------------------
+@testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in
+                                          Iterators.product(blockszs, eltypes)
+    a = BlockDiagonal([rand(T, i, j) for (i, j) in zip(m, n)])
+    usv = svd(a)
+    # TODO: `BlockDiagonal * Adjoint` errors
+    test_svd(a, usv)
+    @test usv.U isa BlockDiagonal
+    @test usv.Vt isa BlockDiagonal
+    @test usv.S isa BlockVector
+
+    usv2 = tsvd(a)
+    test_svd(a, usv2)
+    @test usv.U isa BlockDiagonal
+    @test usv.Vt isa BlockDiagonal
+    @test usv.S isa BlockVector
+
+    # TODO: need to find a slicing fix to make this work
+    # usv3 = tsvd(a; trunc=truncdim(2))
+    # @test length(usv3.S) == 2
+    # @test usv3.U' * usv3.U ≈ LinearAlgebra.I
+    # @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
+    # @test usv.U isa BlockDiagonal
+    # @test usv.Vt isa BlockDiagonal
+    # @test usv.S isa BlockVector
+
+    # @show s = usv3.S[end]
+    # usv4 = tsvd(a; trunc=truncbelow(s))
+    # @test length(usv4.S) == 2
+    # @test usv4.U' * usv4.U ≈ LinearAlgebra.I
+    # @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
+    # @test usv.U isa BlockDiagonal
+    # @test usv.Vt isa BlockDiagonal
+    # @test usv.S isa BlockVector
+end
+
+a = mortar([rand(2, 2) for i in 1:2, j in 1:3])
+usv = svd(a)
+test_svd(a, usv)
+
+a = mortar([rand(2, 2) for i in 1:3, j in 1:2])
+usv = svd(a)
+test_svd(a, usv)
+
+# blocksparse 
+# -----------
+@testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in
+                                          Iterators.product(blockszs, eltypes)
+    a = BlockSparseArray{T}(m, n)
+    for i in LinearAlgebra.diagind(blocks(a))
+        I = CartesianIndices(blocks(a))[i]
+        a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
+    end
+    perm = Random.randperm(length(m))
+    a = a[Block.(perm), Block.(1:length(n))]
+
+    # errors because `blocks(a)[CartesianIndex.(...)]` is not implemented
+    usv = svd(a)
+    # TODO: `BlockDiagonal * Adjoint` errors
+    # test_svd(a, usv)
+    @test usv.U isa BlockDiagonal
+    @test usv.Vt isa BlockDiagonal
+    @test usv.S isa BlockVector
+
+    # usv2 = tsvd(a)
+    test_svd(a, usv2)
+    @test usv.U isa BlockDiagonal
+    @test usv.Vt isa BlockDiagonal
+    @test usv.S isa BlockVector
+
+    # TODO: need to find a slicing fix to make this work
+    # usv3 = tsvd(a; trunc=truncdim(2))
+    # @test length(usv3.S) == 2
+    # @test usv3.U' * usv3.U ≈ LinearAlgebra.I
+    # @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
+    # @test usv.U isa BlockDiagonal
+    # @test usv.Vt isa BlockDiagonal
+    # @test usv.S isa BlockVector
+
+    # @show s = usv3.S[end]
+    # usv4 = tsvd(a; trunc=truncbelow(s))
+    # @test length(usv4.S) == 2
+    # @test usv4.U' * usv4.U ≈ LinearAlgebra.I
+    # @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
+    # @test usv.U isa BlockDiagonal
+    # @test usv.Vt isa BlockDiagonal
+    # @test usv.S isa BlockVector
+end

From aaf8e697802c652305db122e01cf00229bfd7905 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 17 Dec 2024 10:01:11 +0100
Subject: [PATCH 02/21] remove `tsvd` tests

---
 test/basics/test_svd.jl | 85 +----------------------------------------
 1 file changed, 2 insertions(+), 83 deletions(-)

diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl
index db9655d..832a03e 100644
--- a/test/basics/test_svd.jl
+++ b/test/basics/test_svd.jl
@@ -1,7 +1,7 @@
 using Test
 using BlockSparseArrays
 using BlockSparseArrays:
-    BlockSparseArray, svd, tsvd, notrunc, truncbelow, truncdim, BlockDiagonal
+    BlockSparseArray, svd, notrunc, truncbelow, truncdim, BlockDiagonal
 using BlockArrays
 using LinearAlgebra: LinearAlgebra, Diagonal, svdvals
 
@@ -21,21 +21,6 @@ eltypes = (Float32, Float64, ComplexF64)
     a = rand(m, n)
     usv = @inferred svd(a)
     test_svd(a, usv)
-
-    # TODO: type unstable?
-    usv2 = tsvd(a)
-    test_svd(a, usv2)
-
-    usv3 = tsvd(a; trunc=truncdim(2))
-    @test length(usv3.S) == 2
-    @test usv3.U' * usv3.U ≈ LinearAlgebra.I
-    @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
-
-    s = usv3.S[end]
-    usv4 = tsvd(a; trunc=truncbelow(s))
-    @test length(usv4.S) == 2
-    @test usv4.U' * usv4.U ≈ LinearAlgebra.I
-    @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
 end
 
 # block matrix
@@ -48,29 +33,6 @@ blockszs = (([2, 2], [2, 2]), ([2, 2], [2, 3]), ([2, 2, 1], [2, 3]), ([2, 3], [2
     @test usv.U isa BlockedMatrix
     @test usv.Vt isa BlockedMatrix
     @test usv.S isa BlockedVector
-
-    usv2 = tsvd(a)
-    test_svd(a, usv2)
-    @test usv.U isa BlockedMatrix
-    @test usv.Vt isa BlockedMatrix
-    @test usv.S isa BlockedVector
-
-    usv3 = tsvd(a; trunc=truncdim(2))
-    @test length(usv3.S) == 2
-    @test usv3.U' * usv3.U ≈ LinearAlgebra.I
-    @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
-    @test usv.U isa BlockedMatrix
-    @test usv.Vt isa BlockedMatrix
-    @test usv.S isa BlockedVector
-
-    s = usv3.S[end]
-    usv4 = tsvd(a; trunc=truncbelow(s))
-    @test length(usv4.S) == 2
-    @test usv4.U' * usv4.U ≈ LinearAlgebra.I
-    @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
-    @test usv.U isa BlockedMatrix
-    @test usv.Vt isa BlockedMatrix
-    @test usv.S isa BlockedVector
 end
 
 # Block-Diagonal matrices
@@ -84,30 +46,6 @@ end
     @test usv.U isa BlockDiagonal
     @test usv.Vt isa BlockDiagonal
     @test usv.S isa BlockVector
-
-    usv2 = tsvd(a)
-    test_svd(a, usv2)
-    @test usv.U isa BlockDiagonal
-    @test usv.Vt isa BlockDiagonal
-    @test usv.S isa BlockVector
-
-    # TODO: need to find a slicing fix to make this work
-    # usv3 = tsvd(a; trunc=truncdim(2))
-    # @test length(usv3.S) == 2
-    # @test usv3.U' * usv3.U ≈ LinearAlgebra.I
-    # @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
-    # @test usv.U isa BlockDiagonal
-    # @test usv.Vt isa BlockDiagonal
-    # @test usv.S isa BlockVector
-
-    # @show s = usv3.S[end]
-    # usv4 = tsvd(a; trunc=truncbelow(s))
-    # @test length(usv4.S) == 2
-    # @test usv4.U' * usv4.U ≈ LinearAlgebra.I
-    # @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
-    # @test usv.U isa BlockDiagonal
-    # @test usv.Vt isa BlockDiagonal
-    # @test usv.S isa BlockVector
 end
 
 a = mortar([rand(2, 2) for i in 1:2, j in 1:3])
@@ -133,32 +71,13 @@ test_svd(a, usv)
     # errors because `blocks(a)[CartesianIndex.(...)]` is not implemented
     usv = svd(a)
     # TODO: `BlockDiagonal * Adjoint` errors
-    # test_svd(a, usv)
+    test_svd(a, usv)
     @test usv.U isa BlockDiagonal
     @test usv.Vt isa BlockDiagonal
     @test usv.S isa BlockVector
 
-    # usv2 = tsvd(a)
     test_svd(a, usv2)
     @test usv.U isa BlockDiagonal
     @test usv.Vt isa BlockDiagonal
     @test usv.S isa BlockVector
-
-    # TODO: need to find a slicing fix to make this work
-    # usv3 = tsvd(a; trunc=truncdim(2))
-    # @test length(usv3.S) == 2
-    # @test usv3.U' * usv3.U ≈ LinearAlgebra.I
-    # @test usv3.Vt * usv3.V ≈ LinearAlgebra.I
-    # @test usv.U isa BlockDiagonal
-    # @test usv.Vt isa BlockDiagonal
-    # @test usv.S isa BlockVector
-
-    # @show s = usv3.S[end]
-    # usv4 = tsvd(a; trunc=truncbelow(s))
-    # @test length(usv4.S) == 2
-    # @test usv4.U' * usv4.U ≈ LinearAlgebra.I
-    # @test usv4.Vt * usv4.V ≈ LinearAlgebra.I
-    # @test usv.U isa BlockDiagonal
-    # @test usv.Vt isa BlockDiagonal
-    # @test usv.S isa BlockVector
 end

From dcd44666cd7f788e4c5c941c675445a993b8ac3c Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 17 Dec 2024 14:32:55 +0100
Subject: [PATCH 03/21] Formatter

---
 test/basics/test_svd.jl | 80 ++++++++++++++++++++---------------------
 1 file changed, 40 insertions(+), 40 deletions(-)

diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl
index 832a03e..00903d4 100644
--- a/test/basics/test_svd.jl
+++ b/test/basics/test_svd.jl
@@ -1,16 +1,16 @@
 using Test
 using BlockSparseArrays
-using BlockSparseArrays:
-    BlockSparseArray, svd, notrunc, truncbelow, truncdim, BlockDiagonal
+using BlockSparseArrays: BlockSparseArray, svd, notrunc, truncbelow, truncdim, BlockDiagonal
 using BlockArrays
 using LinearAlgebra: LinearAlgebra, Diagonal, svdvals
+using Random
 
 function test_svd(a, usv)
-    U, S, V = usv
+  U, S, V = usv
 
-    @test U * Diagonal(S) * V' ≈ a
-    @test U' * U ≈ LinearAlgebra.I
-    @test V' * V ≈ LinearAlgebra.I
+  @test U * Diagonal(S) * V' ≈ a
+  @test U' * U ≈ LinearAlgebra.I
+  @test V' * V ≈ LinearAlgebra.I
 end
 
 # regular matrix
@@ -18,34 +18,34 @@ end
 sizes = ((3, 3), (4, 3), (3, 4))
 eltypes = (Float32, Float64, ComplexF64)
 @testset "($m, $n) Matrix{$T}" for ((m, n), T) in Iterators.product(sizes, eltypes)
-    a = rand(m, n)
-    usv = @inferred svd(a)
-    test_svd(a, usv)
+  a = rand(m, n)
+  usv = @inferred svd(a)
+  test_svd(a, usv)
 end
 
 # block matrix
 # ------------
 blockszs = (([2, 2], [2, 2]), ([2, 2], [2, 3]), ([2, 2, 1], [2, 3]), ([2, 3], [2]))
 @testset "($m, $n) BlockMatrix{$T}" for ((m, n), T) in Iterators.product(blockszs, eltypes)
-    a = mortar([rand(T, i, j) for i in m, j in n])
-    usv = svd(a)
-    test_svd(a, usv)
-    @test usv.U isa BlockedMatrix
-    @test usv.Vt isa BlockedMatrix
-    @test usv.S isa BlockedVector
+  a = mortar([rand(T, i, j) for i in m, j in n])
+  usv = svd(a)
+  test_svd(a, usv)
+  @test usv.U isa BlockedMatrix
+  @test usv.Vt isa BlockedMatrix
+  @test usv.S isa BlockedVector
 end
 
 # Block-Diagonal matrices
 # -----------------------
 @testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in
                                           Iterators.product(blockszs, eltypes)
-    a = BlockDiagonal([rand(T, i, j) for (i, j) in zip(m, n)])
-    usv = svd(a)
-    # TODO: `BlockDiagonal * Adjoint` errors
-    test_svd(a, usv)
-    @test usv.U isa BlockDiagonal
-    @test usv.Vt isa BlockDiagonal
-    @test usv.S isa BlockVector
+  a = BlockDiagonal([rand(T, i, j) for (i, j) in zip(m, n)])
+  usv = svd(a)
+  # TODO: `BlockDiagonal * Adjoint` errors
+  test_svd(a, usv)
+  @test usv.U isa BlockDiagonal
+  @test usv.Vt isa BlockDiagonal
+  @test usv.S isa BlockVector
 end
 
 a = mortar([rand(2, 2) for i in 1:2, j in 1:3])
@@ -60,24 +60,24 @@ test_svd(a, usv)
 # -----------
 @testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in
                                           Iterators.product(blockszs, eltypes)
-    a = BlockSparseArray{T}(m, n)
-    for i in LinearAlgebra.diagind(blocks(a))
-        I = CartesianIndices(blocks(a))[i]
-        a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
-    end
-    perm = Random.randperm(length(m))
-    a = a[Block.(perm), Block.(1:length(n))]
+  a = BlockSparseArray{T}(m, n)
+  for i in LinearAlgebra.diagind(blocks(a))
+    I = CartesianIndices(blocks(a))[i]
+    a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
+  end
+  perm = Random.randperm(length(m))
+  a = a[Block.(perm), Block.(1:length(n))]
 
-    # errors because `blocks(a)[CartesianIndex.(...)]` is not implemented
-    usv = svd(a)
-    # TODO: `BlockDiagonal * Adjoint` errors
-    test_svd(a, usv)
-    @test usv.U isa BlockDiagonal
-    @test usv.Vt isa BlockDiagonal
-    @test usv.S isa BlockVector
+  # errors because `blocks(a)[CartesianIndex.(...)]` is not implemented
+  usv = svd(a)
+  # TODO: `BlockDiagonal * Adjoint` errors
+  test_svd(a, usv)
+  @test usv.U isa BlockDiagonal
+  @test usv.Vt isa BlockDiagonal
+  @test usv.S isa BlockVector
 
-    test_svd(a, usv2)
-    @test usv.U isa BlockDiagonal
-    @test usv.Vt isa BlockDiagonal
-    @test usv.S isa BlockVector
+  test_svd(a, usv2)
+  @test usv.U isa BlockDiagonal
+  @test usv.Vt isa BlockDiagonal
+  @test usv.S isa BlockVector
 end

From 1a1aef95d804641a6524b9a49988ebe7c227214d Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Thu, 9 Jan 2025 06:55:12 -0500
Subject: [PATCH 04/21] Rename variable + remove duplicate `Block` call

---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index cb5e374..109f0ac 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -7,9 +7,9 @@ function svd(
   A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A)
 )
   T = LinearAlgebra.eigtype(eltype(A))
-  A′ = try_to_blockdiagonal(A)
+  A′_and_blockperms = try_to_blockdiagonal(A)
 
-  if isnothing(A′)
+  if isnothing(A′_and_blockperms)
     # not block-diagonal, fall back to dense case
     Adense = eigencopy_oftype(A, T)
     return svd!(Adense; full, alg)
@@ -18,5 +18,5 @@ function svd(
   # compute block-by-block and permute back
   A″, (I, J) = A′
   F = svd!(eigencopy_oftype(A″, T); full, alg)
-  return SVD(F.U[Block.(I), Block.(J)], F.S, F.Vt)
+  return SVD(F.U[I, J], F.S, F.Vt)
 end

From 61977e4f643f0b82754320c907eaef524733b900 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Thu, 9 Jan 2025 12:45:12 -0500
Subject: [PATCH 05/21] Rewrite `svd` to work with rectangular matrices

---
 .../abstractblocksparsematrix.jl              | 80 ++++++++++++++++---
 src/blocksparsearray/blockdiagonalarray.jl    | 33 +-------
 2 files changed, 69 insertions(+), 44 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 109f0ac..f9305bf 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -3,20 +3,74 @@ const AbstractBlockSparseMatrix{T} = AbstractBlockSparseArray{T,2}
 # SVD is implemented by trying to
 # 1. Attempt to find a block-diagonal implementation by permuting
 # 2. Fallback to AbstractBlockArray implementation via BlockedArray
-function svd(
-  A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A)
-)
-  T = LinearAlgebra.eigtype(eltype(A))
-  A′_and_blockperms = try_to_blockdiagonal(A)
 
-  if isnothing(A′_and_blockperms)
-    # not block-diagonal, fall back to dense case
-    Adense = eigencopy_oftype(A, T)
-    return svd!(Adense; full, alg)
+function eigencopy_oftype(A::AbstractBlockSparseMatrix, T)
+  if is_block_permutation_matrix(A)
+    Acopy = similar(A, T)
+    for bI in eachblockstoredindex(A)
+      Acopy[bI] = eigencopy_oftype(A[bI], T)
+    end
+    return Acopy
+  else
+    return BlockedMatrix{T}(A)
+  end
+end
+
+function is_block_permutation_matrix(a::AbstractBlockSparseMatrix)
+  return allunique(first ∘ Tuple, eachblockstoredindex(a)) &&
+         allunique(last ∘ Tuple, eachblockstoredindex(a))
+end
+
+function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algorithm)
+  @assert !full "TODO"
+  bm, bn = blocksize(A)
+  bmn = min(bm, bn)
+
+  brows = blocklengths(axes(A, 1))
+  bcols = blocklengths(axes(A, 2))
+  slengths = Vector{Int}(undef, 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))
+    nrows = brows[row]
+    ncols = bcols[col]
+    slengths[col] = min(nrows, ncols)
   end
 
-  # compute block-by-block and permute back
-  A″, (I, J) = A′
-  F = svd!(eigencopy_oftype(A″, T); full, alg)
-  return SVD(F.U[I, J], F.S, F.Vt)
+  # 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 = findall(∉(browIs), 1:bmn)
+  emptycols = findall(∉(bcolIs), 1:bmn)
+  for (row, col) in zip(emptyrows, emptycols)
+    slengths[col] = min(brows[row], bcols[col])
+  end
+
+  U = similar(A, axes(A, 1), blockedrange(slengths))
+  S = similar(A, real(eltype(A)), blockedrange(slengths))
+  Vt = similar(A, blockedrange(slengths), axes(A, 2))
+
+  return U, S, Vt
+end
+
+function svd(A::AbstractBlockSparseMatrix; kwargs...)
+  return svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...)
+end
+
+function svd!(
+  A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A)
+)
+  @assert is_block_permutation_matrix(A) "Cannot keep sparsity: use `svd` to convert to `BlockedMatrix"
+  U, S, Vt = _allocate_svd_output(A, full, alg)
+  for bI in eachblockstoredindex(A)
+    bUSV = svd!(A[bI]; full, alg)
+    brow, bcol = Int.(Tuple(bI))
+    U[Block(brow, bcol)] = bUSV.U
+    S[Block(bcol)] = bUSV.S
+    Vt[Block(bcol, bcol)] = bUSV.Vt
+  end
+  return SVD(U, S, Vt)
 end
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index b900101..1aa0bac 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -11,38 +11,9 @@ function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
   )
 end
 
-# Cast to block-diagonal implementation if permuted-blockdiagonal
-function try_to_blockdiagonal_perm(A)
-  inds = map(x -> Int.(Tuple(x)), vec(collect(block_stored_indices(A))))
-  I = first.(inds)
-  allunique(I) || return nothing
-  J = last.(inds)
-  p = sortperm(J)
-  Jsorted = J[p]
-  allunique(Jsorted) || return nothing
-  return Block.(I[p], Jsorted)
-end
-
-"""
-    try_to_blockdiagonal(A)
-
-Attempt to find a permutation of blocks that makes `A` blockdiagonal. If unsuccesful,
-returns nothing, otherwise returns both the blockdiagonal `B` as well as the permutation `I, J`.
-"""
-function try_to_blockdiagonal(A::AbstractBlockSparseMatrix)
-  perm = try_to_blockdiagonal_perm(A)
-  isnothing(perm) && return perm
-  I = first.(Tuple.(perm))
-  J = last.(Tuple.(perm))
-  diagblocks = map(invperm(I), J) do i, j
-    return A[Block(i, j)]
-  end
-  return BlockDiagonal(diagblocks), perm
-end
-
 # SVD implementation
-function eigencopy_oftype(A::BlockDiagonal, S)
-  diag = map(Base.Fix2(eigencopy_oftype, S), A.blocks.diag)
+function eigencopy_oftype(A::BlockDiagonal, T)
+  diag = map(Base.Fix2(eigencopy_oftype, T), A.blocks.diag)
   return BlockDiagonal(diag)
 end
 

From a16ae7210d06149d515e61229a9edc9394f59891 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Thu, 9 Jan 2025 15:55:15 -0500
Subject: [PATCH 06/21] Update abstractblocksparsematrix.jl

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index f9305bf..933c814 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -67,10 +67,10 @@ function svd!(
   U, S, Vt = _allocate_svd_output(A, full, alg)
   for bI in eachblockstoredindex(A)
     bUSV = svd!(A[bI]; full, alg)
-    brow, bcol = Int.(Tuple(bI))
-    U[Block(brow, bcol)] = bUSV.U
-    S[Block(bcol)] = bUSV.S
-    Vt[Block(bcol, bcol)] = bUSV.Vt
+    brow, bcol = Tuple(bI)
+    U[brow, bcol] = bUSV.U
+    S[bcol] = bUSV.S
+    Vt[bcol, bcol] = bUSV.Vt
   end
   return SVD(U, S, Vt)
 end

From 3c061ecb3ac0c34cc041430931ff0f5a5b043a26 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Thu, 9 Jan 2025 12:46:45 -0500
Subject: [PATCH 07/21] update `BlockDiagonal` data access

---
 src/blocksparsearray/blockdiagonalarray.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 1aa0bac..120a62b 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -13,7 +13,7 @@ end
 
 # SVD implementation
 function eigencopy_oftype(A::BlockDiagonal, T)
-  diag = map(Base.Fix2(eigencopy_oftype, T), A.blocks.diag)
+  diag = map(Base.Fix2(eigencopy_oftype, T), LinearAlgebra.diag(blocks(A)))
   return BlockDiagonal(diag)
 end
 

From 3ecd6c8a8687787abfc45b2776d5707de15e864f Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Thu, 9 Jan 2025 15:56:12 -0500
Subject: [PATCH 08/21] Update abstractblocksparsematrix.jl

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 933c814..0b568bf 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -49,9 +49,10 @@ function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algori
     slengths[col] = min(brows[row], bcols[col])
   end
 
-  U = similar(A, axes(A, 1), blockedrange(slengths))
-  S = similar(A, real(eltype(A)), blockedrange(slengths))
-  Vt = similar(A, blockedrange(slengths), axes(A, 2))
+  s_axis = blockedrange(slengths)
+  U = similar(A, axes(A, 1), s_axis)
+  S = similar(A, real(eltype(A)), s_axis)
+  Vt = similar(A, s_axis, axes(A, 2))
 
   return U, S, Vt
 end

From 26390049068736f82e85ba03faa7c6179dadac29 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 08:49:21 -0500
Subject: [PATCH 09/21] remove `svd(::BlockDiagonal)` specialization

---
 src/blocksparsearray/blockdiagonalarray.jl | 23 +++++-----------------
 1 file changed, 5 insertions(+), 18 deletions(-)

diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 120a62b..63a4306 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -4,27 +4,14 @@ using LinearAlgebra: Diagonal
 const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{
   T,A,Diagonal{A,V},Axes
 }
+const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A}
+
+@interface interface::BlockSparseArrayInterface function blocks(a::BlockSparseDiagonal)
+  return Diagonal(Diagonal.(blocks(a.diag)))
+end
 
 function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
   return BlockSparseArray(
     Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
   )
 end
-
-# SVD implementation
-function eigencopy_oftype(A::BlockDiagonal, T)
-  diag = map(Base.Fix2(eigencopy_oftype, T), LinearAlgebra.diag(blocks(A)))
-  return BlockDiagonal(diag)
-end
-
-function svd(A::BlockDiagonal; kwargs...)
-  return svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...)
-end
-function svd!(A::BlockDiagonal; full::Bool=false, alg::Algorithm=default_svd_alg(A))
-  # TODO: handle full
-  F = map(a -> svd!(a; full, alg), blocks(A).diag)
-  Us = map(Base.Fix2(getproperty, :U), F)
-  Ss = map(Base.Fix2(getproperty, :S), F)
-  Vts = map(Base.Fix2(getproperty, :Vt), F)
-  return SVD(BlockDiagonal(Us), mortar(Ss), BlockDiagonal(Vts))
-end

From 054ddf60937ed832d4777186ef5376b32b13c4d6 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 08:50:14 -0500
Subject: [PATCH 10/21] small update to tests

---
 test/basics/test_svd.jl | 14 +++-----------
 1 file changed, 3 insertions(+), 11 deletions(-)

diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl
index 00903d4..934028e 100644
--- a/test/basics/test_svd.jl
+++ b/test/basics/test_svd.jl
@@ -1,6 +1,6 @@
 using Test
 using BlockSparseArrays
-using BlockSparseArrays: BlockSparseArray, svd, notrunc, truncbelow, truncdim, BlockDiagonal
+using BlockSparseArrays: BlockSparseArray, svd, BlockDiagonal
 using BlockArrays
 using LinearAlgebra: LinearAlgebra, Diagonal, svdvals
 using Random
@@ -58,8 +58,8 @@ test_svd(a, usv)
 
 # blocksparse 
 # -----------
-@testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in
-                                          Iterators.product(blockszs, eltypes)
+@testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in
+                                              Iterators.product(blockszs, eltypes)
   a = BlockSparseArray{T}(m, n)
   for i in LinearAlgebra.diagind(blocks(a))
     I = CartesianIndices(blocks(a))[i]
@@ -72,12 +72,4 @@ test_svd(a, usv)
   usv = svd(a)
   # TODO: `BlockDiagonal * Adjoint` errors
   test_svd(a, usv)
-  @test usv.U isa BlockDiagonal
-  @test usv.Vt isa BlockDiagonal
-  @test usv.S isa BlockVector
-
-  test_svd(a, usv2)
-  @test usv.U isa BlockDiagonal
-  @test usv.Vt isa BlockDiagonal
-  @test usv.S isa BlockVector
 end

From 9073fbf264c72d86223ea5cdc83d58972c9df895 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 08:52:13 -0500
Subject: [PATCH 11/21] avoid unnecesary copy

---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 0b568bf..b746b64 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -67,7 +67,7 @@ function svd!(
   @assert is_block_permutation_matrix(A) "Cannot keep sparsity: use `svd` to convert to `BlockedMatrix"
   U, S, Vt = _allocate_svd_output(A, full, alg)
   for bI in eachblockstoredindex(A)
-    bUSV = svd!(A[bI]; full, alg)
+    bUSV = svd!(@view!(A[bI]); full, alg)
     brow, bcol = Tuple(bI)
     U[brow, bcol] = bUSV.U
     S[bcol] = bUSV.S

From 2e97b97dc9083a8e174d1d2d5c4fc688a3addb3d Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 11:59:49 -0500
Subject: [PATCH 12/21] Fix docs

---
 docs/make.jl              | 10 ++++++++--
 docs/src/library.md       |  5 +++++
 src/factorizations/svd.jl | 24 ++++++++++++------------
 3 files changed, 25 insertions(+), 14 deletions(-)
 create mode 100644 docs/src/library.md

diff --git a/docs/make.jl b/docs/make.jl
index d24e6c3..18b0932 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,7 +2,13 @@ using BlockSparseArrays: BlockSparseArrays
 using Documenter: Documenter, DocMeta, deploydocs, makedocs
 
 DocMeta.setdocmeta!(
-  BlockSparseArrays, :DocTestSetup, :(using BlockSparseArrays); recursive=true
+  BlockSparseArrays,
+  :DocTestSetup,
+  quote
+    using BlockSparseArrays
+    using LinearAlgebra: Diagonal
+  end;
+  recursive=true,
 )
 
 include("make_index.jl")
@@ -16,7 +22,7 @@ makedocs(;
     edit_link="main",
     assets=String[],
   ),
-  pages=["Home" => "index.md"],
+  pages=["Home" => "index.md", "Library" => "library.md"],
 )
 
 deploydocs(;
diff --git a/docs/src/library.md b/docs/src/library.md
new file mode 100644
index 0000000..63af0e8
--- /dev/null
+++ b/docs/src/library.md
@@ -0,0 +1,5 @@
+# Library
+
+```@autodocs
+Modules = [BlockSparseArrays]
+```
diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 9e46dc6..00d7287 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -33,14 +33,14 @@ julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
  0.0  0.0  0.0  0.0  0.0
  0.0  2.0  0.0  0.0  0.0
 
-julia> F = svd(A)
-SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
+julia> F = BlockSparseArrays.svd(A)
+BlockSparseArrays.SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}, Matrix{Float64}}
 U factor:
 4×4 Matrix{Float64}:
- 0.0  1.0   0.0  0.0
- 1.0  0.0   0.0  0.0
- 0.0  0.0   0.0  1.0
- 0.0  0.0  -1.0  0.0
+ 0.0  1.0  0.0   0.0
+ 1.0  0.0  0.0   0.0
+ 0.0  0.0  0.0  -1.0
+ 0.0  0.0  1.0   0.0
 singular values:
 4-element Vector{Float64}:
  3.0
@@ -49,10 +49,10 @@ singular values:
  0.0
 Vt factor:
 4×5 Matrix{Float64}:
- -0.0        0.0  1.0  -0.0  0.0
-  0.447214   0.0  0.0   0.0  0.894427
-  0.0       -1.0  0.0   0.0  0.0
-  0.0        0.0  0.0   1.0  0.0
+ -0.0       0.0  1.0  -0.0  0.0
+  0.447214  0.0  0.0   0.0  0.894427
+ -0.0       1.0  0.0  -0.0  0.0
+  0.0       0.0  0.0   1.0  0.0
 
 julia> F.U * Diagonal(F.S) * F.Vt
 4×5 Matrix{Float64}:
@@ -181,7 +181,7 @@ number of singular values.
 ```jldoctest
 julia> A = rand(4,3);
 
-julia> F = svd(A); # Store the Factorization Object
+julia> F = BlockSparseArrays.svd(A); # Store the Factorization Object
 
 julia> A ≈ F.U * Diagonal(F.S) * F.Vt
 true
@@ -191,7 +191,7 @@ julia> U, S, V = F; # destructuring via iteration
 julia> A ≈ U * Diagonal(S) * V'
 true
 
-julia> Uonly, = svd(A); # Store U only
+julia> Uonly, = BlockSparseArrays.svd(A); # Store U only
 
 julia> Uonly == U
 true

From 21e6d5d8eab6b5a09372be159fbdb28de222e7d8 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 12:19:00 -0500
Subject: [PATCH 13/21] Apply code suggestion

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/factorizations/svd.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index 00d7287..bdd9a66 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -200,7 +200,7 @@ true
 svd(A; kwargs...) =
   SVD(svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...))
 
-LinearAlgebra.svdvals(usv::SVD{<:Any,T}) where {T} = (usv.S)::Vector{T}
+LinearAlgebra.svdvals(usv::SVD{<:Any,T}) where {T} = (usv.S)::AbstractVector{T}
 
 # Added here to avoid type-piracy
 eigencopy_oftype(A, S) = LinearAlgebra.eigencopy_oftype(A, S)

From 73e5c58d2a66632753b53ec4d779605fef9c0381 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Fri, 10 Jan 2025 12:22:38 -0500
Subject: [PATCH 14/21] Apply code suggestion

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index b746b64..6fd57d7 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -8,7 +8,7 @@ function eigencopy_oftype(A::AbstractBlockSparseMatrix, T)
   if is_block_permutation_matrix(A)
     Acopy = similar(A, T)
     for bI in eachblockstoredindex(A)
-      Acopy[bI] = eigencopy_oftype(A[bI], T)
+      Acopy[bI] = eigencopy_oftype(@view!(A[bI]), T)
     end
     return Acopy
   else

From 9169c5264d2602d77e574261fd6bc5e5e2ee7624 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 13 Jan 2025 16:25:46 -0500
Subject: [PATCH 15/21] catch non-full sparse matrices and simplify diagonal
 implementation

---
 Project.toml                                  |  2 +
 .../abstractblocksparsematrix.jl              | 19 ++++++++-
 src/blocksparsearray/blockdiagonalarray.jl    |  9 ++++
 test/Project.toml                             |  1 +
 test/basics/test_svd.jl                       | 41 ++++++++++---------
 5 files changed, 51 insertions(+), 21 deletions(-)

diff --git a/Project.toml b/Project.toml
index 056f7f3..d6a8a57 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DerivableInterfaces = "6c5e35bf-e59e-4898-b73c-732dcc4ba65f"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
@@ -32,6 +33,7 @@ Aqua = "0.8.9"
 ArrayLayouts = "1.10.4"
 BlockArrays = "1.2.0"
 DerivableInterfaces = "0.3.7"
+DiagonalArrays = "0.2.2"
 Dictionaries = "0.4.3"
 FillArrays = "1.13.0"
 GPUArraysCore = "0.1.0, 0.2"
diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 6fd57d7..0481355 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -43,8 +43,8 @@ function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algori
 
   # 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 = findall(∉(browIs), 1:bmn)
-  emptycols = findall(∉(bcolIs), 1:bmn)
+  emptyrows = findall(∉(browIs), 1:bm)
+  emptycols = findall(∉(bcolIs), 1:bn)
   for (row, col) in zip(emptyrows, emptycols)
     slengths[col] = min(brows[row], bcols[col])
   end
@@ -73,5 +73,20 @@ function svd!(
     S[bcol] = bUSV.S
     Vt[bcol, bcol] = bUSV.Vt
   end
+
+  # fill in values for blocks that aren't present, pairing them in order of occurence
+  block_inds_S = eachblockstoredindex(S)
+  i = findfirst(∉(block_inds_S), blockaxes(S, 1))
+  bIs = collect(eachblockstoredindex(U))
+  browIs = Int.(first.(Tuple.(bIs)))
+  emptyrows = findall(∉(browIs), 1:blocksize(U, 1))
+  j = 0
+  while !isnothing(i)
+    copyto!(@view!(Vt[Block(i), Block(i)]), LinearAlgebra.I)
+    j += 1
+    copyto!(@view!(U[Block(emptyrows[j]), Block(i)]), LinearAlgebra.I)
+    i = findnext(∉(block_inds_S), blockaxes(S, 1), i + 1)
+  end
+
   return SVD(U, S, Vt)
 end
diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl
index 63a4306..7d0f281 100644
--- a/src/blocksparsearray/blockdiagonalarray.jl
+++ b/src/blocksparsearray/blockdiagonalarray.jl
@@ -1,5 +1,6 @@
 # type alias for block-diagonal
 using LinearAlgebra: Diagonal
+using DiagonalArrays: DiagonalArrays, diagonal
 
 const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{
   T,A,Diagonal{A,V},Axes
@@ -15,3 +16,11 @@ function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix})
     Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2)))
   )
 end
+
+function DiagonalArrays.diagonal(S::BlockSparseVector)
+  D = similar(S, (axes(S, 1), axes(S, 1)))
+  for bI in eachblockstoredindex(S)
+    D[bI, bI] = diagonal(@view!(S[bI]))
+  end
+  return D
+end
diff --git a/test/Project.toml b/test/Project.toml
index 4268ad3..0cc3624 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl
index 934028e..70dc984 100644
--- a/test/basics/test_svd.jl
+++ b/test/basics/test_svd.jl
@@ -1,14 +1,15 @@
 using Test
 using BlockSparseArrays
-using BlockSparseArrays: BlockSparseArray, svd, BlockDiagonal
+using BlockSparseArrays: BlockSparseArray, svd, BlockDiagonal, eachblockstoredindex
 using BlockArrays
-using LinearAlgebra: LinearAlgebra, Diagonal, svdvals
 using Random
+using DiagonalArrays: diagonal
+using LinearAlgebra: LinearAlgebra
 
 function test_svd(a, usv)
   U, S, V = usv
 
-  @test U * Diagonal(S) * V' ≈ a
+  @test U * diagonal(S) * V' ≈ a
   @test U' * U ≈ LinearAlgebra.I
   @test V' * V ≈ LinearAlgebra.I
 end
@@ -43,33 +44,35 @@ end
   usv = svd(a)
   # TODO: `BlockDiagonal * Adjoint` errors
   test_svd(a, usv)
-  @test usv.U isa BlockDiagonal
-  @test usv.Vt isa BlockDiagonal
-  @test usv.S isa BlockVector
 end
 
-a = mortar([rand(2, 2) for i in 1:2, j in 1:3])
-usv = svd(a)
-test_svd(a, usv)
-
-a = mortar([rand(2, 2) for i in 1:3, j in 1:2])
-usv = svd(a)
-test_svd(a, usv)
-
 # blocksparse 
 # -----------
 @testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in
                                               Iterators.product(blockszs, eltypes)
   a = BlockSparseArray{T}(m, n)
+
+  # test empty matrix
+  usv_empty = svd(a)
+  test_svd(a, usv_empty)
+
+  # test blockdiagonal
   for i in LinearAlgebra.diagind(blocks(a))
     I = CartesianIndices(blocks(a))[i]
     a[Block(I.I...)] = rand(T, size(blocks(a)[i]))
   end
-  perm = Random.randperm(length(m))
-  a = a[Block.(perm), Block.(1:length(n))]
-
-  # errors because `blocks(a)[CartesianIndex.(...)]` is not implemented
   usv = svd(a)
-  # TODO: `BlockDiagonal * Adjoint` errors
   test_svd(a, usv)
+
+  perm = Random.randperm(length(m))
+  b = a[Block.(perm), Block.(1:length(n))]
+  usv = svd(b)
+  test_svd(b, usv)
+
+  # test permuted blockdiagonal with missing row/col
+  I_removed = rand(eachblockstoredindex(b))
+  c = copy(b)
+  delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed))))
+  usv = svd(c)
+  test_svd(c, usv)
 end

From f03bed2c5fb92e3638119d3920fb9ab6dee80f92 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 13 Jan 2025 16:33:09 -0500
Subject: [PATCH 16/21] Bump minimal version `DerivableInterfaces`

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d6a8a57..1a36e22 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ Adapt = "4.1.1"
 Aqua = "0.8.9"
 ArrayLayouts = "1.10.4"
 BlockArrays = "1.2.0"
-DerivableInterfaces = "0.3.7"
+DerivableInterfaces = "0.3.8"
 DiagonalArrays = "0.2.2"
 Dictionaries = "0.4.3"
 FillArrays = "1.13.0"

From d9b307ff508c9aacb6a08857832c1dbec48cc618 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 13 Jan 2025 16:38:57 -0500
Subject: [PATCH 17/21] Fix docstring tests (CI runs on different LAPACK)

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

diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl
index bdd9a66..3c82f23 100644
--- a/src/factorizations/svd.jl
+++ b/src/factorizations/svd.jl
@@ -37,10 +37,10 @@ julia> F = BlockSparseArrays.svd(A)
 BlockSparseArrays.SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}, Matrix{Float64}}
 U factor:
 4×4 Matrix{Float64}:
- 0.0  1.0  0.0   0.0
- 1.0  0.0  0.0   0.0
- 0.0  0.0  0.0  -1.0
- 0.0  0.0  1.0   0.0
+ 0.0  1.0   0.0  0.0
+ 1.0  0.0   0.0  0.0
+ 0.0  0.0   0.0  1.0
+ 0.0  0.0  -1.0  0.0
 singular values:
 4-element Vector{Float64}:
  3.0
@@ -49,10 +49,10 @@ singular values:
  0.0
 Vt factor:
 4×5 Matrix{Float64}:
- -0.0       0.0  1.0  -0.0  0.0
-  0.447214  0.0  0.0   0.0  0.894427
- -0.0       1.0  0.0  -0.0  0.0
-  0.0       0.0  0.0   1.0  0.0
+ -0.0        0.0  1.0  -0.0  0.0
+  0.447214   0.0  0.0   0.0  0.894427
+  0.0       -1.0  0.0   0.0  0.0
+  0.0        0.0  0.0   1.0  0.0
 
 julia> F.U * Diagonal(F.S) * F.Vt
 4×5 Matrix{Float64}:

From 857f8c04a8201ba5d3701f50f676b13b9ec1cad3 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Jan 2025 12:54:06 -0500
Subject: [PATCH 18/21] Apply suggestions from code review

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 0481355..b87d774 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -82,9 +82,9 @@ function svd!(
   emptyrows = findall(∉(browIs), 1:blocksize(U, 1))
   j = 0
   while !isnothing(i)
-    copyto!(@view!(Vt[Block(i), Block(i)]), LinearAlgebra.I)
+    copyto!(@view!(Vt[Block(i, i)]), LinearAlgebra.I)
     j += 1
-    copyto!(@view!(U[Block(emptyrows[j]), Block(i)]), LinearAlgebra.I)
+    copyto!(@view!(U[Block(emptyrows[j], i)]), LinearAlgebra.I)
     i = findnext(∉(block_inds_S), blockaxes(S, 1), i + 1)
   end
 

From da3296b238aa9b079865584fb0572e4b0ce9a1b6 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Jan 2025 13:09:09 -0500
Subject: [PATCH 19/21] reorganize allocation of identities

---
 .../abstractblocksparsematrix.jl              | 20 ++++++-------------
 1 file changed, 6 insertions(+), 14 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index b87d774..1a1ebce 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -54,6 +54,12 @@ function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algori
   S = similar(A, real(eltype(A)), s_axis)
   Vt = similar(A, s_axis, axes(A, 2))
 
+  # also fill in identities for blocks that aren't present
+  for (row, col) in zip(emptyrows, emptycols)
+    copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I)
+    copyto!(@view!(Vt[Block(col, col)]), LinearAlgebra.I)
+  end
+
   return U, S, Vt
 end
 
@@ -74,19 +80,5 @@ function svd!(
     Vt[bcol, bcol] = bUSV.Vt
   end
 
-  # fill in values for blocks that aren't present, pairing them in order of occurence
-  block_inds_S = eachblockstoredindex(S)
-  i = findfirst(∉(block_inds_S), blockaxes(S, 1))
-  bIs = collect(eachblockstoredindex(U))
-  browIs = Int.(first.(Tuple.(bIs)))
-  emptyrows = findall(∉(browIs), 1:blocksize(U, 1))
-  j = 0
-  while !isnothing(i)
-    copyto!(@view!(Vt[Block(i, i)]), LinearAlgebra.I)
-    j += 1
-    copyto!(@view!(U[Block(emptyrows[j], i)]), LinearAlgebra.I)
-    i = findnext(∉(block_inds_S), blockaxes(S, 1), i + 1)
-  end
-
   return SVD(U, S, Vt)
 end

From f54c269f5f7b6ce1c7853eafd83eaa677514b18d Mon Sep 17 00:00:00 2001
From: Matt Fishman <mtfishman@users.noreply.github.com>
Date: Tue, 14 Jan 2025 14:03:51 -0500
Subject: [PATCH 20/21] Bump to v0.2.7

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 1a36e22..2a2d8b2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.2.6"
+version = "0.2.7"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From af4a5954501627b2d23144d196c23dcbbddafefe Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Jan 2025 14:15:19 -0500
Subject: [PATCH 21/21] Refactor `findall` in terms of `setdiff`

Co-authored-by: Matt Fishman <mtfishman@users.noreply.github.com>
---
 src/abstractblocksparsearray/abstractblocksparsematrix.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
index 1a1ebce..376f408 100644
--- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl
+++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl
@@ -43,8 +43,8 @@ function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algori
 
   # 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 = findall(∉(browIs), 1:bm)
-  emptycols = findall(∉(bcolIs), 1:bn)
+  emptyrows = setdiff(1:bm, browIs)
+  emptycols = setdiff(1:bn, bcolIs)
   for (row, col) in zip(emptyrows, emptycols)
     slengths[col] = min(brows[row], bcols[col])
   end