diff --git a/docs/make.jl b/docs/make.jl
index aa8d4c8..1cb651b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -12,10 +12,10 @@ makedocs(;
         edit_link="main",
         assets=String[]),
     pages=[
-        "Home" => "index.md",
+        "Home" => "index.md"
     ])
 
 deploydocs(;
     repo="github.com/tensor4all/Quantics.jl.git",
-    devbranch="main",
+    devbranch="main"
 )
diff --git a/src/binaryop.jl b/src/binaryop.jl
index 13ac71a..21004e0 100644
--- a/src/binaryop.jl
+++ b/src/binaryop.jl
@@ -14,8 +14,8 @@
   ``T_{x, y, \mathrm{out}, \mathrm{cin}, \mathrm{cout}} = 1`` if ``a  x + b  y + \mathrm{cin} = \mathrm{cout}``, ``=0`` otherwise (`out` is the output bit).
 """
 function _binaryop_tensor(a::Int, b::Int, site_x::Index{T}, site_y::Index{T},
-    site_out::Index{T},
-    cin_on::Bool, cout_on::Bool, bc::Int) where {T}
+        site_out::Index{T},
+        cin_on::Bool, cout_on::Bool, bc::Int) where {T}
     abs(a) <= 1 || error("a must be either 0, 1, -1")
     abs(b) <= 1 || error("b must be either 0, 1, -1")
     abs(bc) == 1 || error("bc must be either 1, -1")
@@ -47,11 +47,11 @@ end
 Create a tensor acting on a vector of sites.
 """
 function binaryop_tensor_multisite(sites::Vector{Index{T}},
-    coeffs::Vector{Tuple{Int,Int}},
-    pos_sites_in::Vector{Tuple{Int,Int}},
-    cin_on::Bool,
-    cout_on::Bool,
-    bc::Vector{Int}) where {T<:Number}
+        coeffs::Vector{Tuple{Int,Int}},
+        pos_sites_in::Vector{Tuple{Int,Int}},
+        cin_on::Bool,
+        cout_on::Bool,
+        bc::Vector{Int}) where {T<:Number}
 
     # Check
     sites = noprime.(sites)
@@ -122,11 +122,11 @@ where a, b, c, d = +/- 1, 0, and s1, s1 are arbitrary integers.
 `bc` is a vector of boundary conditions for each arguments of `g` (not of `f`).
 """
 function affinetransform(M::MPS,
-    tags::AbstractVector{String},
-    coeffs_dic::AbstractVector{Dict{String,Int}},
-    shift::AbstractVector{Int},
-    bc::AbstractVector{Int};
-    kwargs...)
+        tags::AbstractVector{String},
+        coeffs_dic::AbstractVector{Dict{String,Int}},
+        shift::AbstractVector{Int},
+        bc::AbstractVector{Int};
+        kwargs...)
     # f(x, y) = g(a * x + b * y + s1, c * x + d * y + s2)
     #         = h(a * x + b * y,      c * x + d * y),
     # where h(x, y) = g(x + s1, y + s2).
@@ -164,10 +164,10 @@ end
 
 # Version without shift
 function affinetransform(M::MPS,
-    tags::AbstractVector{String},
-    coeffs_dic::AbstractVector{Dict{String,Int}},
-    bc::AbstractVector{Int};
-    kwargs...)
+        tags::AbstractVector{String},
+        coeffs_dic::AbstractVector{Dict{String,Int}},
+        bc::AbstractVector{Int};
+        kwargs...)
 
     # f(x, y) = g(a * x + b * y + s1, c * x + d * y + s2)
     #         = h(a * x + b * y,      c * x + d * y),
@@ -297,10 +297,10 @@ f(x_R, y_R, ..., x_1, y_1) = M(x_R, y_R, ...; x'_R, y'_R, ...) f(x'_R, y'_R, ...
 `bc` is a vector of boundary conditions for each arguments of `g` (not of `f`).
 """
 function _binaryop_mpo(sites::Vector{Index{T}},
-    coeffs::Vector{Tuple{Int,Int}},
-    pos_sites_in::Vector{Tuple{Int,Int}};
-    rev_carrydirec=false,
-    bc::Union{Nothing,Vector{Int}}=nothing) where {T<:Number}
+        coeffs::Vector{Tuple{Int,Int}},
+        pos_sites_in::Vector{Tuple{Int,Int}};
+        rev_carrydirec=false,
+        bc::Union{Nothing,Vector{Int}}=nothing) where {T<:Number}
     # Number of variables involved in transformation
     nsites_bop = length(coeffs)
 
