diff --git a/.gitignore b/.gitignore
index 21d120e..ea57bb3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ Manifest.toml
 *.DS_Store
 *.pb.gz
 *.pdf
+.history
diff --git a/src/Lehmann.jl b/src/Lehmann.jl
index e739a4e..103a4dd 100644
--- a/src/Lehmann.jl
+++ b/src/Lehmann.jl
@@ -8,16 +8,97 @@ include("utility/chebyshev.jl")
 include("spectral.jl")
 export Spectral
 
-include("sample/sample.jl")
-export Sample
-
 include("discrete/builder.jl")
 include("functional/builder.jl")
 
 include("dlr.jl")
 export DLRGrid
 
+include("sample/sample.jl")
+export Sample
+
 include("operation.jl")
 export tau2dlr, dlr2tau, matfreq2dlr, dlr2matfreq, tau2matfreq, matfreq2tau, tau2tau, matfreq2matfreq
 
+##################### precompile #######################
+# precompile as the final step of the module definition:
+if ccall(:jl_generating_output, Cint, ()) == 1   # if we're precompiling the package
+    let
+        #cover data type and symmetry
+        # Float32, :none
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :none; dtype=Float32)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+
+        # Float32, :pha
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :pha; dtype=Float32)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+        # Float32, :ph
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :ph; dtype=Float32)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+        # Float64, :none
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :none; dtype=Float64)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+        # Float64, :pha
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :pha; dtype=Float64)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+        # Float64, :ph
+        dlr = DLRGrid(1.0, 1.0, 1e-6, true, :ph; dtype=Float64)
+        g = Sample.SemiCircle(dlr, :τ, dlr.τ)
+        coeff = tau2dlr(dlr, g)
+        _g = dlr2tau(dlr, coeff)
+        _g = tau2tau(dlr, g, dlr.τ)
+
+        g = Sample.SemiCircle(dlr, :n, dlr.n)
+        coeff = matfreq2dlr(dlr, g)
+        _g = dlr2matfreq(dlr, coeff)
+        _g = matfreq2matfreq(dlr, g, dlr.n)
+
+    end
+end
+
 end
diff --git a/src/discrete/kernel.jl b/src/discrete/kernel.jl
index 5fa1913..3c9cf07 100644
--- a/src/discrete/kernel.jl
+++ b/src/discrete/kernel.jl
@@ -103,11 +103,11 @@ function preciseKernelT(dlrGrid, τ, ω, print::Bool = true)
         #symmetrize K(τ, ω)=K(β-τ, -ω) for τ>0 
         @assert isodd(τ.np) #symmetrization is only possible for odd τ panels
         halfτ = ((τ.np - 1) ÷ 2) * τ.degree
-        kernel[1:halfτ, :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid, 1.0, true)
-        kernel[end:-1:(halfτ+1), :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid[end:-1:1], 1.0, true)
+        kernel[1:halfτ, :] = Spectral.kernelT(Float64, Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid, 1.0, true)
+        kernel[end:-1:(halfτ+1), :] = Spectral.kernelT(Float64, Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid[end:-1:1], 1.0, true)
         # use the fermionic kernel for both the fermionic and bosonic propagators
     else
-        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(symmetry), τGrid, ωGrid, 1.0, true)
+        kernel = Spectral.kernelT(Float64, Val(dlrGrid.isFermi), Val(symmetry), τGrid, ωGrid, 1.0, true)
     end
 
     # print && println("=====  Kernel Discretization =====")
@@ -191,8 +191,8 @@ function preciseKernelΩn(dlrGrid, ω, print::Bool = true)
     symmetry = dlrGrid.symmetry
     n = nGrid(dlrGrid.isFermi, symmetry, dlrGrid.Λ)
 
