Change block diagonal form structure (#349)

* Change block diagonal form structure

* Add changelog

Author: albertomercurio

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

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

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")

src/permutation.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

src/qobj/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

test/core-test/permutation.jl renamed to test/core-test/block_diagonal_form.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",