From 94c50149db617fd9a60dcae7c5a5ce0fa8fcb4e0 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 28 Apr 2025 16:12:47 -0400 Subject: [PATCH 1/3] Add MatrixAlgebraKit factorizations --- Project.toml | 2 + src/BlockSparseArrays.jl | 5 +- .../getunstoredblock.jl | 15 + src/factorizations/svd.jl | 394 ++++++++++-------- test/Project.toml | 1 + test/test_factorizations.jl | 85 ++++ test/test_svd.jl | 75 ---- 7 files changed, 319 insertions(+), 258 deletions(-) create mode 100644 test/test_factorizations.jl delete mode 100644 test/test_svd.jl diff --git a/Project.toml b/Project.toml index 1fdf66a..3d99835 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" MapBroadcast = "ebd9b9da-f48d-417c-9660-449667d60261" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" @@ -38,6 +39,7 @@ GPUArraysCore = "0.1.0, 0.2" LinearAlgebra = "1.10" MacroTools = "0.5.13" MapBroadcast = "0.1.5" +MatrixAlgebraKit = "0.1.2" SparseArraysBase = "0.5" SplitApplyCombine = "1.2.3" TensorAlgebra = "0.3.2" diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl index ff2ed6e..f261233 100644 --- a/src/BlockSparseArrays.jl +++ b/src/BlockSparseArrays.jl @@ -8,8 +8,6 @@ export BlockSparseArray, eachstoredblock, sparsemortar -# factorizations -include("factorizations/svd.jl") # possible upstream contributions include("BlockArraysExtensions/blockedunitrange.jl") @@ -45,4 +43,7 @@ include("blocksparsearray/blockdiagonalarray.jl") include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl") +# factorizations +include("factorizations/svd.jl") + end diff --git a/src/blocksparsearrayinterface/getunstoredblock.jl b/src/blocksparsearrayinterface/getunstoredblock.jl index 43b58c0..b8c2425 100644 --- a/src/blocksparsearrayinterface/getunstoredblock.jl +++ b/src/blocksparsearrayinterface/getunstoredblock.jl @@ -23,3 +23,18 @@ end ) where {N} return f(a, Tuple(I)...) end + +# TODO: this is a hack and is also type-unstable +function (f::GetUnstoredBlock)( + a::AbstractMatrix{LinearAlgebra.Diagonal{T,V}}, I::Vararg{Int,2} +) where {T,V} + b_size = ntuple(ndims(a)) do d + return length(f.axes[d][Block(I[d])]) + end + if I[1] == I[2] + diag = similar(V, b_size[1]) + return LinearAlgebra.Diagonal{T,V}(diag) + else + return zeros(T, b_size...) + end +end diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index f25a50c..1826a47 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -1,207 +1,239 @@ -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. +using MatrixAlgebraKit + """ - 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 = 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 -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 -``` + BlockDiagonalAlgorithm(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. """ -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), - ) +struct BlockDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: + MatrixAlgebraKit.AbstractAlgorithm + alg::A 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), - ) +function MatrixAlgebraKit.default_svd_algorithm(A::AbstractBlockSparseMatrix; kwargs...) + @assert blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} "unsupported type: + $(blocktype(A))" + alg = MatrixAlgebraKit.LAPACK_DivideAndConquer(; kwargs...) + return BlockDiagonalAlgorithm(alg) 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) +#= +Note: here I'm being generic about the matrix type, which effectively means that I'm making +some assumptions about the output type of the algorithm, ie that this will return +Matrix{T},Diagonal{real(T)},Matrix{T}. In principle this is not guaranteed, although it +should cover most cases. The alternative is to be more specific about the matrix type and +replace the `similar` calls with explicit `BlockSparseArray` constructor calls. In that case +I can simply call `initialize_output` on the input blocks, and take whatever is returned. +We should probably discuss which way to go. +=# + +function MatrixAlgebraKit.initialize_output( + ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm +) + 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 -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)) + # 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) + slengths[col] = min(brows[row], bcols[col]) + end -# 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) + s_axis = blockedrange(slengths) + U = similar(A, axes(A, 1), s_axis) + Tr = real(eltype(A)) + S = BlockSparseArray{Tr,2,Diagonal{Tr,Vector{Tr}}}(undef, (s_axis, s_axis)) + Vt = similar(A, s_axis, axes(A, 2)) + + # allocate output + for bI in eachblockstoredindex(A) + brow, bcol = Tuple(bI) + U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output( + svd_compact!, @view!(A[bI]), alg.alg + ) + end -# functions default to LinearAlgebra -# ---------------------------------- -""" - svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD + # 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 -`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...)) + return U, S, Vt +end -""" - svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD +function MatrixAlgebraKit.initialize_output( + ::typeof(svd_full!), A::AbstractBlockSparseMatrix, alg::BlockDiagonalAlgorithm +) + bm, bn = blocksize(A) + + brows = blocklengths(axes(A, 1)) + slengths = 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)) + nrows = brows[row] + slengths[col] = nrows + end -Compute the singular value decomposition (SVD) of `A` and return an `SVD` object. + # 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) + slengths[col] = brows[row] + end + for (i, k) in enumerate((length(emptycols) + 1):length(emptyrows)) + slengths[bn + i] = brows[emptyrows[k]] + end -`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. + s_axis = blockedrange(slengths) + U = similar(A, axes(A, 1), s_axis) + Tr = real(eltype(A)) + S = similar(A, Tr, (s_axis, axes(A, 2))) + Vt = similar(A, axes(A, 2), axes(A, 2)) + + # allocate output + for bI in eachblockstoredindex(A) + brow, bcol = Tuple(bI) + U[brow, bcol], S[bcol, bcol], Vt[bcol, bcol] = MatrixAlgebraKit.initialize_output( + svd_full!, @view!(A[bI]), alg.alg + ) + end -Iterating the decomposition produces the components `U`, `S`, and `V`. + # 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 -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. + return U, S, Vt +end -`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) . +function MatrixAlgebraKit.check_input( + ::typeof(svd_compact!), A::AbstractBlockSparseMatrix, USVᴴ +) + U, S, Vt = USVᴴ + @assert isa(U, AbstractBlockSparseMatrix) && + isa(S, AbstractBlockSparseMatrix) && + isa(Vt, AbstractBlockSparseMatrix) + @assert eltype(A) == eltype(U) == eltype(Vt) + @assert real(eltype(A)) == eltype(S) + @assert axes(A, 1) == axes(U, 1) && axes(A, 2) == axes(Vt, 2) + @assert axes(S, 1) == axes(S, 2) + + # TODO: implement checks on axes of S, or find better way to do this without recomputing + # pairing + return nothing +end -!!! compat "Julia 1.3" - The `alg` keyword argument requires Julia 1.3 or later. +function MatrixAlgebraKit.check_input( + ::typeof(svd_full!), A::AbstractBlockSparseMatrix, USVᴴ +) + U, S, Vt = USVᴴ + @assert isa(U, AbstractBlockSparseMatrix) && + isa(S, AbstractBlockSparseMatrix) && + isa(Vt, AbstractBlockSparseMatrix) + @assert eltype(A) == eltype(U) == eltype(Vt) + @assert real(eltype(A)) == eltype(S) + @assert axes(A, 1) == axes(U, 1) && axes(A, 2) == axes(Vt, 1) == axes(Vt, 2) + @assert axes(S, 2) == axes(A, 2) + + # TODO: implement checks on axes of S, or find better way to do this without recomputing + # pairing + return nothing +end -# Examples -```jldoctest -julia> A = rand(4,3); +function MatrixAlgebraKit.svd_compact!( + A::AbstractBlockSparseMatrix, USVᴴ, alg::BlockDiagonalAlgorithm +) + MatrixAlgebraKit.check_input(svd_compact!, A, USVᴴ) + U, S, Vt = USVᴴ + + # do decomposition on each block + for bI in eachblockstoredindex(A) + brow, bcol = Tuple(bI) + usvᴴ = (@view!(U[brow, bcol]), @view!(S[bcol, bcol]), @view!(Vt[bcol, bcol])) + usvᴴ′ = svd_compact!(@view!(A[bI]), usvᴴ, alg.alg) + @assert usvᴴ === usvᴴ′ "svd_compact! might not be in-place" + end -julia> F = BlockSparseArrays.svd(A); # Store the Factorization Object + # 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) + for (row, col) in zip(emptyrows, emptycols) + copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I) + copyto!(@view!(Vt[Block(col, col)]), LinearAlgebra.I) + end -julia> A ≈ F.U * Diagonal(F.S) * F.Vt -true + return USVᴴ +end -julia> U, S, V = F; # destructuring via iteration +function MatrixAlgebraKit.svd_full!( + A::AbstractBlockSparseMatrix, USVᴴ, alg::BlockDiagonalAlgorithm +) + MatrixAlgebraKit.check_input(svd_full!, A, USVᴴ) + U, S, Vt = USVᴴ + + # do decomposition on each block + for bI in eachblockstoredindex(A) + brow, bcol = Tuple(bI) + usvᴴ = (@view!(U[brow, bcol]), @view!(S[bcol, bcol]), @view!(Vt[bcol, bcol])) + usvᴴ′ = svd_full!(@view!(A[bI]), usvᴴ, alg.alg) + @assert usvᴴ === usvᴴ′ "svd_full! might not be in-place" + end -julia> A ≈ U * Diagonal(S) * V' -true + # 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) + for (row, col) in zip(emptyrows, emptycols) + copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I) + copyto!(@view!(Vt[Block(col, col)]), LinearAlgebra.I) + end -julia> Uonly, = BlockSparseArrays.svd(A); # Store U only + # also handle extra rows/cols + for i in (length(emptyrows) + 1):length(emptycols) + copyto!(@view!(Vt[Block(emptycols[i], emptycols[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) + end -julia> Uonly == U -true -``` -""" -function svd(A; kwargs...) - return SVD(svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...)) + return USVᴴ end - -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) diff --git a/test/Project.toml b/test/Project.toml index 455fe69..5516f04 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,6 +8,7 @@ DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl new file mode 100644 index 0000000..7a28fdb --- /dev/null +++ b/test/test_factorizations.jl @@ -0,0 +1,85 @@ +using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar +using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex +using MatrixAlgebraKit: svd_compact, svd_full +using LinearAlgebra: LinearAlgebra +using Random: Random +using Test: @inferred, @testset, @test + +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 + + # Check factors are unitary + (U * U' ≈ LinearAlgebra.I) || return false + (Vᴴ' * Vᴴ ≈ LinearAlgebra.I) || return false + return true +end + +blockszs = ( + ([2, 2], [2, 2]), ([2, 2], [2, 3]), ([2, 2, 1], [2, 3]), ([2, 3], [2]), ([2], [2, 3]) +) +eltypes = (Float32, Float64, ComplexF64) +test_params = Iterators.product(blockszs, eltypes) + +# svd_compact! +# ------------ +@testset "svd_compact ($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in test_params + a = BlockSparseArray{T}(undef, m, n) + + # test empty matrix + usv_empty = svd_compact(a) + @test 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 + usv = svd_compact(a) + @test test_svd(a, usv) + + perm = Random.randperm(length(m)) + b = a[Block.(perm), Block.(1:length(n))] + usv = svd_compact(b) + @test 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_compact(c) + @test test_svd(c, usv) +end + +# svd_full! +# --------- +@testset "svd_full ($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in test_params + a = BlockSparseArray{T}(undef, m, n) + + # test empty matrix + usv_empty = svd_full(a) + @test test_svd(a, usv_empty; full=true) + + # 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 + usv = svd_full(a) + @test test_svd(a, usv; full=true) + + perm = Random.randperm(length(m)) + b = a[Block.(perm), Block.(1:length(n))] + usv = svd_full(b) + @test test_svd(b, usv; full=true) + + # 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_full(c) + @test test_svd(c, usv; full=true) +end diff --git a/test/test_svd.jl b/test/test_svd.jl deleted file mode 100644 index 73efa16..0000000 --- a/test/test_svd.jl +++ /dev/null @@ -1,75 +0,0 @@ -using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar -using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex, svd -using DiagonalArrays: diagonal -using LinearAlgebra: LinearAlgebra -using Random: Random -using Test: @inferred, @testset, @test - -function test_svd(a, usv) - U, S, V = usv - return (U * diagonal(S) * V' ≈ a) && - (U' * U ≈ LinearAlgebra.I) && - (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 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 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) - @test test_svd(a, usv) -end - -# blocksparse -# ----------- -@testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in - Iterators.product(blockszs, eltypes) - a = BlockSparseArray{T}(undef, m, n) - - # test empty matrix - usv_empty = svd(a) - @test 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 - usv = svd(a) - @test test_svd(a, usv) - - perm = Random.randperm(length(m)) - b = a[Block.(perm), Block.(1:length(n))] - usv = svd(b) - @test 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 test_svd(c, usv) -end From 8b9f91d89dd40d22c9d8c06cdd9731c23633f564 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 28 Apr 2025 16:30:01 -0400 Subject: [PATCH 2/3] formatter --- src/BlockSparseArrays.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl index f261233..78ce823 100644 --- a/src/BlockSparseArrays.jl +++ b/src/BlockSparseArrays.jl @@ -8,7 +8,6 @@ export BlockSparseArray, eachstoredblock, sparsemortar - # possible upstream contributions include("BlockArraysExtensions/blockedunitrange.jl") include("BlockArraysExtensions/BlockArraysExtensions.jl") From 7840264b694011f3d7693949c9849569a1e5cc35 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 28 Apr 2025 16:33:15 -0400 Subject: [PATCH 3/3] Bump version --- Project.toml | 2 +- docs/Project.toml | 2 +- examples/Project.toml | 2 +- test/Project.toml | 3 ++- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 3d99835..b13be4a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.4.6" +version = "0.5.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/Project.toml b/docs/Project.toml index a5c0d2e..2bf7adb 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" [compat] BlockArrays = "1" -BlockSparseArrays = "0.4" +BlockSparseArrays = "0.5" Documenter = "1" Literate = "2" diff --git a/examples/Project.toml b/examples/Project.toml index a847c67..08e9378 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] BlockArrays = "1" -BlockSparseArrays = "0.4" +BlockSparseArrays = "0.5" Test = "1" diff --git a/test/Project.toml b/test/Project.toml index 5516f04..7d8eb2e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,11 +23,12 @@ Adapt = "4" Aqua = "0.8" ArrayLayouts = "1" BlockArrays = "1" -BlockSparseArrays = "0.4" +BlockSparseArrays = "0.5" DiagonalArrays = "0.3" GPUArraysCore = "0.2" JLArrays = "0.2" LinearAlgebra = "1" +MatrixAlgebraKit = "0.1" Random = "1" SafeTestsets = "0.1" SparseArraysBase = "0.5"