diff --git a/example/demoSumRule.jl b/example/demoSumRule.jl
new file mode 100644
index 0000000..24d8e68
--- /dev/null
+++ b/example/demoSumRule.jl
@@ -0,0 +1,36 @@
+using Lehmann
+using Printf
+using Gaston
+β = 40.0
+rtol = 1e-8
+eta = 1e-2
+N = 128
+
+dlr = DLRGrid(β = β, Euv = 1.0, isFermi = true, rtol = rtol, symmetry = :none)
+τgrid = collect(LinRange(0.0, β, N))
+println(τgrid)
+G = Sample.SemiCircle(dlr, :τ, τgrid)
+G += rand(length(G)) / 2.0 * eta
+coeff = tau2dlr(dlr, G, τgrid)
+coeff_sumrule = tau2dlr(dlr, G, τgrid, sumrule = π / 2)
+println("sum rule: ", abs(sum(coeff) - π / 2), " versus ", abs(sum(coeff_sumrule) - π / 2))
+
+
+dlr10 = DLRGrid(β = β, Euv = 100.0, isFermi = true, rtol = rtol, symmetry = :none)
+
+Gtrue = Sample.SemiCircle(dlr, :τ, dlr10.τ)
+Gfit = real.(dlr2tau(dlr, coeff, dlr10.τ))
+Gfit_sumrule = real.(dlr2tau(dlr, coeff_sumrule, dlr10.τ))
+
+@printf("%15s%30s%30s%30s\n", "tau", "true", "no sumrule", "with sumrule")
+for i in 1:dlr10.size
+    @printf("%15.6f%30.15f%30.15e%30.15e\n", dlr10.τ[i], Gtrue[i], abs(Gfit[i] - Gtrue[i]), abs(Gfit_sumrule[i] - Gtrue[i]))
+end
+
+set(term = "qt")
+# p = plot(dlr10.τ, Gtrue, label = "original", axis = "semilogy")
+xrange = (0.0, 40)
+p = plot(dlr10.τ, Gfit - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR no sum rule"), Axes(xrange = xrange))
+plot!(dlr10.τ, Gfit_sumrule - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR with sum rule"))
+display(p)
+readline()
\ No newline at end of file
diff --git a/example/test.jl b/example/test.jl
deleted file mode 100644
index 59e4b28..0000000
--- a/example/test.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-using Lehmann
-using Gaston
-# dlr = DLRGrid(β = 1000.0, isFermi = true, rtol = 1e-6, symmetry = :ph)
-# G = Sample.SemiCircle(dlr, :τ)
-# coeff2 = tau2dlr(dlr, G, sumrule = π / 2)
-# println(coeff2)
-# println("diff: ", sum(coeff2), ", ", sum(coeff2) - π / 2)
-
-dlr = DLRGrid(β = 1000.0, isFermi = true, rtol = 1e-6, symmetry = :none)
-G = Sample.SemiCircle(dlr, :τ)
-coeff2 = tau2dlr(dlr, G, sumrule = π / 2)
-println(coeff2)
-println("diff2: ", sum(coeff2), ", ", sum(coeff2) - π / 2)
-
-# set(term = "qt")
-# p = plot(dlr.τ, G, label = "original", axis = "semilogy")
-# plot!(dlr.τ, real.(dlr2tau(dlr, G)), label = "DLR")
-# display(p)
-# readline()
-Gfit = real.(dlr2tau(dlr, coeff2))
-for i in 1:dlr.size
-    println("$(dlr.τ[i])   $(G[i])  $(Gfit[i])   $(G[i]-Gfit[i])")
-end
\ No newline at end of file
diff --git a/src/operation.jl b/src/operation.jl
index 5911619..b0eea2c 100644
--- a/src/operation.jl
+++ b/src/operation.jl
@@ -122,9 +122,14 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
     Nτ, Nω = size(kernel)
     @assert size(Gτ)[1] == Nτ
     if isnothing(sumrule) == false
+        @assert dlrGrid.symmetry == :none && dlrGrid.isFermi "only unsymmetrized ferminoic sum rule has been implemented!"
         # println(size(Gτ))
-        kernel_m0 = kernel[:, end]
-        kernel = kernel[:, 1:Nω-1] #a copy of kernel submatrix will be created
+        M = Int(floor(dlrGrid.size / 2))
+        # M = dlrGrid.size
+
+        kernel_m0 = kernel[:, M]
+        # kernel = kernel[:, 1:Nω-1] #a copy of kernel submatrix will be created
+        kernel = hcat(kernel[:, 1:M-1], kernel[:, M+1:end])
 
         for i in 1:Nτ
             Gτ[i, :] .-= kernel_m0[i] * sumrule
@@ -161,8 +166,9 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
         #add back the coeff that are fixed by the sum rule
         coeffmore = sumrule' .- sum(coeff, dims = 1)
         cnew = zeros(eltype(coeff), size(coeff)[1] + 1, size(coeff)[2])
-        cnew[1:end-1, :] = coeff
-        cnew[end, :] = coeffmore
+        cnew[1:M-1, :] = coeff[1:M-1, :]
+        cnew[M+1:end, :] = coeff[M:end, :]
+        cnew[M, :] = coeffmore
         return cnew
     else
         return coeff