-    nkernelFermi = Spectral.kernelΩ(Val(true), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
-    nkernelBose = Spectral.kernelΩ(Val(false), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
+    nkernelFermi = Spectral.kernelΩ(Float64, Val(true), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
+    nkernelBose = Spectral.kernelΩ(Float64, Val(false), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
 
     return n, nkernelFermi, nkernelBose
 end
\ No newline at end of file
diff --git a/src/dlr.jl b/src/dlr.jl
index 4347c30..92d64c9 100644
--- a/src/dlr.jl
+++ b/src/dlr.jl
@@ -26,145 +26,32 @@ struct DLRGrid
 - `ωn` or `omegaN` : (2n+1)π/β
 - `τ` or `tau` : selected representative imaginary-time grid
 """
-mutable struct DLRGrid
+mutable struct DLRGrid{T<:Real,S}
     isFermi::Bool
     symmetry::Symbol
-    Euv::Float64
-    β::Float64
-    Λ::Float64
-    rtol::Float64
+    Euv::T
+    β::T
+    Λ::T
+    rtol::T
 
     # dlr grids
     # size::Int # rank of the dlr representation
-    ω::Vector{Float64}
+    ω::Vector{T}
     n::Vector{Int} # integers, (2n+1)π/β gives the Matsubara frequency
-    ωn::Vector{Float64} # (2n+1)π/β
-    τ::Vector{Float64}
-
-    kernel_τ::Any
-    kernel_n::Any
-
-    """
-    function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
-    function DLRGrid(; isFermi::Bool, β = -1.0, beta = -1.0, Euv = 1.0, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
-
-    Create DLR grids
-
-    #Arguments:
-    - `Euv`         : the UV energy scale of the spectral density 
-    - `β` or `beta` : inverse temeprature
-    - `isFermi`     : bool is fermionic or bosonic
-    - `symmetry`    : particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
-    - `rtol`        : tolerance absolute error
-    - `rebuild`     : set false to load DLR basis from the file, set true to recalculate the DLR basis on the fly
-    - `folder`      : if rebuild is true and folder is set, then dlrGrid will be rebuilt and saved to the specified folder
-                      if rebuild is false and folder is set, then dlrGrid will be loaded from the specified folder
-    - `algorithm`   : if rebuild = true, then set :functional to use the functional algorithm to generate the DLR basis, or set :discrete to use the matrix algorithm.
-    - `verbose`     : false not to print DLRGrid to terminal, or true to print
-    """
-    function DLRGrid(Euv, β, rtol, isFermi::Bool, symmetry::Symbol = :none; rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
-        Λ = Euv * β # dlr only depends on this dimensionless scale
-        # println("Get $Λ")
-        @assert rtol > 0.0 "rtol=$rtol is not positive and nonzero!"
-        @assert Λ > 0 "Energy scale $Λ must be positive!"
-        @assert symmetry == :ph || symmetry == :pha || symmetry == :none "symmetry must be :ph, :pha or nothing"
-        @assert algorithm == :functional || algorithm == :discrete "Algorithm is either :functional or :discrete"
-        @assert β > 0.0 "Inverse temperature must be temperature."
-        @assert Euv > 0.0 "Energy cutoff must be positive."
-
-        if Λ > 1e8 && symmetry == :none
-            @warn("Current DLR without symmetry may cause ~ 3-4 digits loss for Λ ≥ 1e8!")
-        end
+    ωn::Vector{T} # (2n+1)π/β
+    τ::Vector{T}
 
-        if rtol > 1e-6
-            @warn("Current implementation may cause ~ 3-4 digits loss for rtol > 1e-6!")
-        end
+    kernel_τ::Matrix{T}
+    kernel_n::Matrix{T}
+    kernel_nc::Matrix{Complex{T}}
 
-        rtolpower = Int(floor(log10(rtol))) # get the biggest n so that rtol>1e-n
-        if abs(rtolpower) < 4
-            rtolpower = -4
-        end
-        rtol = 10.0^float(rtolpower)
-
-        function finddlr(folder, filename)
-            searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
-            for dir in folder
-                if length(searchdir(dir, filename)) > 0
-                    #dlr file found
-                    return joinpath(dir, filename)
-                end
-            end
-            @warn("Cann't find the DLR file $filename in the folders $folder. Regenerating DLR...")
-            return nothing
-        end
-
-        function filename(lambda, errpower)
-            lambda = Int128(floor(lambda))
-            errstr = "1e$errpower"
-
-            if symmetry == :none
-                return "universal_$(lambda)_$(errstr).dlr"
-            elseif symmetry == :ph
-                return "ph_$(lambda)_$(errstr).dlr"
-            elseif symmetry == :pha
-                return "pha_$(lambda)_$(errstr).dlr"
-            else
-                error("$symmetry is not implemented!")
-            end
-        end
-
-
-        if rebuild == false
-            if isnothing(folder)
-                Λ = Λ < 100 ? Int(100) : 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n
-                folderList = [string(@__DIR__, "/../basis/"),]
-            else
-                folderList = [folder,]
-            end
-
-            file = filename(Λ, rtolpower)
-            dlrfile = finddlr(folderList, file)
-
-            if isnothing(dlrfile) == false
-                dlr = new(isFermi, symmetry, Euv, β, Λ, rtol, [], [], [], [], nothing, nothing)
-                _load!(dlr, dlrfile, verbose)
-                dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
-                dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
-                return dlr
-            else
-                @warn("No DLR is found in the folder $folder, try to rebuild instead.")
-            end
-
-        end
-
-        # try to rebuild the dlrGrid
-        dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], [], nothing, nothing)
-        file2save = filename(Euv * β, rtolpower)
-        _build!(dlr, folder, file2save, algorithm, verbose)
-
-        dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
-        dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
-        return dlr
-    end
-
-    function DLRGrid(; isFermi::Bool, β = -1.0, beta = -1.0, Euv = 1.0, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
-        if β <= 0.0 && beta > 0.0
-            β = beta
-        elseif β > 0.0 && beta <= 0.0
-            beta = β
-        elseif β < 0.0 && beta < 0.0
-            error("Either β or beta needs to be initialized with a positive value!")
-        end
-        @assert β ≈ beta
-        return DLRGrid(Euv, β, rtol, isFermi, symmetry; rebuild = rebuild, folder = folder, algorithm = algorithm, verbose = verbose)
-    end
 end
 
 function Base.getproperty(obj::DLRGrid, sym::Symbol)
     # if sym === :hasTau
     #     return obj.totalTauNum > 0
     if sym == :size
-        return size(obj)
+        return length(obj)
     elseif sym == :tau
         return obj.τ
     elseif sym == :beta
@@ -180,21 +67,142 @@ function Base.getproperty(obj::DLRGrid, sym::Symbol)
     end
 end
 
+"""
+function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
+function DLRGrid(; isFermi::Bool, β = -1.0, beta = -1.0, Euv = 1.0, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
+
+Create DLR grids
+
+#Arguments:
+- `Euv`         : the UV energy scale of the spectral density 
+- `β` or `beta` : inverse temeprature
+- `isFermi`     : bool is fermionic or bosonic
+- `symmetry`    : particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
+- `rtol`        : tolerance absolute error
+- `rebuild`     : set false to load DLR basis from the file, set true to recalculate the DLR basis on the fly
+- `folder`      : if rebuild is true and folder is set, then dlrGrid will be rebuilt and saved to the specified folder
+                    if rebuild is false and folder is set, then dlrGrid will be loaded from the specified folder
+- `algorithm`   : if rebuild = true, then set :functional to use the functional algorithm to generate the DLR basis, or set :discrete to use the matrix algorithm.
+- `verbose`     : false not to print DLRGrid to terminal, or true to print
+"""
+function DLRGrid(Euv, β, rtol, isFermi::Bool, symmetry::Symbol=:none;
+    rebuild=false, folder=nothing, algorithm=:functional, verbose=false, dtype=Float64)
+
+    T = dtype
+
+    Λ = Euv * β # dlr only depends on this dimensionless scale
+    # println("Get $Λ")
+    @assert rtol > 0.0 "rtol=$rtol is not positive and nonzero!"
+    @assert Λ > 0 "Energy scale $Λ must be positive!"
+    @assert symmetry == :ph || symmetry == :pha || symmetry == :none "symmetry must be :ph, :pha or nothing"
+    @assert algorithm == :functional || algorithm == :discrete "Algorithm is either :functional or :discrete"
+    @assert β > 0.0 "Inverse temperature must be temperature."
+    @assert Euv > 0.0 "Energy cutoff must be positive."
+
+    if Λ > 1e8 && symmetry == :none
+        @warn("Current DLR without symmetry may cause ~ 3-4 digits loss for Λ ≥ 1e8!")
+    end
+
+    if rtol > 1e-6
+        @warn("Current implementation may cause ~ 3-4 digits loss for rtol > 1e-6!")
+    end
+
+    rtolpower = Int(floor(log10(rtol))) # get the biggest n so that rtol>1e-n
+    if abs(rtolpower) < 4
+        rtolpower = -4
+    end
+    rtol = T(10.0)^T(rtolpower)
+
+    function finddlr(folder, filename)
+        searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
+        for dir in folder
+            if length(searchdir(dir, filename)) > 0
+                #dlr file found
+                return joinpath(dir, filename)
+            end
+        end
+        @warn("Cann't find the DLR file $filename in the folders $folder. Regenerating DLR...")
+        return nothing
+    end
+
+    function filename(lambda, errpower)
+        lambda = Int128(floor(lambda))
+        errstr = "1e$errpower"
+
+        if symmetry == :none
+            return "universal_$(lambda)_$(errstr).dlr"
+        elseif symmetry == :ph
+            return "ph_$(lambda)_$(errstr).dlr"
+        elseif symmetry == :pha
+            return "pha_$(lambda)_$(errstr).dlr"
+        else
+            error("$symmetry is not implemented!")
+        end
+    end
+
+
+    if rebuild == false
+        if isnothing(folder)
+            Λ = Λ < 100 ? Int(100) : 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n
+            folderList = [string(@__DIR__, "/../basis/"),]
+        else
+            folderList = [folder,]
+        end
+
+        file = filename(Λ, rtolpower)
+        dlrfile = finddlr(folderList, file)
+
+        if isnothing(dlrfile) == false
+            dlr = DLRGrid{T,symmetry}(isFermi, symmetry, Euv, β, Λ, rtol, [], [], [], [], zeros(T, 1, 1), zeros(T, 1, 1), zeros(Complex{T}, 1, 1))
+            _load!(dlr, dlrfile, verbose)
+            # dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
+            # dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
+            return dlr
+        else
+            @warn("No DLR is found in the folder $folder, try to rebuild instead.")
+        end
+
+    end
+
+    # try to rebuild the dlrGrid
+    dlr = DLRGrid{T,symmetry}(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], [], zeros(T, 1, 1), zeros(T, 1, 1), zeros(Complex{T}, 1, 1))
+    file2save = filename(Euv * β, rtolpower)
+    _build!(dlr, folder, file2save, algorithm, verbose)
+
+    # dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
+    # dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
+    return dlr
+end
+function DLRGrid(; isFermi::Bool, β=-1.0, beta=-1.0, Euv=1.0, symmetry::Symbol=:none,
+    rtol=1e-14, rebuild=false, folder=nothing, algorithm=:functional, verbose=false, dtype=Float64)
+    T = dtype
+    if β <= T(0) && beta > T(0)
+        β = beta
+    elseif β > T(0) && beta <= T(0)
+        beta = β
+    elseif β < T(0) && beta < T(0)
+        error("Either β or beta needs to be initialized with a positive value!")
+    end
+    @assert β ≈ beta
+    return DLRGrid(Euv, β, rtol, isFermi, symmetry; rebuild=rebuild, folder=folder, algorithm=algorithm, verbose=verbose, dtype=dtype)
+end
+
 
 """
-Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω)
+Base.size(dlrGrid::DLRGrid) = (length(dlrGrid.ω),)
 Base.length(dlrGrid::DLRGrid) = length(dlrGrid.ω)
 rank(dlrGrid::DLRGrid) = length(dlrGrid.ω)
 
 get the rank of the DLR grid, namely the number of the DLR coefficients.
 """
-Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω)
+Base.size(dlrGrid::DLRGrid) = (length(dlrGrid.ω),) # following the Julia convention: size(vector) returns (length(vector),)
+# Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω) # following the Julia convention: size(vector) returns (length(vector),)
 Base.length(dlrGrid::DLRGrid) = length(dlrGrid.ω)
 rank(dlrGrid::DLRGrid) = length(dlrGrid.ω)
 
-function _load!(dlrGrid::DLRGrid, dlrfile, verbose = false)
+function _load!(dlrGrid::DLRGrid, dlrfile, verbose=false)
 
-    grid = readdlm(dlrfile, comments = true, comment_char = '#')
+    grid = readdlm(dlrfile, comments=true, comment_char='#')
     # println("reading $filename")
 
     β = dlrGrid.β
@@ -216,7 +224,7 @@ function _load!(dlrGrid::DLRGrid, dlrfile, verbose = false)
     verbose && println(dlrGrid)
 end
 
-function _build!(dlrGrid::DLRGrid, folder, filename, algorithm, verbose = false)
+function _build!(dlrGrid::DLRGrid, folder, filename, algorithm, verbose=false)
     isFermi = dlrGrid.isFermi
     β = dlrGrid.β
     if algorithm == :discrete || dlrGrid.symmetry == :none
diff --git a/src/operation.jl b/src/operation.jl
index b0eea2c..dfa5b14 100644
--- a/src/operation.jl
+++ b/src/operation.jl
@@ -1,39 +1,140 @@
-function _tensor2matrix(tensor, axis)
+# function _tensor2matrix(tensor::AbstractArray{T,N}, axis) where {T,N}
+#     # internal function to move the axis dim to the first index, then reshape the tensor into a matrix
+#     dim = N
+#     n1 = size(tensor)[axis]
+#     partialsize = deleteat!(collect(size(tensor)), axis) # the size of the tensor except the axis-th dimension
+#     n2 = reduce(*, partialsize)
+
+#     if axis == 1 #no need to permutate the axis
+#         return reshape(tensor, (n1, n2)), partialsize
+#     elseif axis == 2 && dim == 2 #for matrix, simply transpose, no copy is created
+#         return transpose(tensor), partialsize
+#     else
+#         permu = [i for i = 1:dim]
+#         permu[1], permu[axis] = axis, 1
+#         partialsize = collect(size(tensor)[permu][2:end])
+#         ntensor = permutedims(tensor, permu) # permutate the axis-th and the 1st dim, a copy of the tensor is created 
+#         # ntensor = nocopy ? PermutedDimsArray(tensor, permu) : permutedims(tensor, permu) # permutate the axis-th and the 1st dim
+#         ntensor = reshape(ntensor, (n1, n2)) # no copy is created
+#         return ntensor, partialsize
+#     end
+# end
+
+# function _matrix2tensor(mat, partialsize, axis)
+#     # internal function to reshape matrix to a tensor, then swap the first index with the axis-th dimension
+#     @assert size(mat)[2] == reduce(*, partialsize) # total number of elements of mat and the tensor must match
+#     tsize = vcat(size(mat)[1], partialsize)
+#     tensor = reshape(mat, Tuple(tsize))
+#     dim = length(partialsize) + 1
+
+#     if axis == 1
+#         return tensor
+#     elseif axis == 2 && dim == 2
+#         return transpose(tensor) #transpose do not create copy
+#     else
+#         permu = [i for i = 1:dim]
+#         permu[1], permu[axis] = axis, 1
+#         return permutedims(tensor, permu) # permutate the axis-th and the 1st dim, a copy of the tensor is created
+#         # ntensor = nocopy ? PermutedDimsArray(tensor, permu) : permutedims(tensor, permu) # permutate the axis-th and the 1st dim
+#         # return ntensor
+#     end
+# end
+
+#replace one of the tuple elements. See https://discourse.julialang.org/t/computing-tuple-replacements/69581/10
+@generated function _reset_tuple(t::NTuple{N,Int}, x, i) where {N}
+    Expr(:tuple, (:(ifelse($j == i, x, t[$j])) for j in 1:N)...)
+end
+
+# @generated function _remove_tuple(t::NTuple{N,Int}, x, i) where {N}
+#     Expr(:tuple, (:(ifelse($j == i, x, t[$j])) for j in 1:N)...)
+# end
+
+function _matrix_tensor_dot(mat::AbstractMatrix{TC}, tensor::AbstractArray{T,N}, axis::Int) where {T,TC,N}
+    #calculate \sum_j mat[i, j]*tensor[..., j, ...]  where j is the axis-th dimension of tensor
+    @assert 0 < axis <= N
+    _n, _m = size(mat)
+    _size = collect(size(tensor))
+    @assert (_m == _size[axis]) "matrix size $(size(mat)) and tensor size ($(_size)) do not match at axis = $axis"
+    _target_size = _reset_tuple(size(tensor), _n, axis)
+    if axis == 1
+        _r = reduce(*, _size[axis+1:end])
+        _tensor = reshape(tensor, (_m, _r))
+        res = mat * _tensor
+        # _target_size = (_n, _size[2:end]...)::NTuple{N,Int}
+        return reshape(res, _target_size)
+    elseif axis == N
+        # _target_size = (_size[1:end-1]..., _n)::NTuple{N,Int}
+        _l = reduce(*, _size[1:axis-1])
+        _tensor = reshape(tensor, (_l, _m))
+        res = _tensor * transpose(mat)
+        return reshape(res, _target_size)
+    else
+        _l = reduce(*, _size[1:axis-1])
+        _r = reduce(*, _size[axis+1:end])
+        _tensor = reshape(tensor, (_l, _m, _r))
+        res = zeros(promote_type(T, TC), _l, _n, _r)
+        @inbounds for j = 1:_r
+            @inbounds for q = 1:_n
+                @inbounds for k = 1:_m
+                    @inbounds for i = 1:_l
+                        res[i, q, j] += _tensor[i, k, j] * mat[q, k]
+                    end
+                end
+            end
+        end
+        return reshape(res, _target_size)
+    end
+end
+
+function _tensor2matrix(tensor::AbstractVector{T}, ::Val{axis}) where {T,axis}
+    return reshape(tensor, length(tensor), 1), nothing
+end
+
+function _tensor2matrix(tensor::AbstractArray{T,N}, ::Val{axis}) where {T,N,axis}
     # internal function to move the axis dim to the first index, then reshape the tensor into a matrix
-    dim = length(size(tensor))
-    n1 = size(tensor)[axis]
-    partialsize = deleteat!(collect(size(tensor)), axis) # the size of the tensor except the axis-th dimension
+    _size = size(tensor)
+    n1 = _size[axis]
+    partialsize = (_size[1:axis-1]..., _size[axis+1:end]...) # the size of the tensor except the axis-th dimension
     n2 = reduce(*, partialsize)
 
     if axis == 1 #no need to permutate the axis
         return reshape(tensor, (n1, n2)), partialsize
-    elseif axis == 2 && dim == 2 #for matrix, simply transpose, no copy is created
-        return transpose(tensor), partialsize
+    elseif axis == N
+        _tensor = reshape(tensor, (n2, n1))
+        return transpose(_tensor), partialsize  #transpose do not create copy
     else
-        permu = [i for i = 1:dim]
-        permu[1], permu[axis] = axis, 1
-        partialsize = collect(size(tensor)[permu][2:end])
-        ntensor = permutedims(tensor, permu) # permutate the axis-th and the 1st dim, a copy of the tensor is created 
+        permu = (axis, 2:axis-1..., 1, axis+1:N...)
+        _tensor = permutedims(tensor, permu) # permutate the axis-th and the 1st dim, a copy of the tensor is created 
         # ntensor = nocopy ? PermutedDimsArray(tensor, permu) : permutedims(tensor, permu) # permutate the axis-th and the 1st dim
-        ntensor = reshape(ntensor, (n1, n2)) # no copy is created
+        ntensor = reshape(_tensor, (n1, n2)) # no copy is created
         return ntensor, partialsize
     end
 end
 
-function _matrix2tensor(mat, partialsize, axis)
+function _matrix2tensor(mat::AbstractMatrix{T}, partialsize::Nothing, ::Val{axis}) where {T,axis}
+    return reshape(mat, length(mat))
+end
+
+function _matrix2tensor(mat::AbstractMatrix{T}, partialsize::NTuple{dim,Int}, ::Val{axis}) where {T,dim,axis}
     # internal function to reshape matrix to a tensor, then swap the first index with the axis-th dimension
-    @assert size(mat)[2] == reduce(*, partialsize) # total number of elements of mat and the tensor must match
-    tsize = vcat(size(mat)[1], partialsize)
-    tensor = reshape(mat, Tuple(tsize))
-    dim = length(partialsize) + 1
+    n1, n2 = size(mat)
+    @assert n2 == reduce(*, partialsize) # total number of elements of mat and the tensor must match
+    # tsize = vcat(n1, partialsize)
+    N = dim + 1
 
     if axis == 1
+        tsize = (n1, partialsize...)
+        tensor = reshape(mat, tsize)
         return tensor
-    elseif axis == 2 && dim == 2
-        return transpose(tensor) #transpose do not create copy
+    elseif axis == N # mat must be transposed to (n2, n1)
+        _tensor = transpose(mat) #transpose do not create copy
+        return reshape(_tensor, (partialsize..., n1))
     else
-        permu = [i for i = 1:dim]
-        permu[1], permu[axis] = axis, 1
+        # permu = [i for i = 1:dim]
+        # permu[1], permu[axis] = axis, 1
+        permu = (axis, 2:axis-1..., 1, axis+1:N...)
+        tsize = (n1, partialsize...)
+        tensor = reshape(mat, tsize)
         return permutedims(tensor, permu) # permutate the axis-th and the 1st dim, a copy of the tensor is created
         # ntensor = nocopy ? PermutedDimsArray(tensor, permu) : permutedims(tensor, permu) # permutate the axis-th and the 1st dim
         # return ntensor
@@ -121,7 +222,7 @@ end
 function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
     Nτ, Nω = size(kernel)
     @assert size(Gτ)[1] == Nτ
-    if isnothing(sumrule) == false
+    if isnothing(sumrule) == false #require sumrule
         @assert dlrGrid.symmetry == :none && dlrGrid.isFermi "only unsymmetrized ferminoic sum rule has been implemented!"
         # println(size(Gτ))
         M = Int(floor(dlrGrid.size / 2))
@@ -149,7 +250,8 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
         w = 1.0 ./ (error .+ 1e-16)
 
         for i = 1:Nτ
-            w[i, :] /= sum(w[i, :]) / length(w[i, :])
+            wview = view(w, i, :)
+            w[i, :] /= sum(wview) / length(wview)
         end
         B = w .* kernel
         C = w .* Gτ
@@ -164,11 +266,13 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
             Gτ[i, :] .+= kernel_m0[i] * sumrule
         end
         #add back the coeff that are fixed by the sum rule
-        coeffmore = sumrule' .- sum(coeff, dims = 1)
+        coeffmore = sumrule' .- sum(coeff, dims=1)
         cnew = zeros(eltype(coeff), size(coeff)[1] + 1, size(coeff)[2])
-        cnew[1:M-1, :] = coeff[1:M-1, :]
-        cnew[M+1:end, :] = coeff[M:end, :]
-        cnew[M, :] = coeffmore
+        cnew[1:M-1, :] .= coeff[1:M-1, :]
+        cnew[M+1:end, :] .= coeff[M:end, :]
+        # println(size(coeffmore), ", ", size(cnew))
+        # println(coeffmore)
+        cnew[M, :] = coeffmore #broadcast cnew[M, :] .= coeffmore doesn't work for Julia 1.6
         return cnew
     else
         return coeff
@@ -189,40 +293,34 @@ function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing,
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing, axis = 1, sumrule = nothing, verbose = true)
+function tau2dlr(dlrGrid::DLRGrid{T,S}, green::AbstractArray{TC,N}, τGrid=dlrGrid.τ; error=nothing, axis=1, sumrule=nothing, verbose=true) where {T,S,TC,N}
     @assert length(size(green)) >= axis "dimension of the Green's function should be larger than axis!"
     @assert size(green)[axis] == length(τGrid)
     ωGrid = dlrGrid.ω
 
-    typ = promote_type(eltype(dlrGrid.kernel_τ), eltype(green))
-
-    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol = 1e-14)
+    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol=10 * eps(T))
+        if length(dlrGrid.kernel_τ) == 1
+            dlrGrid.kernel_τ = Spectral.kernelT(T, Val(dlrGrid.isFermi), Val(S), τGrid, ωGrid, dlrGrid.β, true)
+        end
         kernel = dlrGrid.kernel_τ
     else
-        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), τGrid, ωGrid, dlrGrid.β, true; type = typ)
+        kernel = Spectral.kernelT(T, Val(dlrGrid.isFermi), Val(S), τGrid, ωGrid, dlrGrid.β, true)
     end
 
-    if typ != eltype(kernel)
-        kernel = convert.(typ, kernel)
-    end
-    if typ != eltype(green)
-        green = convert.(typ, green)
-    end
-
-    g, partialsize = _tensor2matrix(green, axis)
+    g, partialsize = _tensor2matrix(green, Val(axis))
 
     if isnothing(sumrule) == false
         # if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha
         #     sumrule = sumrule ./ 2.0
         # end
-        if isempty(partialsize) == false
+        if isnothing(partialsize) == false
             sumrule = reshape(sumrule, size(g)[2])
         end
     end
 
     if isnothing(error) == false
         @assert size(error) == size(green)
-        error, psize = _tensor2matrix(error, axis)
+        error, partialsize = _tensor2matrix(error, Val(axis))
     end
 
     coeff = _weightedLeastSqureFit(dlrGrid, g, error, kernel, sumrule)
@@ -233,13 +331,13 @@ function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing,
 
     if isnothing(sumrule) == false
         #check how exact is the sum rule
-        coeffsum = sum(coeff, dims = 1) .- sumrule
+        coeffsum = sum(coeff, dims=1) .- sumrule
         if verbose && all(x -> abs(x) < 1000 * dlrGrid.rtol * max(maximum(abs.(green)), 1.0), coeffsum) == false
             @warn("Sumrule error $(maximum(abs.(coeffsum))) is larger than the DLRGrid error threshold.")
         end
     end
 
-    return _matrix2tensor(coeff, partialsize, axis)
+    return _matrix2tensor(coeff, partialsize, Val(axis))
 end
 
 """
@@ -254,24 +352,28 @@ function dlr2tau(dlrGrid::DLRGrid, dlrcoeff, τGrid = dlrGrid.τ; axis = 1, verb
 - `axis`     : imaginary-time axis in the data `dlrcoeff`
 - `verbose`  : true to print warning information
 """
-function dlr2tau(dlrGrid::DLRGrid, dlrcoeff, τGrid = dlrGrid.τ; axis = 1, verbose = true)
+function dlr2tau(dlrGrid::DLRGrid{T,S}, dlrcoeff::AbstractArray{TC,N}, τGrid=dlrGrid.τ; axis=1, verbose=true) where {T,S,TC,N}
     @assert length(size(dlrcoeff)) >= axis "dimension of the dlr coefficients should be larger than axis!"
-    @assert size(dlrcoeff)[axis] == size(dlrGrid)
+    @assert size(dlrcoeff)[axis] == length(dlrGrid)
 
     β = dlrGrid.β
     ωGrid = dlrGrid.ω
 
-    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol = 1e-14)
+    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol=10 * eps(T))
+        if length(dlrGrid.kernel_τ) == 1
+            dlrGrid.kernel_τ = Spectral.kernelT(T, Val(dlrGrid.isFermi), Val(S), τGrid, ωGrid, dlrGrid.β, true)
+        end
         kernel = dlrGrid.kernel_τ
     else
-        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), τGrid, ωGrid, dlrGrid.β, true)
+        kernel = Spectral.kernelT(T, Val(dlrGrid.isFermi), Val(S), τGrid, ωGrid, dlrGrid.β, true)
     end
 
-    coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
+    # coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
 
-    G = kernel * coeff # tensor dot product: \sum_i kernel[..., i]*coeff[i, ...]
+    # G = kernel * coeff # tensor dot product: \sum_i kernel[..., i]*coeff[i, ...]
 
-    return _matrix2tensor(G, partialsize, axis)
+    # return _matrix2tensor(G, partialsize, axis)
+    return _matrix_tensor_dot(kernel, dlrcoeff, axis)
 end
 
 """
@@ -288,42 +390,56 @@ function matfreq2dlr(dlrGrid::DLRGrid, green, nGrid = dlrGrid.n; error = nothing
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function matfreq2dlr(dlrGrid::DLRGrid, green, nGrid = dlrGrid.n; error = nothing, axis = 1, sumrule = nothing, verbose = true)
+function matfreq2dlr(dlrGrid::DLRGrid{T,S}, green::AbstractArray{TC,N}, nGrid=dlrGrid.n; error=nothing, axis=1, sumrule=nothing, verbose=true) where {T,S,TC,N}
     @assert length(size(green)) >= axis "dimension of the Green's function should be larger than axis!"
     @assert size(green)[axis] == length(nGrid)
     @assert eltype(nGrid) <: Integer
     ωGrid = dlrGrid.ω
 
-    typ = promote_type(eltype(dlrGrid.kernel_n), eltype(green))
+    # typ = promote_type(eltype(dlrGrid.kernel_n), eltype(green))
 
-    if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol = 1e-14)
-        kernel = dlrGrid.kernel_n
+    if (S == :ph && dlrGrid.isFermi == false) || (S == :pha && dlrGrid.isFermi == true)
+        if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol=10 * eps(T))
+            if length(dlrGrid.kernel_n) == 1
+                dlrGrid.kernel_n = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+            end
+            kernel = dlrGrid.kernel_n
+        else
+            kernel = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+        end
     else
-        kernel = Spectral.kernelΩ(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), nGrid, ωGrid, dlrGrid.β, true; type = typ)
+        if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol=10 * eps(T))
+            if length(dlrGrid.kernel_n) == 1
+                dlrGrid.kernel_nc = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+            end
+            kernel = dlrGrid.kernel_nc
+        else
+            kernel = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+        end
     end
 
-    if typ != eltype(green)
-        green = convert.(typ, green)
-    end
+    # if typ != eltype(green)
+    #     green = convert.(typ, green)
+    # end
 
-    if typ != eltype(kernel)
-        dlrGrid.kernel_n = convert.(typ, dlrGrid.kernel_n)
-    end
+    # if typ != eltype(kernel)
+    #     dlrGrid.kernel_n = convert.(typ, dlrGrid.kernel_n)
+    # end
 
-    g, partialsize = _tensor2matrix(green, axis)
+    g, partialsize = _tensor2matrix(green, Val(axis))
 
     if isnothing(sumrule) == false
         # if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha
         #     sumrule = sumrule ./ 2.0
         # end
-        if isempty(partialsize) == false
+        if isnothing(partialsize) == false
             sumrule = reshape(sumrule, size(g)[2])
         end
     end
 
     if isnothing(error) == false
         @assert size(error) == size(green)
-        error, psize = _tensor2matrix(error, axis)
+        error, partialsize = _tensor2matrix(error, Val(axis))
     end
     coeff = _weightedLeastSqureFit(dlrGrid, g, error, kernel, sumrule)
     if verbose && all(x -> abs(x) < 1e16, coeff) == false
@@ -332,12 +448,12 @@ function matfreq2dlr(dlrGrid::DLRGrid, green, nGrid = dlrGrid.n; error = nothing
 
     if isnothing(sumrule) == false
         #check how exact is the sum rule
-        coeffsum = sum(coeff, dims = 1) .- sumrule
+        coeffsum = sum(coeff, dims=1) .- sumrule
         if verbose && all(x -> abs(x) < 1000 * dlrGrid.rtol * max(maximum(abs.(green)), 1.0), coeffsum) == false
             @warn("Sumrule error $(maximum(abs.(coeffsum))) is larger than the DLRGrid error threshold.")
         end
     end
-    return _matrix2tensor(coeff, partialsize, axis)
+    return _matrix2tensor(coeff, partialsize, Val(axis))
 end
 
 """
@@ -352,23 +468,39 @@ function dlr2matfreq(dlrGrid::DLRGrid, dlrcoeff, nGrid = dlrGrid.n; axis = 1, ve
 - `axis`     : Matsubara-frequency axis in the data `dlrcoeff`
 - `verbose`  : true to print warning information
 """
-function dlr2matfreq(dlrGrid::DLRGrid, dlrcoeff, nGrid = dlrGrid.n; axis = 1, verbose = true)
+function dlr2matfreq(dlrGrid::DLRGrid{T,S}, dlrcoeff::AbstractArray{TC,N}, nGrid=dlrGrid.n; axis=1, verbose=true) where {T,S,TC,N}
     @assert length(size(dlrcoeff)) >= axis "dimension of the dlr coefficients should be larger than axis!"
-    @assert size(dlrcoeff)[axis] == size(dlrGrid)
+    @assert size(dlrcoeff)[axis] == length(dlrGrid)
     @assert eltype(nGrid) <: Integer
     ωGrid = dlrGrid.ω
 
-    if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol = 1e-14)
-        kernel = dlrGrid.kernel_n
+    if (S == :ph && dlrGrid.isFermi == false) || (S == :pha && dlrGrid.isFermi == true)
+        if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol=10 * eps(T))
+            if length(dlrGrid.kernel_n) == 1
+                dlrGrid.kernel_n = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+            end
+            kernel = dlrGrid.kernel_n
+        else
+            kernel = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+        end
     else
-        kernel = Spectral.kernelΩ(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), nGrid, ωGrid, dlrGrid.β, true)
+        if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol=10 * eps(T))
+            if length(dlrGrid.kernel_n) == 1
+                dlrGrid.kernel_nc = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+            end
+            kernel = dlrGrid.kernel_nc
+        else
+            kernel = Spectral.kernelΩ(T, Val(dlrGrid.isFermi), Val(S), nGrid, ωGrid, dlrGrid.β, true)
+        end
     end
 
