add permute systems for hermitian/symmetric matrices and in-place ver…

…sion for vectors

docs/src/api.md

@@ -13,6 +13,7 @@ gell_mann
 gell_mann!
 partial_trace
 partial_transpose
+permute_systems!
 permute_systems
 cleanup!
 ```
@@ -106,4 +107,5 @@ Ket._partition
 Ket._fiducial_WH
 Ket._idx
 Ket._tidx
+Ket._idxperm
 ```

src/multilinear.jl

@@ -202,47 +202,75 @@ Takes the partial transpose of matrix `X` on the subsystems in `transp` assuming
 partial_transpose(X::AbstractMatrix, transp) = partial_transpose(X, transp, _equal_sizes(X))
 
 """
-    permute_systems(X::AbstractVector, perm::Vector, dims::Vector)
+    _idxperm(perm::Vector, dims::Vector)
 
-Permutes the order of the subsystems of vector `X` with subsystem dimensions `dims` according to the permutation `perm`.
+Computes the index permutation associated with permuting the subsystems of a vector with subsystem dimensions `dims` according to `perm`.
 """
-function permute_systems(X::AbstractVector{T}, perm::Vector{<:Integer}, dims::Vector{<:Integer}) where {T}
-    perm == 1:length(perm) && return X
-
-    dX = length(X)
-    dimsY = dims[perm]
+function _idxperm(perm::Vector{<:Integer}, dims::Vector{<:Integer})
+    p = Vector{Int64}(undef, prod(dims))
+    _idxperm!(p, perm, dims)
+    return p
+end
 
-    Y = Vector{T}(undef, dX)
+function _idxperm!(p::Vector{<:Integer}, perm::Vector{<:Integer}, dims::Vector{<:Integer})
+    pdims = dims[perm]
 
     ti = Vector{Int64}(undef, length(dims))
 
-    for i in 1:dX
+    for i in eachindex(p)
         _tidx!(ti, i, dims)
         permute!(ti, perm)
-        Yi = _idx(ti, dimsY)
-        Y[Yi] = X[i]
+        j = _idx(ti, pdims)
+        p[j] = i
     end
-    return Y
+    return p
 end
 
 """
-    permute_systems(X::AbstractMatrix, perm::Vector, dims::Vector)
+    permute_systems!(X::AbstractVector, perm::Vector, dims::Vector)
 
-Permutes the order of the subsystems of the square matrix `X`, which is composed by square subsystems of dimensions `dims`, according to the permutation `perm`.
+Permutes the order of the subsystems of vector `X` with subsystem dimensions `dims` in-place according to the permutation `perm`.
 """
-function permute_systems(X::AbstractMatrix, perm::Vector{<:Integer}, dims::Vector{<:Integer})
+function permute_systems!(X::AbstractVector{T}, perm::Vector{<:Integer}, dims::Vector{<:Integer}) where {T}
     perm == 1:length(perm) && return X
+    X == 1:length(X) && return _idxperm!(X, perm, dims)
 
-    idxperm = permute_systems(axes(X, 1), perm, dims)
-    return X[idxperm, idxperm]
+    p = _idxperm(perm, dims)
+    permute!(X, p)
+    return X 
 end
+export permute_systems!
 
 """
+    permute_systems!(X::AbstractVector, perm::Vector)
+
+Permutes the order of the subsystems of vector `X`, which is composed by two subsystems of equal dimensions, in-place according to the permutation `perm`.
+"""
+permute_systems!(X::AbstractVector, perm::Vector{<:Integer}) = permute_systems!(X, perm, _equal_sizes(X))
+
+@doc"""
     permute_systems(X::AbstractMatrix, perm::Vector, dims::Vector)
 
+Permutes the order of the subsystems of the square matrix `X`, which is composed by square subsystems of dimensions `dims`, according to the permutation `perm`.
+""" permute_systems(X::AbstractMatrix, perm::Vector, dims::Vector)
+for (T, wrapper) in
+    [(:AbstractMatrix, :identity), (:(LA.Hermitian), :(LA.Hermitian)), (:(LA.Symmetric), :(LA.Symmetric))]
+    @eval begin
+        function permute_systems(X::$T, perm::Vector{<:Integer}, dims::Vector{<:Integer})
+            perm == 1:length(perm) && return X
+
+	        p = _idxperm(perm, dims)
+	        return $wrapper(X[p, p])
+	    end
+    end
+end
+
+"""
+    permute_systems(X::AbstractMatrix, perm::Vector)
+
 Permutes the order of the subsystems of the square matrix `X`, which is composed by two square subsystems of equal dimensions, according to the permutation `perm`.
 """