@@ -336,10 +336,10 @@ end
 
 # Limitation: a = -1 and b = -1 not supported. The same applies to (c, d).
 function _binaryop_mpo_backend(sites::Vector{Index{T}},
-    coeffs::Vector{Tuple{Int,Int}},
-    pos_sites_in::Vector{Tuple{Int,Int}};
-    rev_carrydirec=false,
-    bc::Union{Nothing,Vector{Int}}=nothing) where {T<:Number}
+        coeffs::Vector{Tuple{Int,Int}},
+        pos_sites_in::Vector{Tuple{Int,Int}};
+        rev_carrydirec=false,
+        bc::Union{Nothing,Vector{Int}}=nothing) where {T<:Number}
     nsites = length(sites)
     nsites_bop = length(coeffs)
     ncsites = nsites ÷ nsites_bop
diff --git a/src/fouriertransform.jl b/src/fouriertransform.jl
index bc1c6b9..a41516e 100644
--- a/src/fouriertransform.jl
+++ b/src/fouriertransform.jl
@@ -209,13 +209,13 @@ Instead of specifying `sitessrc`, one can specify the source sites by setting `t
 If `tag` = `x`, all sites with tags `x=1`, `x=2`, ... are used as `sitessrc`.
 """
 function fouriertransform(M::MPS;
-    sign::Int=1,
-    tag::String="",
-    sitessrc=nothing,
-    sitesdst=nothing,
-    originsrc::Float64=0.0,
-    origindst::Float64=0.0,
-    cutoff_MPO=1e-25, kwargs...)
+        sign::Int=1,
+        tag::String="",
+        sitessrc=nothing,
+        sitesdst=nothing,
+        originsrc::Float64=0.0,
+        origindst::Float64=0.0,
+        cutoff_MPO=1e-25, kwargs...)
     sites = siteinds(M)
     sitepos, target_sites = _find_target_sites(M; sitessrc=sitessrc, tag=tag)
 
diff --git a/src/imaginarytime.jl b/src/imaginarytime.jl
index 55fd37c..39827d5 100644
--- a/src/imaginarytime.jl
+++ b/src/imaginarytime.jl
@@ -17,7 +17,7 @@ _stat_sign(::Fermionic) = -1
 _stat_sign(::Bosonic) = 1
 
 function to_wn(stat::Statistics, gtau::MPS, beta::Float64; sitessrc=nothing, tag="",
-    sitesdst=nothing, kwargs...)::MPS
+        sitesdst=nothing, kwargs...)::MPS
     sitepos, _ = _find_target_sites(gtau; sitessrc=sitessrc, tag=tag)
     nqbit_t = length(sitepos)
     originwn = 0.5 * (-2.0^nqbit_t + _stat_shift(stat))
@@ -28,7 +28,7 @@ function to_wn(stat::Statistics, gtau::MPS, beta::Float64; sitessrc=nothing, tag
 end
 
 function to_tau(stat::Statistics, giv::MPS, beta::Float64; sitessrc=nothing, tag="",
-    sitesdst=nothing, kwargs...)::MPS
+        sitesdst=nothing, kwargs...)::MPS
     sitepos, _ = _find_target_sites(giv; sitessrc=sitessrc, tag=tag)
     nqbit_t = length(sitepos)
     originwn = 0.5 * (-2.0^nqbit_t + _stat_shift(stat))
diff --git a/src/mul.jl b/src/mul.jl
index 87f5f94..af11a6b 100644
--- a/src/mul.jl
+++ b/src/mul.jl
@@ -9,15 +9,15 @@ struct MatrixMultiplier{T} <: AbstractMultiplier where {T}
     sites_col::Vector{Index{T}}
 
     function MatrixMultiplier(sites_row::Vector{Index{T}},
-        sites_shared::Vector{Index{T}},
-        sites_col::Vector{Index{T}}) where {T}
+            sites_shared::Vector{Index{T}},
+            sites_col::Vector{Index{T}}) where {T}
         new{T}(sites_row, sites_shared, sites_col)
     end
 end
 
 function MatrixMultiplier(site_row::Index{T},
-    site_shared::Index{T},
-    site_col::Index{T}) where {T}
+        site_shared::Index{T},
+        site_col::Index{T}) where {T}
     return MatrixMultiplier([site_row], [site_shared], [site_col])
 end
 
@@ -128,7 +128,7 @@ end
 By default, elementwise multiplication will be performed.
 """
 function automul(M1::MPS, M2::MPS; tag_row::String="", tag_shared::String="",
-    tag_col::String="", alg="naive", kwargs...)
+        tag_col::String="", alg="naive", kwargs...)
     if in(:maxbonddim, keys(kwargs))
         error("Illegal keyward parameter: maxbonddim. Use maxdim instead!")
     end
diff --git a/src/tag.jl b/src/tag.jl
index 4b47cfe..e17a163 100644
--- a/src/tag.jl
+++ b/src/tag.jl
@@ -12,7 +12,7 @@ If not, the function seach for all Index objects with tags `x=1`, `x=2`, ..., an
 If no Index object is found, an empty vector will be returned.
 """
 function findallsites_by_tag(sites::Vector{Index{T}}; tag::String="x",
-    maxnsites::Int=1000)::Vector{Int} where {T}
+        maxnsites::Int=1000)::Vector{Int} where {T}
     _valid_tag(tag) || error("Invalid tag: $tag")
     result = Int[]
     for n in 1:maxnsites
@@ -28,8 +28,36 @@ function findallsites_by_tag(sites::Vector{Index{T}}; tag::String="x",
     return result
 end
 
-function findallsiteinds_by_tag(sites; tag::String="x", maxnsites::Int=1000)
+function findallsiteinds_by_tag(
+        sites::AbstractVector{Index{T}}; tag::String="x", maxnsites::Int=1000) where {T}
     _valid_tag(tag) || error("Invalid tag: $tag")
     positions = findallsites_by_tag(sites; tag=tag, maxnsites=maxnsites)
     return [sites[p] for p in positions]
 end
+
+function findallsites_by_tag(sites::Vector{Vector{Index{T}}}; tag::String="x",
+        maxnsites::Int=1000)::Vector{NTuple{2,Int}} where {T}
+    _valid_tag(tag) || error("Invalid tag: $tag")
+
+    sites_dict = Dict{Index{T},NTuple{2,Int}}()
+    for i in 1:length(sites)
+        for j in 1:length(sites[i])
+            sites_dict[sites[i][j]] = (i, j)
+        end
+    end
+
+    result = NTuple{2,Int}[]
+    sitesflatten = collect(Iterators.flatten(sites))
+    for n in 1:maxnsites
+        tag_ = tag * "=$n"
+        idx = findall(i -> hastags(i, tag_) && hasplev(i, 0), sitesflatten)
+        if length(idx) == 0
+            break
+        elseif length(idx) > 1
+            error("Found more than one site indices with $(tag_)!")
+        end
+
+        push!(result, sites_dict[sitesflatten[only(idx)]])
+    end
+    return result
+end
diff --git a/src/transformer.jl b/src/transformer.jl
index 1a1a588..c8e8640 100644
--- a/src/transformer.jl
+++ b/src/transformer.jl
@@ -19,7 +19,7 @@ f(x) = g(-x)
 where f(x) = M * g(x) for x = 0, 1, ..., 2^R-1.
 """
 function flipop_to_negativedomain(sites::Vector{Index{T}}; rev_carrydirec=false,
-    bc::Int=1)::MPO where {T}
+        bc::Int=1)::MPO where {T}
     return flipop(sites; rev_carrydirec=rev_carrydirec, bc=bc) * bc
 end
 
@@ -146,7 +146,7 @@ end
 Create QTT for a upper/lower triangle matrix filled with one except the diagonal line
 """
 function upper_lower_triangle_matrix(sites::Vector{Index{T}}, value::S;
-    upper_or_lower::Symbol=:upper)::MPO where {T,S}
+        upper_or_lower::Symbol=:upper)::MPO where {T,S}
     upper_or_lower ∈ [:upper, :lower] || error("Invalid upper_or_lower $(upper_or_lower)")
     N = length(sites)
 
diff --git a/src/util.jl b/src/util.jl
index 0c084be..bbac3d5 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -93,7 +93,7 @@ new_sites: Vector of vectors of new siteinds
 When splitting MPS tensors, the column major is assumed.
 """
 function unfuse_siteinds(M::MPS, targetsites::Vector{Index{T}},
-    newsites::AbstractVector{Vector{Index{T}}})::MPS where {T}
+        newsites::AbstractVector{Vector{Index{T}}})::MPS where {T}
     length(targetsites) == length(newsites) || error("Length mismatch")
     links = linkinds(M)
     L = length(M)
@@ -437,7 +437,8 @@ end
 function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MPS where {T}
     sitesold = siteinds(MPO(collect(M)))
 
-    Set(Iterators.flatten(sites)) == Set(Iterators.flatten(sitesold)) || error("siteinds do not match $(sites) != $(sitesold)")
+    Set(Iterators.flatten(sites)) == Set(Iterators.flatten(sitesold)) ||
+        error("siteinds do not match $(sites) != $(sitesold)")
 
     t = ITensor(1)
     tensors = Vector{ITensor}(undef, length(sites))
@@ -447,7 +448,7 @@ function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MP
             if ind ∈ inds(t)
                 continue
             end
-            contract_until = findfirst(x->ind ∈ Set(collect(x)), inds.(tensors_old))
+            contract_until = findfirst(x -> ind ∈ Set(collect(x)), inds.(tensors_old))
             contract_until !== nothing || error("ind $ind not found")
             for j in 1:contract_until
                 t *= tensors_old[j]
@@ -457,9 +458,8 @@ function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MP
             end
         end
 
-        linds = 
-        if i > 1
-            vcat(only(commoninds(t, tensors[i-1])), sites[i])
+        linds = if i > 1
+            vcat(only(commoninds(t, tensors[i - 1])), sites[i])
         else
             sites[i]
         end
@@ -467,4 +467,59 @@ function rearrange_siteinds(M::AbstractMPS, sites::Vector{Vector{Index{T}}})::MP
     end
     tensors[end] *= t
     cleanup_linkinds!(MPS(tensors))
-end
\ No newline at end of file
+end
+
+"""
+Makes an MPS/MPO diagonal for a specified a site index `s`.
+On return, the data will be deep copied and the target core tensor will be diagonalized with an additional site index `s'`.
+"""
+function makesitediagonal(M::AbstractMPS, site::Index{T})::MPS where {T}
+    M_ = deepcopy(MPO(collect(M)))
+
+    target_site::Int = only(findsites(M_, site))
+    M_[target_site] = _asdiagonal(M_[target_site], site)
+
+    return MPS(collect(M_))
+end
+
+function makesitediagonal(M::AbstractMPS, tag::String)::MPS
+    M_ = deepcopy(MPO(collect(M)))
+    sites = siteinds(M_)
+
+    target_positions = findallsites_by_tag(siteinds(M_); tag=tag)
+
+    for t in eachindex(target_positions)
+        i, j = target_positions[t]
+        M_[i] = _asdiagonal(M_[i], sites[i][j])
+    end
+
+    return MPS(collect(M_))
+end
+
+"""
+Extract diagonal components
+"""
+function extractdiagonal(M::AbstractMPS, tag::String)::MPS
+    M_ = deepcopy(MPO(collect(M)))
+    sites = siteinds(M_)
+
+    target_positions = findallsites_by_tag(siteinds(M_); tag=tag)
+
+    for t in eachindex(target_positions)
+        i, j = target_positions[t]
+        M_[i] = _extract_diagonal(M_[i], sites[i][j], sites[i][j]')
+    end
+
+    return MPS(collect(M_))
+end
+
+function _extract_diagonal(t, site::Index{T}, site2::Index{T}) where {T<:Number}
+    dim(site) == dim(site2) || error("Dimension mismatch")
+    restinds = uniqueinds(inds(t), site, site2)
+    newdata = zeros(eltype(t), dim.(restinds)..., dim(site))
+    olddata = Array(t, restinds..., site, site2)
+    for i in 1:dim(site)
+        newdata[.., i] = olddata[.., i, i]
+    end
+    return ITensor(newdata, restinds..., site)
+end
diff --git a/test/test_with_aqua.jl b/test/test_with_aqua.jl
index 61db0de..0a8888c 100644
--- a/test/test_with_aqua.jl
+++ b/test/test_with_aqua.jl
@@ -5,4 +5,4 @@
     @testset "Aqua" begin
         Aqua.test_stale_deps(Quantics)
     end
-end
\ No newline at end of file
+end
diff --git a/test/transformer_tests.jl b/test/transformer_tests.jl
index 272009a..6096f68 100644
--- a/test/transformer_tests.jl
+++ b/test/transformer_tests.jl
@@ -76,6 +76,8 @@ end
             f_ref[i] = g_reconst[mod(2^nbit - (i - 1), 2^nbit) + 1]
         end
         f_ref[1] *= bc
+
+        @test f_reconst ≈ f_ref
     end
 
     @testset "reverseaxis2" for nbit in 2:3, rev_carrydirec in [true, false]
diff --git a/test/util_tests.jl b/test/util_tests.jl
index 649fc87..6dc5325 100644
--- a/test/util_tests.jl
+++ b/test/util_tests.jl
@@ -224,4 +224,25 @@
 
         @test Ψ ≈ Ψ_reconst
     end
-end
\ No newline at end of file
+
+    @testset "makesitediagonal" begin
+        L = 2
+
+        sitesx = [Index(2, "x=$n") for n in 1:L]
+
+        Ψ = random_mps(sitesx)
+
+        M = Quantics.makesitediagonal(Ψ, "x")
+
+        Ψ_recost = Array(prod(Ψ), sitesx...)
+        M_recost = Array(prod(M), prime.(sitesx)..., sitesx...)
+
+        for i in 1:2, i2 in 1:2, j in 1:2, j2 in 1:2
+            if i != i2 || j != j2
+                @test M_recost[i, j, i2, j2] ≈ 0.0
+            else
+                @test M_recost[i, j, i2, j2] ≈ Ψ_recost[i, j]
+            end
+        end
+    end
+end