-    coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
+    # coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
 
-    G = kernel * coeff # tensor dot product: \sum_i kernel[..., i]*coeff[i, ...]
+    # G = kernel * coeff # tensor dot product: \sum_i kernel[..., i]*coeff[i, ...]
 
-    return _matrix2tensor(G, partialsize, axis)
+    # return _matrix2tensor(G, partialsize, axis)
+
+    return _matrix_tensor_dot(kernel, dlrcoeff, axis)
 end
 
 """
@@ -386,9 +518,10 @@ function tau2matfreq(dlrGrid, green, nNewGrid = dlrGrid.n, τGrid = dlrGrid.τ;
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function tau2matfreq(dlrGrid, green, nNewGrid = dlrGrid.n, τGrid = dlrGrid.τ; error = nothing, axis = 1, sumrule = nothing, verbose = true)
-    coeff = tau2dlr(dlrGrid, green, τGrid; error = error, axis = axis, sumrule = sumrule, verbose = verbose)
-    return dlr2matfreq(dlrGrid, coeff, nNewGrid, axis = axis, verbose = verbose)
+function tau2matfreq(dlrGrid::DLRGrid{T,S}, green::AbstractArray{TC,N}, nNewGrid::AbstractVector{Int}=dlrGrid.n, τGrid=dlrGrid.τ;
+    error=nothing, axis=1, sumrule=nothing, verbose=true) where {T,S,TC,N}
+    coeff = tau2dlr(dlrGrid, green, τGrid; error=error, axis=axis, sumrule=sumrule, verbose=verbose)
+    return dlr2matfreq(dlrGrid, coeff, nNewGrid, axis=axis, verbose=verbose)
 end
 
 """
