Fix bug when generating DLR grid at small Lambda

Author: fsxbhyy

example/hightemperature.jl

@@ -63,9 +63,11 @@ function τnQR(kernel, rank, print::Bool=true)
     @assert rank == size(kernel)[2] #the ω dimension of the τkernel should be the effective rank
 
     τnqr = qr(transpose(kernel), Val(true)) # julia qr has a strange, Val(true) will do a pivot QR
-
-    return τnqr.p[1:rank]
-
+    if rank > length(τnqr.p)
+        return τnqr.p
+    else
+        return τnqr.p[1:rank]
+    end
 end
 
 function buildωn(dlrGrid, print::Bool=true)
@@ -120,7 +122,6 @@ function build(dlrGrid, print::Bool=true)
     print && println("ω grid size = $(ω.ngrid)")
 
     kernel = preciseKernelT(dlrGrid, τ, ω, print)
-    println("test $(size(kernel))  $(length(τ.grid)) $(length(ω.grid ))\n")
     if dlrGrid.symmetry != :sym
         testInterpolation(dlrGrid, τ, ω, kernel, print)
     end

src/discrete/builder.jl

@@ -25,7 +25,9 @@ function ωChebyGrid(dlrGrid, degree, print=true)
     Λ, rtol = dlrGrid.Λ, dlrGrid.rtol
 
     npo = Int(ceil(log(Λ) / log(2.0))) # subintervals on [0,lambda] in omega space (subintervals on [-lambda,lambda] is 2*npo)
-
+    if npo < 1
+        npo = 1
+    end
     if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha
         # Panel break points for the real frequency ∈ [0, Λ]
         # get exponentially dense near 0⁺
@@ -52,8 +54,8 @@ function τChebyGrid(dlrGrid, degree, print=true)
     Λ, rtol = dlrGrid.Λ, dlrGrid.rtol
 
     npt = Int(ceil(log(Λ) / log(2.0))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt)
-    if npt < 2
-        npt = 2
+    if npt < 1
+        npt = 1
     end
 
     if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha

src/discrete/kernel.jl

@@ -343,20 +343,40 @@ function _build!(dlrGrid::DLRGrid, folder, filename, algorithm, verbose=false)
     end
     rank = length(ω)
     if isnothing(folder) == false
-        open(joinpath(folder, filename), "w") do io
-            @printf(io, "# %5s  %25s  %25s  %25s  %20s\n", "index", "freq", "tau", "fermi n", "bose n")
-            for r = 1:rank
-                @printf(io, "%5i  %32.17g  %32.17g  %16i  %16i\n", r, ω[r], τ[r], nF[r], nB[r])
-            end
+        nan = "NAN"
+        print("$(filename)\n")
+        file = open(joinpath(folder, filename), "w")
+        #open(joinpath(folder, filename), "w") do io
+        @printf(file, "# %5s  %25s  %25s  %25s  %20s\n", "index", "freq", "tau", "fermi n", "bose n")
+        for r = 1:rank
+            s0 = "%5i "
+            s1 = r > length(ω) ? "%48s " : "%48.40g "
+            s2 = r > length(τ) ? "%48s " : "%48.40g "
+            s3 = r > length(nF) ? "%16s " : "%16i "
+            s4 = r > length(nB) ? "%16s\n" : "%16i\n"
+            f = Printf.Format(s0 * s1 * s2 * s3 * s4)
+            Printf.format(file, f, r, r > length(ω) ? nan : ω[r],
+                r > length(τ) ? nan : τ[r],
+                r > length(nF) ? nan : nF[r],
+                r > length(nB) ? nan : nB[r])
         end
+        # for r = 1:rank
+        #     @printf(io, "%5i  %32.17g  %32.17g  %16i  %16i\n", r, ω[r], τ[r], nF[r], nB[r])
+        # end
+        close(file)
     end
