From da3296b238aa9b079865584fb0572e4b0ce9a1b6 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 14 Jan 2025 13:09:09 -0500 Subject: [PATCH] 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