@@ -406,9 +539,9 @@ function matfreq2tau(dlrGrid, green, τNewGrid = dlrGrid.τ, nGrid = dlrGrid.n;
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function matfreq2tau(dlrGrid, green, τNewGrid = dlrGrid.τ, nGrid = dlrGrid.n; error = nothing, axis = 1, sumrule = nothing, verbose = true)
-    coeff = matfreq2dlr(dlrGrid, green, nGrid; error = error, axis = axis, sumrule = sumrule, verbose = verbose)
-    return dlr2tau(dlrGrid, coeff, τNewGrid, axis = axis, verbose = verbose)
+function matfreq2tau(dlrGrid, green, τNewGrid=dlrGrid.τ, nGrid=dlrGrid.n; error=nothing, axis=1, sumrule=nothing, verbose=true)
+    coeff = matfreq2dlr(dlrGrid, green, nGrid; error=error, axis=axis, sumrule=sumrule, verbose=verbose)
+    return dlr2tau(dlrGrid, coeff, τNewGrid, axis=axis, verbose=verbose)
 end
 
 """
@@ -426,9 +559,9 @@ function tau2tau(dlrGrid, green, τNewGrid, τGrid = dlrGrid.τ; error = nothing
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function tau2tau(dlrGrid, green, τNewGrid, τGrid = dlrGrid.τ; error = nothing, axis = 1, sumrule = nothing, verbose = true)
-    coeff = tau2dlr(dlrGrid, green, τGrid; error = error, axis = axis, sumrule = sumrule, verbose = verbose)
-    return dlr2tau(dlrGrid, coeff, τNewGrid, axis = axis, verbose = verbose)
+function tau2tau(dlrGrid, green, τNewGrid, τGrid=dlrGrid.τ; error=nothing, axis=1, sumrule=nothing, verbose=true)
+    coeff = tau2dlr(dlrGrid, green, τGrid; error=error, axis=axis, sumrule=sumrule, verbose=verbose)
+    return dlr2tau(dlrGrid, coeff, τNewGrid, axis=axis, verbose=verbose)
 end
 
 """
@@ -446,9 +579,9 @@ function matfreq2matfreq(dlrGrid, green, nNewGrid, nGrid = dlrGrid.n; error = no
 - `sumrule`  : enforce the sum rule 
 - `verbose`  : true to print warning information
 """
-function matfreq2matfreq(dlrGrid, green, nNewGrid, nGrid = dlrGrid.n; error = nothing, axis = 1, sumrule = nothing, verbose = true)
-    coeff = matfreq2dlr(dlrGrid, green, nGrid; error = error, axis = axis, sumrule = sumrule, verbose = verbose)
-    return dlr2matfreq(dlrGrid, coeff, nNewGrid, axis = axis, verbose = verbose)
+function matfreq2matfreq(dlrGrid, green, nNewGrid, nGrid=dlrGrid.n; error=nothing, axis=1, sumrule=nothing, verbose=true)
+    coeff = matfreq2dlr(dlrGrid, green, nGrid; error=error, axis=axis, sumrule=sumrule, verbose=verbose)
+    return dlr2matfreq(dlrGrid, coeff, nNewGrid, axis=axis, verbose=verbose)
 end
 
 # function convolution(dlrGrid, green1, green2; axis = 1)
diff --git a/src/sample/sample.jl b/src/sample/sample.jl
index 75176e2..d9b8d23 100644
--- a/src/sample/sample.jl
+++ b/src/sample/sample.jl
@@ -1,5 +1,6 @@
 module Sample
 using FastGaussQuadrature
+import ..DLRGrid
 using ..Spectral
 """
     SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol = :none; rtol = nothing, degree = 24, regularized::Bool = true)
@@ -20,7 +21,7 @@ Return the function on Grid and the systematic error.
 - `degree`: polynomial degree for integral
 - `regularized`: use regularized bosonic kernel if symmetry = :none
 """
-function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol = :none; rtol = nothing, degree = 24, regularized::Bool = true)
+function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol=:none; rtol=nothing, degree=24, regularized::Bool=true, dtype=Float64)
     # calculate Green's function defined by the spectral density
     # S(ω) = sqrt(1 - (ω / Euv)^2) / Euv # semicircle -1<ω<1
     if type == :τ
@@ -30,7 +31,19 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     else
         error("$type is not implemented!")
     end
-
+    return SemiCircle(Val(isFermi), Val(IsMatFreq), Val(symmetry), Euv, β, Grid; rtol=rtol, degree=degree, regularized=true, dtype=Float64)
+end
+function SemiCircle(dlr::DLRGrid{T,S}, type::Symbol, Grid=dlrGrid(dlr, type); degree=24, regularized::Bool=true, dtype=Float64) where {T,S}
+    if type == :τ
+        IsMatFreq = false
+    elseif type == :n
+        IsMatFreq = true
+    else
+        error("$type is not implemented!")
+    end
+    return SemiCircle(Val(dlr.isFermi), Val(IsMatFreq), Val(S), dlr.Euv, dlr.β, Grid; rtol=dlr.rtol, degree=degree, regularized=regularized, dtype=T)
+end
+function SemiCircle(::Val{isFermi}, ::Val{IsMatFreq}, ::Val{symmetry}, Euv, β, Grid; rtol=nothing, degree=24, regularized::Bool=true, dtype=Float64) where {isFermi,IsMatFreq,symmetry}
     ##### Panels endpoints for composite quadrature rule ###
     npo = Int(ceil(log(β * Euv) / log(2.0)))
     pbp = zeros(Float64, 2npo + 1)
@@ -48,19 +61,16 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     # else
     # end
 
-    g1 = _Green(Val(IsMatFreq), Euv, β, isFermi, Grid, symmetry, degree, pbp, npo, regularized)
+    g1 = _Green(dtype, Val(IsMatFreq), dtype(Euv), dtype(β), Val(isFermi), Grid, Val(symmetry), degree, pbp, npo, regularized)
     # g1 = _Green(IsMatFreq, Euv, β, isFermi, Grid, symmetry, degree, pbp, npo, regularized)
 
     if isnothing(rtol) == false
-        g2 = _Green(Val(IsMatFreq), Euv, β, isFermi, Grid, symmetry, degree * 2, pbp, npo, regularized)
+        g2 = _Green(dtype, Val(IsMatFreq), dtype(Euv), dtype(β), Val(isFermi), Grid, Val(symmetry), degree * 2, pbp, npo, regularized)
         err = abs.(g1 - g2)
         @assert maximum(err) < rtol "Systematic error $(maximum(err)) is larger than $rtol, increase degree for the integral!"
     end
     return g1
 end
-function SemiCircle(dlr, type::Symbol, Grid = dlrGrid(dlr, type); degree = 24, regularized::Bool = true)
-    return SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, Grid, type, dlr.symmetry; rtol = dlr.rtol, degree = degree, regularized = regularized)
-end
 
 # @inline function kernelΩ(isFermi, symmetry, regularized::Bool = false) where {T<:AbstractFloat,isFermi,symmetry}
 #     if symmetry == :none
@@ -95,51 +105,88 @@ end
 # function getG(::Val{false}, Grid)
 #     return zeros(Float64, length(Grid))
 # end
+function _Green(::Type{T}, ::Val{true}, Euv::T, β::T, ::Val{isFermi}, Grid::AbstractVector, ::Val{sym}, n, pbp, npo, regularized) where {T<:Real,IsMatFreq,isFermi,sym}
+    #n: polynomial order
+    xl, wl = gausslegendre(n)
+    xj, wj = gaussjacobi(n, 1 / 2, T(0))
+
+    xl, wl = T.(xl), T.(wl)
+    xj, wj = T.(xj), T.(wj)
 
-function _Green(::Val{IsMatFreq}, Euv, β, isFermi, Grid, symmetry, n, pbp, npo, regularized) where {IsMatFreq}
+    G = zeros(Complex{T}, length(Grid))
+    # G = getG(isMatFreq, Grid)
+    for (τi, τ) in enumerate(Grid)
+        for ii = 2:2npo-1
+            a, b = pbp[ii], pbp[ii+1]
+            for jj = 1:n
+                x = (a + b) / 2 + (b - a) / 2 * xl[jj]
+                if x < T(0) && (sym == :ph || sym == :pha)
+                    #spectral density is defined for positivie frequency only for correlation functions
+                    continue
+                end
+                ker = Spectral.kernelΩ(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
+                G[τi] += (b - a) / 2 * wl[jj] * ker * sqrt(1 - x^2)
+            end
+        end
+
+        a, b = T(1) / 2, T(1)
+        for jj = 1:n
+            x = (a + b) / 2 + (b - a) / 2 * xj[jj]
+            ker = Spectral.kernelΩ(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
+            G[τi] += ((b - a) / 2)^T(1.5) * wj[jj] * ker * sqrt(1 + x)
+        end
+
+        if sym != :ph && sym != :pha
+            #spectral density is defined for positivie frequency only for correlation functions
+            a, b = -T(1), -T(1) / 2
+            for jj = 1:n
+                x = (a + b) / 2 + (b - a) / 2 * (-xj[n-jj+1])
+                ker = Spectral.kernelΩ(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
+                G[τi] += ((b - a) / 2)^T(1.5) * wj[n-jj+1] * ker * sqrt(1 - x)
+            end
+        end
+    end
+    return G
+end
+
+function _Green(::Type{T}, ::Val{false}, Euv::T, β::T, ::Val{isFermi}, Grid::AbstractVector, ::Val{sym}, n, pbp, npo, regularized) where {T<:Real,IsMatFreq,isFermi,sym}
     #n: polynomial order
     xl, wl = gausslegendre(n)
-    xj, wj = gaussjacobi(n, 1 / 2, 0.0)
-    # println(IsMatFreq)
-    type = Val(isFermi)
-    sym = Val(symmetry)
+    xj, wj = gaussjacobi(n, 1 / 2, T(0))
+
+    xl, wl = T.(xl), T.(wl)
+    xj, wj = T.(xj), T.(wj)
 
-    G = IsMatFreq ? zeros(ComplexF64, length(Grid)) : zeros(Float64, length(Grid))
+    G = zeros(T, length(Grid))
     # G = getG(isMatFreq, Grid)
     for (τi, τ) in enumerate(Grid)
         for ii = 2:2npo-1
             a, b = pbp[ii], pbp[ii+1]
             for jj = 1:n
                 x = (a + b) / 2 + (b - a) / 2 * xl[jj]
-                if (symmetry == :ph || symmetry == :pha) && x < 0.0
+                if x < T(0) && (sym == :ph || sym == :pha)
                     #spectral density is defined for positivie frequency only for correlation functions
                     continue
                 end
-                ker = IsMatFreq ?
-                      Spectral.kernelΩ(type, sym, τ, Euv * x, β, regularized) :
-                      Spectral.kernelT(type, sym, τ, Euv * x, β, regularized)
+                ker = Spectral.kernelT(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
                 G[τi] += (b - a) / 2 * wl[jj] * ker * sqrt(1 - x^2)
             end
         end
 
-        a, b = 1.0 / 2, 1.0
+        a, b = T(1) / 2, T(1)
         for jj = 1:n
             x = (a + b) / 2 + (b - a) / 2 * xj[jj]
-            ker = IsMatFreq ?
-                  Spectral.kernelΩ(type, sym, τ, Euv * x, β, regularized) :
-                  Spectral.kernelT(type, sym, τ, Euv * x, β, regularized)
-            G[τi] += ((b - a) / 2)^1.5 * wj[jj] * ker * sqrt(1 + x)
+            ker = Spectral.kernelT(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
+            G[τi] += ((b - a) / 2)^T(1.5) * wj[jj] * ker * sqrt(1 + x)
         end
 
-        if symmetry != :ph && symmetry != :pha
+        if sym != :ph && sym != :pha
             #spectral density is defined for positivie frequency only for correlation functions
-            a, b = -1.0, -1.0 / 2
+            a, b = -T(1), -T(1) / 2
             for jj = 1:n
                 x = (a + b) / 2 + (b - a) / 2 * (-xj[n-jj+1])
-                ker = IsMatFreq ?
-                      Spectral.kernelΩ(type, sym, τ, Euv * x, β, regularized) :
-                      Spectral.kernelT(type, sym, τ, Euv * x, β, regularized)
-                G[τi] += ((b - a) / 2)^1.5 * wj[n-jj+1] * ker * sqrt(1 - x)
+                ker = Spectral.kernelT(Val(isFermi), Val(sym), τ, Euv * x, β, regularized)
+                G[τi] += ((b - a) / 2)^T(1.5) * wj[n-jj+1] * ker * sqrt(1 - x)
             end
         end
     end
@@ -163,7 +210,7 @@ Return the function on Grid and the systematic error.
 - `poles`: a list of frequencies for the delta functions
 - `regularized`: use regularized bosonic kernel if symmetry = :none
 """
-function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbol = :none; regularized::Bool = true)
+function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbol=:none; regularized::Bool=true)
     # poles = [-Euv, -0.2 * Euv, 0.0, 0.8 * Euv, Euv]
     # poles=[0.8Euv, 1.0Euv]
     # poles = [0.0]