-    for r = 1:rank
-        push!(dlrGrid.ω, ω[r] / β)
-        push!(dlrGrid.τ, τ[r] * β)
-        n = isFermi ? nF[r] : nB[r]
-        push!(dlrGrid.n, n)
-        push!(dlrGrid.ωn, isFermi ? (2n + 1.0) * π / β : 2n * π / β)
-    end
+    dlrGrid.ω = ω / β
+    dlrGrid.τ = τ * β
+    n = isFermi ? copy(nF) : copy(nB)
+    dlrGrid.n = n
+    dlrGrid.ωn = isFermi ? (2n .+ 1.0) * π / β : 2n * π / β
+    # for r = 1:rank
+    #     push!(dlrGrid.ω, ω[r] / β)
+    #     push!(dlrGrid.τ, τ[r] * β)
+    #     n = isFermi ? nF[r] : nB[r]
+    #     push!(dlrGrid.n, n)
+    #     push!(dlrGrid.ωn, isFermi ? (2n + 1.0) * π / β : 2n * π / β)
+    # end
     # println(rank)
 end
 

src/dlr.jl

@@ -6,8 +6,8 @@ SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=24, regu
 #using DoubleFloats
 #rtol(x, y) = maximum(abs.(x - y)) / maximum(abs.(x))
 
-function L2normτ(value_dlr, dlr, case, poles=nothing, weight = nothing)
-    function fine_τGrid(Λ::Float,degree,ratio::Float) where {Float}
+function L2normτ(value_dlr, dlr, case, poles=nothing, weight=nothing)
+    function fine_τGrid(Λ::Float, degree, ratio::Float) where {Float}
         ############## use composite grid #############################################
         # Generating a log densed composite grid with LogDensedGrid()
         npo = Int(ceil(log(Λ) / log(ratio))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt)
@@ -16,7 +16,7 @@ function L2normτ(value_dlr, dlr, case, poles=nothing, weight = nothing)
             [0.0, 1.0],# The grid is defined on [0.0, β]
             [0.0, 1.0],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter.
             npo,# N of log grid
-            0.5 / ratio^(npo-1), # minimum interval length of log grid
+            0.5 / ratio^(npo - 1), # minimum interval length of log grid
             degree, # N of bottom layer
             Float
         )
@@ -25,7 +25,7 @@ function L2normτ(value_dlr, dlr, case, poles=nothing, weight = nothing)
         # println("Composite expoential grid size: $(length(grid))")
         #println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[end])]")
         return grid
-    
+
         ############# DLR based fine grid ##########################################
         # dlr = DLRGrid(Euv=Float64(Λ), beta=1.0, rtol=Float64(rtol) / 100, isFermi=true, symmetry=:ph, rebuild=true)
         # # println("fine basis number: $(dlr.size)\n", dlr.ω)
@@ -36,22 +36,22 @@ function L2normτ(value_dlr, dlr, case, poles=nothing, weight = nothing)
         #     uniform = [panel[i] + (panel[i+1] - panel[i]) / degree * j for j in 0:degree-1]
         #     append!(grid, uniform)
         # end
-    
+
         # println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[2])]")
         # return grid
     end
     fineGrid = fine_τGrid(dlr.Euv, 12, typeof(dlr.Euv)(1.5))
-    value = real(tau2tau(dlr, value_dlr,  fineGrid, dlr.τ ))
+    value = real(tau2tau(dlr, value_dlr, fineGrid, dlr.τ))
     if case == MultiPole
         value_analy = case(dlr, fineGrid, :τ, poles, weight)
     else
         value_analy = case(dlr, fineGrid, :τ)
     end
     #print("value_analy $(value_analy[1:10])\n" )
-    interp = Interp.integrate1D( value , fineGrid)
-    interp_analy = Interp.integrate1D( value_analy , fineGrid)
+    interp = Interp.integrate1D(value, fineGrid)
+    interp_analy = Interp.integrate1D(value_analy, fineGrid)
     #print("$(interp_analy) $(interp_analy)\n")
