diff --git a/CHANGELOG.md b/CHANGELOG.md index a9938d93a..d61a642fb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main) +- Change the structure of block diagonalization functions, using `BlockDiagonalForm` struct and changing the function name from `bdf` to `block_diagonal_form`. ([#349]) + ## [v0.24.0] Release date: 2024-12-13 @@ -50,6 +52,7 @@ Release date: 2024-11-13 [v0.22.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.22.0 [v0.23.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.0 [v0.23.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.23.1 +[v0.24.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.24.0 [#86]: https://github.com/qutip/QuantumToolbox.jl/issues/86 [#139]: https://github.com/qutip/QuantumToolbox.jl/issues/139 [#292]: https://github.com/qutip/QuantumToolbox.jl/issues/292 @@ -66,3 +69,4 @@ Release date: 2024-11-13 [#342]: https://github.com/qutip/QuantumToolbox.jl/issues/342 [#346]: https://github.com/qutip/QuantumToolbox.jl/issues/346 [#347]: https://github.com/qutip/QuantumToolbox.jl/issues/347 +[#349]: https://github.com/qutip/QuantumToolbox.jl/issues/349 diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md index 900fed82d..b701d8e7f 100644 --- a/docs/src/resources/api.md +++ b/docs/src/resources/api.md @@ -258,6 +258,13 @@ MultiSiteOperator DissipativeIsing ``` +## [Symmetries and Block Diagonalization](@id doc-API:Symmetries-and-Block-Diagonalization) + +```@docs +block_diagonal_form +BlockDiagonalForm +``` + ## [Miscellaneous](@id doc-API:Miscellaneous) ```@docs diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index ad387e261..28a50f251 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -90,6 +90,7 @@ include("qobj/states.jl") include("qobj/operators.jl") include("qobj/superoperators.jl") include("qobj/synonyms.jl") +include("qobj/block_diagonal_form.jl") # time evolution include("time_evolution/time_evolution.jl") @@ -105,7 +106,6 @@ include("time_evolution/ssesolve.jl") include("time_evolution/time_evolution_dynamical.jl") # Others -include("permutation.jl") include("correlations.jl") include("spectrum.jl") include("wigner.jl") diff --git a/src/permutation.jl b/src/permutation.jl deleted file mode 100644 index 32715790f..000000000 --- a/src/permutation.jl +++ /dev/null @@ -1,38 +0,0 @@ -export bdf, get_bdf_blocks - -function bdf(A::SparseMatrixCSC{T,M}) where {T,M} - n = LinearAlgebra.checksquare(A) - - G = DiGraph(abs.(A) .> 0) - idxs = connected_components(G) - P = sparse(1:n, reduce(vcat, idxs), ones(n), n, n) - block_sizes = map(length, idxs) - - return P, P * A * P', block_sizes -end - -function bdf( - A::QuantumObject{SparseMatrixCSC{T,M},OpType}, -) where {T,M,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} - P, A_bd, block_sizes = bdf(A.data) - return P, QuantumObject(A_bd, A.type, A.dims), block_sizes -end - -function get_bdf_blocks(A::SparseMatrixCSC{T,M}, block_sizes::Vector{Int}) where {T,M} - num_blocks = length(block_sizes) - block_indices = M[1] - block_list = [A[1:block_sizes[1], 1:block_sizes[1]]] - for i in 2:num_blocks - idx = sum(view(block_sizes, 1:i-1)) + 1 - push!(block_indices, idx) - push!(block_list, A[idx:idx-1+block_sizes[i], idx:idx-1+block_sizes[i]]) - end - return block_list, block_indices -end - -function get_bdf_blocks( - A::QuantumObject{SparseMatrixCSC{T,M},OpType}, - block_sizes::Vector{Int}, -) where {T,M,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} - return get_bdf_blocks(A.data, block_sizes) -end diff --git a/src/qobj/block_diagonal_form.jl b/src/qobj/block_diagonal_form.jl new file mode 100644 index 000000000..cc1f97eaa --- /dev/null +++ b/src/qobj/block_diagonal_form.jl @@ -0,0 +1,60 @@ +export BlockDiagonalForm, block_diagonal_form + +@doc raw""" + struct BlockDiagonalForm + +A type for storing a block-diagonal form of a matrix. + +# Fields +- `B::DT`: The block-diagonal matrix. It can be a sparse matrix or a [`QuantumObject`](@ref). +- `P::DT`: The permutation matrix. It can be a sparse matrix or a [`QuantumObject`](@ref). +- `blocks::AbstractVector`: The blocks of the block-diagonal matrix. +- `block_sizes::AbstractVector`: The sizes of the blocks. +""" +struct BlockDiagonalForm{DT,BT<:AbstractVector,BST<:AbstractVector} + B::DT + P::DT + blocks::BT + block_sizes::BST +end + +function block_diagonal_form(A::MT) where {MT<:AbstractSparseMatrix} + n = LinearAlgebra.checksquare(A) + + G = DiGraph(abs.(A)) + idxs_list = connected_components(G) + block_sizes = length.(idxs_list) + + P = MT(sparse(1:n, reduce(vcat, idxs_list), ones(n), n, n)) + + blocks = map(eachindex(idxs_list)) do i + m = block_sizes[i] + idxs = idxs_list[i] + P_i = MT(sparse(1:m, idxs, ones(m), m, n)) + return P_i * A * P_i' + end + + B = P * A * P' + + return BlockDiagonalForm(B, P, blocks, block_sizes) +end + +@doc raw""" + block_diagonal_form(A::QuantumObject) + +Return the block-diagonal form of a [`QuantumObject`](@ref). This is very useful in the presence of symmetries. + +# Arguments +- `A::QuantumObject`: The quantum object. + +# Returns +The [`BlockDiagonalForm`](@ref) of `A`. +""" +function block_diagonal_form( + A::QuantumObject{DT,OpType}, +) where {DT<:AbstractSparseMatrix,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} + bdf = block_diagonal_form(A.data) + B = QuantumObject(bdf.B, type = A.type, dims = A.dims) + P = QuantumObject(bdf.P, type = A.type, dims = A.dims) + return BlockDiagonalForm(B, P, bdf.blocks, bdf.block_sizes) +end diff --git a/test/core-test/permutation.jl b/test/core-test/block_diagonal_form.jl similarity index 62% rename from test/core-test/permutation.jl rename to test/core-test/block_diagonal_form.jl index 6f0892548..cc96a48d9 100644 --- a/test/core-test/permutation.jl +++ b/test/core-test/block_diagonal_form.jl @@ -1,4 +1,4 @@ -@testset "Permutation" begin +@testset "Block Diagonal Form" begin # Block Diagonal Form N = 20 Δ = 0 @@ -16,15 +16,17 @@ c_ops = [√(κ2) * a^2, √(κϕ) * ad * a] L = liouvillian(H, c_ops) - P, L_bd, block_sizes = bdf(L) - blocks_list, block_indices = get_bdf_blocks(L_bd, block_sizes) + bdf = block_diagonal_form(L) + L_bd = bdf.B + block_sizes = bdf.block_sizes + blocks = bdf.blocks + @test size(L_bd) == size(L) @test length(block_sizes) == 4 - @test length(blocks_list) == 4 - @test length(block_indices) == 4 + @test length(blocks) == 4 @test sum(block_sizes .== 100) == 4 - @testset "Type Inference (bdf)" begin - @inferred bdf(L) + @testset "Type Inference (block_diagonal_form)" begin + @inferred block_diagonal_form(L) end end diff --git a/test/runtests.jl b/test/runtests.jl index 502aa43cc..9081ad842 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ const testdir = dirname(@__FILE__) # Put core tests in alphabetical order core_tests = [ + "block_diagonal_form.jl", "correlations_and_spectrum.jl", "dynamical_fock_dimension_mesolve.jl", "dynamical-shifted-fock.jl", @@ -15,7 +16,6 @@ core_tests = [ "generalized_master_equation.jl", "low_rank_dynamics.jl", "negativity_and_partial_transpose.jl", - "permutation.jl", "progress_bar.jl", "quantum_objects.jl", "quantum_objects_evo.jl",