@@ -193,8 +240,8 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
     end
     return g
 end
-function MultiPole(dlr, type::Symbol, poles, Grid = dlrGrid(dlr, type); regularized::Bool = true)
-    return MultiPole(dlr.β, dlr.isFermi, Grid, type, poles, dlr.symmetry; regularized = regularized)
+function MultiPole(dlr, type::Symbol, poles, Grid=dlrGrid(dlr, type); regularized::Bool=true)
+    return MultiPole(dlr.β, dlr.isFermi, Grid, type, poles, dlr.symmetry; regularized=regularized)
 end
 
 end
\ No newline at end of file
diff --git a/src/spectral.jl b/src/spectral.jl
index 6e8252e..ddd4dbe 100644
--- a/src/spectral.jl
+++ b/src/spectral.jl
@@ -18,7 +18,7 @@ Compute the imaginary-time kernel of different type.
 - `β`: the inverse temperature 
 - `regularized`: use regularized kernel or not
 """
-@inline function kernelT(::Val{isFermi}, ::Val{symmetry}, τ::T, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat,isFermi,symmetry}
+@inline function kernelT(::Val{isFermi}, ::Val{symmetry}, τ::T, ω::T, β::T, regularized::Bool=false) where {T<:AbstractFloat,isFermi,symmetry}
     if symmetry == :none
         if regularized
             return isFermi ? kernelFermiT(τ, ω, β) : kernelBoseT_regular(τ, ω, β)
@@ -38,11 +38,11 @@ end
 
 Compute kernel with given τ and ω grids.
 """