-    return abs(interp_analy), abs(interp-interp_analy), maximum(abs.(value - value_analy))/ maximum(abs.(value_analy))
+    return abs(interp_analy), abs(interp - interp_analy), maximum(abs.(value - value_analy)) / maximum(abs.(value_analy))
 end
 #function MultiPole(dlr, grid, type)
 #    Euv = dlr.Euv
@@ -60,7 +60,7 @@ end
 #    return Sample.MultiPole(dlr, type, poles, grid; regularized=true)
 #end
 
-function MultiPole(dlr, grid, type, coeff, weight = nothing)
+function MultiPole(dlr, grid, type, coeff, weight=nothing)
     Euv = dlr.Euv
     poles = coeff * Euv
     # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true)
@@ -71,32 +71,32 @@ function MultiPole(dlr, grid, type, coeff, weight = nothing)
     end
 end
 
-function test_dlr_coeff(case, isFermi, symmetry, Euv, β, eps, eff_poles, weight; dtype=Float64, output = false)
+function test_dlr_coeff(case, isFermi, symmetry, Euv, β, eps, eff_poles, weight; dtype=Float64, output=false)
     para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps"
     dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis
     #dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose
     dlr10 = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct denser dlr basis for benchmark purpose
-    
+
     N_poles = 1000
     N = 1000
-   
+
     Gndlr = case(dlr, dlr.n, :n, eff_poles, weight)
     τSample = dlr10.τ
     Gsample = case(dlr, τSample, :τ, eff_poles, weight)
 
     Gfourier = matfreq2tau(dlr, Gndlr, τSample)
-    dlreff = matfreq2dlr(dlr,Gndlr)
+    dlreff = matfreq2dlr(dlr, Gndlr)
     dlreff = imag(dlreff)
     print("$(symmetry) $(Euv)  $(eps)  max $(maximum(abs.(dlreff) ))  min $(minimum(abs.(dlreff )))\n")
 end
-function test_err(case, isFermi, symmetry, Euv, β, eps, poles, weights; dtype=Float64, output = false)
+function test_err(case, isFermi, symmetry, Euv, β, eps, poles, weights; dtype=Float64, output=false)
     # println("Test $case with isFermi=$isFermi, Symmetry = $symmetry, Euv=$Euv, β=$β, rtol=$eps")
     #N_poles = 100
     para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps"
     dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis
     #dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose
     dlr10 = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct denser dlr basis for benchmark purpose
-    
+
     N_poles = size(poles)[2]
     N = size(poles)[1]
     value_sum = 0.0