-permute_systems(X::AbstractVecOrMat, perm::Vector{<:Integer}) = permute_systems(X, perm, _equal_sizes(X))
+permute_systems(X::AbstractMatrix, perm::Vector{<:Integer}) = permute_systems(X, perm, _equal_sizes(X))
 
 """
     permute_systems(X::AbstractMatrix, perm::Vector, dims::Matrix)
@@ -253,9 +281,8 @@ Permutes the order of the subsystems of the matrix `X`, which is composed by sub
 function permute_systems(X::AbstractMatrix, perm::Vector{<:Integer}, dims::Matrix{<:Integer})
     perm == 1:length(perm) && return X
 
-    rowperm = permute_systems(axes(X, 1), perm, dims[:, 1])
-    colperm = permute_systems(axes(X, 2), perm, dims[:, 2])
-
-    return X[rowperm, colperm]
+    rowp = _idxperm(perm, dims[:, 1])
+    colp = _idxperm(perm, dims[:, 2])
+    return X[rowp, colp]
 end
 export permute_systems

test/multilinear.jl

@@ -76,17 +76,16 @@
                     uv = kron(u, v)
                     uw = kron(u, w)
                     vw = kron(v, w)
-                    uvw = kron(u, v, w)
-                    @test permute_systems(uv, [1, 2]) ≈ kron(u, v)
-                    @test permute_systems(uv, [2, 1]) ≈ kron(v, u)
-                    @test permute_systems(uw, [2, 1], [d1, d3]) ≈ kron(w, u)
-                    @test permute_systems(vw, [2, 1], [d2, d3]) ≈ kron(w, v)
-                    @test permute_systems(uvw, [1, 2, 3], [d1, d2, d3]) ≈ kron(u, v, w)
-                    @test permute_systems(uvw, [2, 3, 1], [d1, d2, d3]) ≈ kron(v, w, u)
-                    @test permute_systems(uvw, [3, 1, 2], [d1, d2, d3]) ≈ kron(w, u, v)
-                    @test permute_systems(uvw, [1, 3, 2], [d1, d2, d3]) ≈ kron(u, w, v)
-                    @test permute_systems(uvw, [2, 1, 3], [d1, d2, d3]) ≈ kron(v, u, w)
-                    @test permute_systems(uvw, [3, 2, 1], [d1, d2, d3]) ≈ kron(w, v, u)
+		            @test permute_systems!(kron(u,v), [1, 2]) ≈ kron(u, v)
+		            @test permute_systems!(kron(u,v), [2, 1]) ≈ kron(v, u)
+                    @test permute_systems!(kron(u,v), [2, 1]) ≈ kron(v, u)
+		            @test permute_systems!(kron(v, w), [2, 1], [d2, d3]) ≈ kron(w, v)
+		            @test permute_systems!(kron(u,v,w), [1, 2, 3], [d1, d2, d3]) ≈ kron(u, v, w)
+                    @test permute_systems!(kron(u,v,w),[2, 3, 1], [d1, d2, d3]) ≈ kron(v, w, u)
+                    @test permute_systems!(kron(u,v,w), [3, 1, 2], [d1, d2, d3]) ≈ kron(w, u, v)
+		            @test permute_systems!(kron(u,v,w), [1, 3, 2], [d1, d2, d3]) ≈ kron(u, w, v)
+		            @test permute_systems!(kron(u,v,w), [2, 1, 3], [d1, d2, d3]) ≈ kron(v, u, w)
+		            @test permute_systems!(kron(u,v,w), [3, 2, 1], [d1, d2, d3]) ≈ kron(w, v, u)
                 end
             end
         end
@@ -114,6 +113,12 @@
                     @test permute_systems(abc, [3, 2, 1], [d1, d2, d3]) ≈ kron(c, b, a)
                 end
             end
+	        for wrapper in [Symmetric, Hermitian]
+		        M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3)))
+		        x = Matrix(M)
+		        @test permute_systems(M, [3, 1, 2], [d1, d2, d3]) ≈ permute_systems(x, [3, 1, 2], [d1, d2, d3])
+		        @test permute_systems(M, [1, 3, 2], [d1, d2, d3]) ≈ permute_systems(x, [1, 3, 2], [d1, d2, d3])
+	        end
         end
 
         @testset "Rectangular matrices" begin