-function kernelT(isFermi, symmetry, τGrid::AbstractVector{T}, ωGrid::AbstractVector{T}, β::T, regularized::Bool = false; type = T) where {T<:AbstractFloat}
-    kernel = zeros(type, (length(τGrid), length(ωGrid)))
+function kernelT(::Type{T}, isFermi, symmetry, τGrid::AbstractVector{T}, ωGrid::AbstractVector{T}, β::T, regularized::Bool=false) where {T<:AbstractFloat}
+    kernel = zeros(T, (length(τGrid), length(ωGrid)))
     for (τi, τ) in enumerate(τGrid)
         for (ωi, ω) in enumerate(ωGrid)
-            kernel[τi, ωi] = kernelT(isFermi, symmetry, τ, ω, β, regularized)
+            kernel[τi, ωi] = kernelT(isFermi, symmetry, T(τ), T(ω), T(β), regularized)
         end
     end
     return kernel
@@ -164,7 +164,7 @@ K(τ) = e^{-ω|τ|}+e^{-ω(β-|τ|)}
 - `ω`: frequency, ω>=0
 - `β`: the inverse temperature 
 """
-@inline function kernelFermiT_PH(τ::T, ω::T, β = T(1)) where {T<:AbstractFloat}
+@inline function kernelFermiT_PH(τ::T, ω::T, β=T(1)) where {T<:AbstractFloat}
     (-β < τ <= β) || error("τ must be (0, β]")
     (ω >= 0) || error("ω must be >=0")
     τ = abs(τ)
@@ -184,7 +184,7 @@ K(τ) = e^{-ω|τ|}+e^{-ω(β-|τ|)}
 - `ω`: frequency, ω>=0
 - `β`: the inverse temperature 
 """