@@ -106,63 +106,64 @@ function test_err(case, isFermi, symmetry, Euv, β, eps, poles, weights; dtype=F
     block = zeros(dtype, 10)
     if case == MultiPole
         for i in 1:N
-            eff_poles = poles[i,:]
-            weight = weights[i,:]
+            eff_poles = poles[i, :]
+            weight = weights[i, :]
             Gndlr = case(dlr, dlr.n, :n, eff_poles, weight)
             τSample = dlr10.τ
             Gsample = case(dlr, τSample, :τ, eff_poles, weight)
 
             Gfourier = matfreq2tau(dlr, Gndlr, τSample)
-            dlreff = matfreq2dlr(dlr,Gndlr)
-            value, err, max_err= L2normτ(Gfourier, dlr, case, eff_poles, weight)
+            dlreff = matfreq2dlr(dlr, Gndlr)
+            value, err, max_err = L2normτ(Gfourier, dlr, case, eff_poles, weight)
             modulus = abs(sum(dlreff))
-            value_sum +=value/modulus
-            err_sum += err/modulus
-            eta += err/value
+            value_sum += value / modulus
+            err_sum += err / modulus
+            eta += err / value
             max_err_sum += max_err
-            block[(i-1) ÷ (N÷10)+1] +=err/value/N*10
+            block[(i-1)÷(N÷10)+1] += err / value / N * 10
         end
     else
         Gndlr = case(dlr, dlr.n, :n)
         τSample = dlr10.τ
         Gsample = case(dlr, τSample, :τ)
         Gfourier = matfreq2tau(dlr, Gndlr, τSample)
-        dlreff = matfreq2dlr(dlr,Gndlr)
+        dlreff = matfreq2dlr(dlr, Gndlr)
         print("max $(maximum(dlreff))  min $(minimum(dlreff))\n")
-        value, err, max_err= L2normτ(Gfourier, dlr, case)
+        value, err, max_err = L2normτ(Gfourier, dlr, case)
         modulus = abs(sum(dlreff))
         print("test Semi: $(modulus)\n")
-        value_sum +=value/modulus
-        err_sum += err/modulus
+        value_sum += value / modulus
+        err_sum += err / modulus
         max_err_sum += max_err
     end
     if output
-        file = open("./accuracy_test1.dat", "a") 
+        file = open("./accuracy_test1.dat", "a")
         #@printf(file, "%48.40g  %48.40g %48.40g\n", eps, abs(b-c),  )
         #@printf(file, "%24.20g  %24.20g %24.20g %24.20g %24.20g %24.20g\n", eps, value_sum/N,  err_sum/N, err_sum /N/eps*Euv, max_err_sum/N, eta/N)
-        @printf(file, "%24.20g  %24.20g %24.20g\n", eps, log10(eta/N/eps),  std( log10.(block/eps) )  )
+        @printf(file, "%24.20g  %24.20g %24.20g\n", eps, log10(eta / N / eps), std(log10.(block / eps)))
         close(file)
     end
 end
 
-cases = [MultiPole]
+#cases = [MultiPole]
+cases = [SemiCircle]
 Λ = [1e4]
 #rtol = [ 1e-12]
-rtol = [ 1e-4, 1e-6, 1e-8, 1e-10, 1e-12]
+rtol = [1e-4, 1e-6, 1e-8, 1e-10, 1e-12]
 N_poles = 1000
 N = 1000
 setprecision(128)
 dtype = Float64
 #dtype = BigFloat
-poles = zeros(dtype,(N,N_poles) )
-weights = zeros(dtype, (N,N_poles))
+poles = zeros(dtype, (N, N_poles))
+weights = zeros(dtype, (N, N_poles))
 for i in 1:N
-    poles[i,:] = 2.0*rand(dtype, N_poles) .- 1.0
-    weights[i,:]  = 2.0 *rand(dtype, N_poles) .- 1.0
+    poles[i, :] = 2.0 * rand(dtype, N_poles) .- 1.0
+    weights[i, :] = 2.0 * rand(dtype, N_poles) .- 1.0
 end
 
 
-          
+
 for case in cases
     for l in Λ
         for r in rtol
@@ -177,11 +178,11 @@ for case in cases
             #     test(case, false, :pha, l, 1.0, r,dtype=BigFloat)
             #     test(case, true, :pha, l, 1.0, r, dtype= BigFloat)
             # end
-            
-            test_err(case, true, :none, l, 1.0, r,  poles ,  weights, dtype = Float64, output = true)
-           
-            test_err(case, true, :sym, l, 1.0, r, poles, weights, dtype = Float64, output = true)
-           # test_err(case, true, :sym, l, 1.0, r, poles, weights, dtype = BigFloat, output = true)
+
+            test_err(case, true, :none, l, 1.0, r, poles, weights, dtype=Float64, output=true)
+
+            test_err(case, true, :sym, l, 1.0, r, poles, weights, dtype=Float64, output=true)
+            # test_err(case, true, :sym, l, 1.0, r, poles, weights, dtype = BigFloat, output = true)
             #test_err(case, false, :sym, l, 1.0, r, dtype = BigFloat, output = true)
             # dtype =BigFloat
             # N_poles = 1000
@@ -193,10 +194,10 @@ for case in cases
             # test_dlr_coeff(case, true, :none, l, 1.0, r, eff_poles, weight, dtype = Float64, output = true)
             # test_dlr_coeff(case, true, :sym, l, 1.0, r, eff_poles , weight, dtype = Float64, output = true)
 
-           
+
         end
     end
 end
-    
+