-@inline function kernelBoseT_PH(τ::T, ω::T, β = T(1)) where {T<:AbstractFloat}
+@inline function kernelBoseT_PH(τ::T, ω::T, β=T(1)) where {T<:AbstractFloat}
     (-β < τ <= β) || error("τ must be (0, β]")
     (ω >= 0) || error("ω must be >=0")
     τ = abs(τ)
@@ -245,7 +245,7 @@ Compute the imaginary-time kernel of different type. Assume ``k_B T/\\hbar=1``
 - `β`: the inverse temperature 
 - `regularized`: use regularized kernel or not
 """
-@inline function kernelΩ(::Val{isFermi}, ::Val{symmetry}, n::Int, ω::T, β::T, regularized::Bool = false) where {T<:AbstractFloat,isFermi,symmetry}
+@inline function kernelΩ(::Val{isFermi}, ::Val{symmetry}, n::Int, ω::T, β::T, regularized::Bool=false) where {T<:AbstractFloat,isFermi,symmetry}
     if symmetry == :none
         if regularized
             return isFermi ? kernelFermiΩ(n, ω, β) : kernelBoseΩ_regular(n, ω, β)
@@ -266,12 +266,17 @@ end
 
 Compute kernel matrix with given ωn (integer!) and ω grids.
 """
-function kernelΩ(isFermi, symmetry, nGrid::Vector{Int}, ωGrid::Vector{T}, β::T, regularized::Bool = false; type = Complex{T}) where {T<:AbstractFloat}
+function kernelΩ(::Type{T}, ::Val{isFermi}, ::Val{symmetry}, nGrid::AbstractVector{Int}, ωGrid::AbstractVector{T}, β::T, regularized::Bool=false) where {T<:AbstractFloat,isFermi,symmetry}
     # println(type)
-    kernel = zeros(type, (length(nGrid), length(ωGrid)))
+    if (symmetry == :none) || (symmetry == :ph && isFermi == true) || (symmetry == :pha && isFermi == false)
+        kernel = zeros(Complex{T}, (length(nGrid), length(ωGrid)))
+    else
+        kernel = zeros(T, (length(nGrid), length(ωGrid)))
+    end
+    # println(symmetry, ", ", isFermi)
     for (ni, n) in enumerate(nGrid)
         for (ωi, ω) in enumerate(ωGrid)
-            kernel[ni, ωi] = kernelΩ(isFermi, symmetry, n, ω, β, regularized)
+            kernel[ni, ωi] = kernelΩ(Val(isFermi), Val(symmetry), n, T(ω), T(β), regularized)
         end
     end
     return kernel
diff --git a/test/dlr.jl b/test/dlr.jl
index db1d360..fe1a1d3 100644
--- a/test/dlr.jl
+++ b/test/dlr.jl
@@ -3,33 +3,78 @@ using FastGaussQuadrature, Printf
 rtol(x, y) = maximum(abs.(x - y)) / maximum(abs.(x))
 
 # SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, grid, type, dlr.symmetry, rtol = dlr.rtol, degree = 24, regularized = true)
-SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree = 24, regularized = true)
+SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=24, regularized=true)
 
 function MultiPole(dlr, grid, type)
     Euv = dlr.Euv
     poles = [-Euv, -0.2 * Euv, 0.0, 0.8 * Euv, Euv]
     # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true)
-    return Sample.MultiPole(dlr, type, poles, grid; regularized = true)
+    return Sample.MultiPole(dlr, type, poles, grid; regularized=true)
 end
 
-function compare(case, a, b, eps, requiredratio, para = "")
+function compare(case, a, b, eps, requiredratio, para="")
     err = rtol(a, b)
     ratio = isfinite(err) ? Int(round(err / eps)) : 0
     if ratio > 50
-        printstyled("$case, $para err: ", color = :white)
-        printstyled("$(round(err, sigdigits=3)) = $ratio x rtol\n", color = :green)
+        printstyled("$case, $para err: ", color=:white)
+        printstyled("$(round(err, sigdigits=3)) = $ratio x rtol\n", color=:green)
     end
     @test rtol(a, b) .< requiredratio * eps # dlr should represent the Green's function up to accuracy of the order eps
 end
 
-function compare_atol(case, a, b, atol, para = "")
+function compare_atol(case, a, b, atol, para="")
     err = rtol(a, b)
-    err = isfinite(err) ? round(err, sigdigits = 3) : 0
+    err = isfinite(err) ? round(err, sigdigits=3) : 0
     if err > 100 * atol
-        printstyled("$case, $para err: ", color = :white)
-        printstyled("$err = $(Int(round(err/atol))) x atol\n", color = :blue)
+        printstyled("$case, $para err: ", color=:white)
+        printstyled("$err = $(Int(round(err/atol))) x atol\n", color=:blue)
     end
-    @test rtol(a, b) .< 5000 * atol # dlr should represent the Green's function up to accuracy of the order eps
+    @test rtol(a, b) .< 6000 * atol # dlr should represent the Green's function up to accuracy of the order eps
+end
+
+@testset "Operation Utility" begin
+    ########## test _matrix_tensor_dot  #################
+    # middle
+    a = rand(5, 4)
+    b = rand(3, 4, 6)
+    c = Lehmann._matrix_tensor_dot(a, b, 2)
+    @inferred Lehmann._matrix_tensor_dot(a, b, 2)
+
+    @test size(c) == (3, 5, 6)
+
+    _b, psize = Lehmann._tensor2matrix(b, Val(2))
+    _c = a * _b
+    _c = Lehmann._matrix2tensor(_c, psize, Val(2))
+
+    @test c ≈ _c
+
+    # left
+    a = rand(5, 4)
+    b = rand(4, 3, 6)
+    c = Lehmann._matrix_tensor_dot(a, b, 1)
+    @inferred Lehmann._matrix_tensor_dot(a, b, 1)
+
+    @test size(c) == (5, 3, 6)
+
+    _b, psize = Lehmann._tensor2matrix(b, Val(1))
+    _c = a * _b
+    _c = Lehmann._matrix2tensor(_c, psize, Val(1))
+
+    @test c ≈ _c
+
+    # right
+    a = rand(5, 4)
+    b = rand(3, 6, 4)
+    c = Lehmann._matrix_tensor_dot(a, b, 3)
+    @inferred Lehmann._matrix_tensor_dot(a, b, 3)
+
+    @test size(c) == (3, 6, 5)
+
+    _b, psize = Lehmann._tensor2matrix(b, Val(3))
+    _c = a * _b
+    _c = Lehmann._matrix2tensor(_c, psize, Val(3))
+
+    @test c ≈ _c
 end
 
 @testset "Correlator Representation" begin
@@ -53,7 +98,9 @@ end
 
         ########################## imaginary-time to dlr #######################################
         coeff = tau2dlr(dlr, Gdlr)
+        @test size(dlr.kernel_τ) == (length(dlr.τ), length(dlr.ω))
         Gfitted = dlr2tau(dlr, coeff, τSample)
+        @test size(dlr.kernel_τ) == (length(dlr.τ), length(dlr.ω))
         compare("dlr τ → dlr → generic τ $case", Gsample, Gfitted, eps, 100, para)
         # for (ti, t) in enumerate(τSample)
         #     @printf("%32.19g    %32.19g   %32.19g   %32.19g\n", t / β, Gsample[1, ti],  Gfitted[1, ti], Gsample[1, ti] - Gfitted[1, ti])
@@ -101,11 +148,11 @@ end
         atol = eps
         noise = atol * rand(eltype(Gsample), length(Gsample))
         GNoisy = Gsample .+ noise
-        compare_atol("noisy generic τ → dlr → τ $case", tau2tau(dlr, GNoisy, dlr.τ, τSample; error = abs.(noise)), Gdlr, atol, para)
+        compare_atol("noisy generic τ → dlr → τ $case", tau2tau(dlr, GNoisy, dlr.τ, τSample; error=abs.(noise)), Gdlr, atol, para)
 
         noise = atol * rand(eltype(Gnsample), length(Gnsample))
         GnNoisy = Gnsample .+ noise
-        compare_atol("noisy generic iω → dlr → iω $case", matfreq2matfreq(dlr, GnNoisy, dlr.n, nSample, error = abs.(noise)), Gndlr, atol, para)
+        compare_atol("noisy generic iω → dlr → iω $case", matfreq2matfreq(dlr, GnNoisy, dlr.n, nSample, error=abs.(noise)), Gndlr, atol, para)
     end
 
     # the accuracy greatly drops beyond Λ >= 1e8 and rtol<=1e-6
@@ -130,25 +177,53 @@ end
 @testset "Tensor ↔ Matrix Mapping" begin
     a = rand(3)
     acopy = deepcopy(a)
-    b, psize = Lehmann._tensor2matrix(a, 1)
-    anew = Lehmann._matrix2tensor(b, psize, 1)
+    b, psize = Lehmann._tensor2matrix(a, Val(1))
+    anew = Lehmann._matrix2tensor(b, psize, Val(1))
     @test acopy ≈ anew
 
     a = rand(3, 4)
     acopy = deepcopy(a)
     for axis = 1:2
-        b, psize = Lehmann._tensor2matrix(a, axis)
-        anew = Lehmann._matrix2tensor(b, psize, axis)
+        b, psize = Lehmann._tensor2matrix(a, Val(axis))
+        anew = Lehmann._matrix2tensor(b, psize, Val(axis))
         @test acopy ≈ anew
     end
 
     a = rand(3, 4, 5)
     acopy = deepcopy(a)
     for axis = 1:3
-        b, psize = Lehmann._tensor2matrix(a, axis)
-        anew = Lehmann._matrix2tensor(b, psize, axis)
+        b, psize = Lehmann._tensor2matrix(a, Val(axis))
+        anew = Lehmann._matrix2tensor(b, psize, Val(axis))
         @test acopy ≈ anew
     end
+
+    #### _tensor2matrix and _matrix2tensor type stability ############
+    a = rand(5)
+    _a, psize = Lehmann._tensor2matrix(a, Val(1))
+    __a = Lehmann._matrix2tensor(_a, psize, Val(1))
+    @test a ≈ __a
+    @inferred Lehmann._tensor2matrix(a, Val(1))
+    @inferred Lehmann._matrix2tensor(_a, psize, Val(1))
+
+    a = rand(3, 4, 5)
+    _a, psize = Lehmann._tensor2matrix(a, Val(2))
+    __a = Lehmann._matrix2tensor(_a, psize, Val(2))
+    @test a ≈ __a
+    @inferred Lehmann._tensor2matrix(a, Val(2))
+    @inferred Lehmann._matrix2tensor(_a, psize, Val(2))
+
+    _a, psize = Lehmann._tensor2matrix(a, Val(1))
+    __a = Lehmann._matrix2tensor(_a, psize, Val(1))
+    @test a ≈ __a
+    @inferred Lehmann._tensor2matrix(a, Val(1))
+    @inferred Lehmann._matrix2tensor(_a, psize, Val(1))
+
+    _a, psize = Lehmann._tensor2matrix(a, Val(3))
+    __a = Lehmann._matrix2tensor(_a, psize, Val(3))
+    @test a ≈ __a
+    @inferred Lehmann._tensor2matrix(a, Val(3))
+    @inferred Lehmann._matrix2tensor(_a, psize, Val(3))
+
 end
 
 @testset "Tensor DLR" begin
@@ -185,25 +260,49 @@ end
     tensorGn_copy = tensorGn
 
     compare("τ ↔ iω tensor", tau2matfreq(dlr, Gτ_copy), Gn, eps, 1000.0, para)
+    @inferred tau2matfreq(dlr, Gτ_copy)
     @test Gτ ≈ Gτ_copy #make sure there is no side effect on G
     compare("iω ↔ τ tensor", matfreq2tau(dlr, Gn_copy), Gτ, eps, 1000.0, para)
+    @inferred matfreq2tau(dlr, Gn_copy)
     @test Gn ≈ Gn_copy #make sure there is no side effect on G
 
-    compare("τ ↔ iω tensor", tau2matfreq(dlr, Gτ_copy; sumrule = weight), Gn, eps, 1000.0, para)
+    compare("τ ↔ iω tensor", tau2matfreq(dlr, Gτ_copy; sumrule=weight), Gn, eps, 1000.0, para)
     @test Gτ ≈ Gτ_copy #make sure there is no side effect on G
-    compare("iω ↔ τ tensor", matfreq2tau(dlr, Gn_copy; sumrule = weight), Gτ, eps, 1000.0, para)
+    compare("iω ↔ τ tensor", matfreq2tau(dlr, Gn_copy; sumrule=weight), Gτ, eps, 1000.0, para)
     @test Gn ≈ Gn_copy #make sure there is no side effect on G
 
-    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis = 3), tensorGn, eps, 1000.0, para)
+    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis=3), tensorGn, eps, 1000.0, para)
+    # println(typeof(tensorGτ_copy))
+    # @inferred tau2matfreq(dlr, tensorGτ_copy; axis=3)
     @test tensorGτ ≈ tensorGτ_copy #make sure there is no side effect on G
-    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis = 3), tensorGτ, eps, 1000.0, para)
+    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis=3), tensorGτ, eps, 1000.0, para)
+    # @inferred matfreq2tau(dlr, tensorGn_copy; axis=3)
     @test tensorGn ≈ tensorGn_copy #make sure there is no side effect on G
 
-    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis = 3, sumrule = sumrule_τ), tensorGn, eps, 1000.0, para)
+
+    compare("τ ↔ iω tensor", tau2matfreq(dlr, tensorGτ_copy; axis=3, sumrule=sumrule_τ), tensorGn, eps, 1000.0, para)
+    # @inferred tau2matfreq(dlr, tensorGτ_copy; axis=3)
     @test tensorGτ ≈ tensorGτ_copy #make sure there is no side effect on G
-    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis = 3, sumrule = sumrule_n), tensorGτ, eps, 1000.0, para)
+    compare("iω ↔ τ tensor", matfreq2tau(dlr, tensorGn_copy; axis=3, sumrule=sumrule_n), tensorGτ, eps, 1000.0, para)
+    # @inferred matfreq2tau(dlr, tensorGn_copy; axis=3)
     @test tensorGn ≈ tensorGn_copy #make sure there is no side effect on G
 
+    tensor = zeros(128, 128, length(dlr))
+    tau2matfreq(dlr, tensor; axis=3)
+    coeff = tau2dlr(dlr, tensor; axis=3)
+    dlr2matfreq(dlr, coeff; axis=3)
+
+    # @inferred tau2dlr(dlr, tensor; axis=3)
+    # @inferred dlr2matfreq(dlr, coeff; axis=3)
+
+    tensor = zeros(128, 128, length(dlr))
+    println("tau2matfreq timing:")
+    @time tau2matfreq(dlr, tensor; axis=3)
+    println("tau2dlr timing:")
+    @time coeff = tau2dlr(dlr, tensor; axis=3)
+    println("dlr2matfreq timing:")
+    @time dlr2matfreq(dlr, coeff; axis=3)
+
 end
 
 @testset "Least square fitting" begin
@@ -211,7 +310,7 @@ end
     kernel = zeros(4, 2)
     kernel[:, 1] = [1.0, 1.0, 1.0, 1.0]
     kernel[:, 2] = [1.0, 2.0, 3.0, 4.0]
-    dlrGrid = DLRGrid(β = 10.0, isFermi = true)
+    dlrGrid = DLRGrid(β=10.0, isFermi=true)
     coeff = Lehmann._weightedLeastSqureFit(dlrGrid, Gτ, nothing, kernel, nothing)
     @test coeff ≈ [3.5, 1.4]
 end
@@ -232,11 +331,11 @@ end
     Euv, β, rtol = 1.5, 110.0, 4.3e-8
 
     #save dlr to a local file
-    dlr = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild = true, folder = folder, verbose = false)
+    dlr = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild=true, folder=folder, verbose=false)
     file = finddlr(folder, ".dlr")
     @test isnothing(file) == false
     #load dlr from the local file
-    dlr_load = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild = false, folder = folder, verbose = false)
+    dlr_load = Lehmann.DLRGrid(Euv, β, rtol, true; rebuild=false, folder=folder, verbose=false)
     @test dlr_load.τ ≈ dlr.τ
     @test dlr_load.n ≈ dlr.n
     @test dlr_load.ω ≈ dlr.ω
@@ -249,8 +348,8 @@ end
 end
 
 @testset "JLD2 IO" begin
-    dlr = DLRGrid(isFermi = true, beta = 10.0)
-    save("test.jld2", Dict("dlr" => dlr), compress = true)
+    dlr = DLRGrid(isFermi=true, beta=10.0)
+    save("test.jld2", Dict("dlr" => dlr), compress=true)
     dlr_new = load("test.jld2")["dlr"]
     println(dlr_new)
     @test dlr.isFermi == dlr_new.isFermi