diff --git a/.gitignore b/.gitignore index 3af67b1..1c02e5e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.jl.*.cov *.jl.cov *.jl.mem +Manifest.toml /docs/build/ diff --git a/Manifest.toml b/Manifest.toml deleted file mode 100644 index 2f54168..0000000 --- a/Manifest.toml +++ /dev/null @@ -1,217 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -[[ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" - -[[Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[ChainRulesCore]] -deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "be770c08881f7bb928dfd86d1ba83798f76cf62a" -uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.10.9" - -[[Compat]] -deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.31.0" - -[[CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" - -[[Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[DelimitedFiles]] -deps = ["Mmap"] -uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" - -[[Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" - -[[DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.5" - -[[Downloads]] -deps = ["ArgTools", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" - -[[FastGaussQuadrature]] -deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "5829b25887e53fb6730a9df2ff89ed24baa6abf6" -uuid = "442a2c76-b920-505d-bb47-c5924d526838" -version = "0.4.7" - -[[InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[JLLWrappers]] -deps = ["Preferences"] -git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.3.0" - -[[LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" - -[[LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" - -[[LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" - -[[Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[LinearAlgebra]] -deps = ["Libdl"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[LogExpFunctions]] -deps = ["DocStringExtensions", "LinearAlgebra"] -git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.2.5" - -[[Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" - -[[Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" - -[[NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" - -[[OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" - -[[Preferences]] -deps = ["TOML"] -git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.2" - -[[Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[Quadmath]] -deps = ["Printf", "Random", "Requires"] -git-tree-sha1 = "5a8f74af8eae654086a1d058b4ec94ff192e3de0" -uuid = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" -version = "0.5.5" - -[[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[Random]] -deps = ["Serialization"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[Requires]] -deps = ["UUIDs"] -git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621" -uuid = "ae029012-a4dd-5104-9daa-d747884805df" -version = "1.1.3" - -[[SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" - -[[Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[SharedArrays]] -deps = ["Distributed", "Mmap", "Random", "Serialization"] -uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" - -[[Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[SparseArrays]] -deps = ["LinearAlgebra", "Random"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" - -[[SpecialFunctions]] -deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"] -git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "1.5.1" - -[[StaticArrays]] -deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "745914ebcd610da69f3cb6bf76cb7bb83dcb8c9a" -uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.2.4" - -[[Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" - -[[TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" - -[[Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" - -[[Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" - -[[nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" - -[[p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" diff --git a/Project.toml b/Project.toml index 03f5db5..fb48cec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Lehmann" uuid = "95bf888a-8996-4655-9f35-1c0506bdfefe" authors = ["Kun Chen", "Tao Wang", "Xiansheng Cai"] -version = "0.2.1" +version = "0.2.2" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/build.jl b/build.jl index c7e6dd6..d226f1a 100644 --- a/build.jl +++ b/build.jl @@ -3,13 +3,14 @@ using Lehmann using Printf # rtol = [-6, -8, -10, -12, -14] -rtol = [-7, -9, -11, -13] -Λ = [100, 1000, 10000, 100000, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14] +# rtol = [-7, -9, -11, -13] +# Λ = [100, 1000, 10000, 100000, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14] # Λ = [1e11, 1e12, 1e13, 1e14] -# rtol = [-12] -# Λ = [1000000,] +rtol = [-12] +Λ = [1000000,] algorithm = :functional -folder = "./basis/" +# folder = "./basis/" +folder = "./" for lambda in Λ for err in rtol diff --git a/example/SYK.jl b/example/SYK.jl index f7f8627..e41d18f 100644 --- a/example/SYK.jl +++ b/example/SYK.jl @@ -17,7 +17,7 @@ A SYK model solver based on a forward fixed-point iteration method. The SYK Green's function has particle-hole symmetry when μ=0. You may enforce such symmetry by setting `symmetry = :ph` when initialize the DLR grids. - A symmetrized solver tends to be more robust than the unsymmetrized counterpart. + A symmetrized solver tends to be more robust than a unsymmetrized one. """ using Lehmann @@ -50,7 +50,7 @@ function dyson(d, sigma_q, mu) end end -function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 1000, G_x = zeros(ComplexF64, length(d))) +function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 5000, G_x = zeros(ComplexF64, length(d)), verbose = true) for iter in 1:maxiter @@ -61,8 +61,10 @@ function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxit G_x_new = matfreq2tau(d, G_q_new) - if iter % (maxiter / 10) == 0 - println("round $iter: change $(diff(G_x_new, G_x))") + if verbose + if iter % (maxiter / 10) == 0 + println("round $iter: change $(diff(G_x_new, G_x))") + end end if maximum(abs.(G_x_new .- G_x)) < tol && iter > 10 break @@ -83,26 +85,51 @@ function printG(d, G_x) println() end -printstyled("===== Symmetrized DLR solver for SYK model =======\n", color = :yellow) -mix = 0.1 -dsym = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-8, symmetry = :ph) # Initialize DLR object -G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix) -printG(dsym, G_x_ph) - -printstyled("===== Unsymmetrized DLR solver for SYK model =======\n", color = :yellow) -mix = 0.02 -dnone = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-8, symmetry = :none) # Initialize DLR object -G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix) -printG(dnone, G_x_none) - -printstyled("===== Unsymmetrized versus Symmetrized DLR solver =======\n", color = :yellow) -@printf("%15s%40s%40s%40s\n", "τ", "sym DLR (interpolated)", "unsym DLR", "difference") -G_x_interp = tau2tau(dsym, G_x_ph, dnone.τ) -for i in 1:dnone.size - if dnone.τ[i] <= dnone.β / 2 - @printf("%15.8f%40.15f%40.15f%40.15f\n", dnone.τ[i], real(G_x_interp[i]), real(G_x_none[i]), abs(real(G_x_interp[i] - G_x_none[i]))) - end +verbose = false + +printstyled("===== Prepare the expected Green's function of the SYK model =======\n", color = :yellow) +dsym_correct = DLRGrid(Euv = 5.0, β = 10000.0, isFermi = true, rtol = 1e-14, symmetry = :ph) # Initialize DLR object +G_x_correct = solve_syk_with_fixpoint_iter(dsym_correct, 0.00, mix = 0.1, verbose = false) +printG(dsym_correct, G_x_correct) + +printstyled("===== Test Symmetrized and Unsymmetrized DLR solver for SYK model =======\n", color = :yellow) + +@printf("%30s%30s%30s%15s\n", "Euv", "symmetrized solver error", "unsymmetrized solver error", "good or bad") +for Euv in LinRange(5.0, 10.0, 50) + + rtol = 1e-10 + β = 10000.0 + # printstyled("===== Symmetrized DLR solver for SYK model =======\n", color = :yellow) + mix = 0.01 + dsym = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :ph, rebuild = true, verbose = false) # Initialize DLR object + G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, verbose = verbose) + # printG(dsym, G_x_ph) + + # printstyled("===== Unsymmetrized DLR solver for SYK model =======\n", color = :yellow) + mix = 0.01 + dnone = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :none, rebuild = true, verbose = false) # Initialize DLR object + G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, verbose = verbose) + # printG(dnone, G_x_none) + + # printstyled("===== Unsymmetrized versus Symmetrized DLR solver =======\n", color = :yellow) + # @printf("%15s%40s%40s%40s\n", "τ", "sym DLR (interpolated)", "unsym DLR", "difference") + # G_x_interp = tau2tau(dsym_correct, G_x_correct, dnone.τ) + # for i in 1:dnone.size + # if dnone.τ[i] <= dnone.β / 2 + # @printf("%15.8f%40.15f%40.15f%40.15f\n", dnone.τ[i], real(G_x_interp[i]), real(G_x_none[i]), abs(real(G_x_interp[i] - G_x_none[i]))) + # end + # end + + G_x_interp_ph = tau2tau(dsym_correct, G_x_correct, dsym.τ) + G_x_interp_none = tau2tau(dsym_correct, G_x_correct, dnone.τ) + d_ph = diff(G_x_interp_ph, G_x_ph) + d_none = diff(G_x_interp_none, G_x_none) + flag = (d_ph < 100rtol) && (d_none < 100rtol) ? "good" : "bad" + + @printf("%30.15f%30.15e%30.15e%20s\n", Euv, d_ph, d_none, flag) + # println("symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_ph)) + # println("non symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_none)) + end -println("maximumal difference: ", diff(G_x_interp, G_x_none)) diff --git a/src/dlr.jl b/src/dlr.jl index 9464aa3..2bdd8ff 100644 --- a/src/dlr.jl +++ b/src/dlr.jl @@ -42,8 +42,8 @@ struct DLRGrid τ::Vector{Float64} """ - function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional, verbose = 0) - function DLRGrid(; Euv, β, isFermi::Bool, symmetry::Symbol = :none, rtol = 1e-8, rebuild = false, folder = nothing, algorithm = :functional, verbose = 0) + function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional, verbose = false) + function DLRGrid(; Euv, β, isFermi::Bool, symmetry::Symbol = :none, rtol = 1e-8, rebuild = false, folder = nothing, algorithm = :functional, verbose = false) Create DLR grids @@ -54,9 +54,10 @@ struct DLRGrid - `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` : the folder to load the DLR file if rebuild = false, or the folder to save the DLR file if rebuild = true + - `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` : 0 not to print DLRGrid to terminal, >0 to print + - `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 @@ -73,39 +74,80 @@ struct DLRGrid @warn("Current implementation may cause ~ 3-4 digits loss for rtol ≥ 1e-6!") end - if Λ < 100 - Λ = Int(100) - else - Λ = 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n - end - 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 - if symmetry == :none - # if isFermi - filename = "universal_$(Λ)_1e$(rtolpower).dlr" - # else - # error("Generic bosonic dlr has not yet been implemented!") - # end - elseif symmetry == :ph - filename = "ph_$(Λ)_1e$(rtolpower).dlr" - elseif symmetry == :pha - filename = "pha_$(Λ)_1e$(rtolpower).dlr" - else - error("$symmetry is not implemented!") + function filename(lambda, errpower) + lambda = Int(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 - dlr = new(isFermi, symmetry, Euv, β, Λ, rtol, [], [], [], []) - if rebuild - _build!(dlr, folder, filename, algorithm, verbose) - else - _load!(dlr, folder, filename, algorithm, verbose) + + if rebuild == false + if isnothing(folder) + Λfloor = Λ < 100 ? Int(100) : 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n + + folderList = [string(@__DIR__, "/../basis/"),] + file = filename(Λfloor, rtolpower) + + dlrfile = finddlr(folderList, file) + + if isnothing(dlrfile) == false + dlr = new(isFermi, symmetry, Euv, β, Λfloor, 10.0^(float(rtolpower)), [], [], [], []) + _load!(dlr, dlrfile, verbose) + return dlr + else + @warn("No DLR is found in the folder $folder, try to rebuild instead.") + end + else + file = filename(Euv * β, rtolpower) + folderList = [folder,] + + dlrfile = finddlr(folderList, file) + + if isnothing(dlrfile) == false + dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], []) + _load!(dlr, dlrfile, verbose) + return dlr + else + @warn("No DLR is found in the folder $folder, try to rebuild instead.") + end + end + end + + # try to rebuild the dlrGrid + dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], []) + file2save = filename(Euv * β, rtolpower) + _build!(dlr, folder, file2save, algorithm, verbose) return dlr end + function DLRGrid(; Euv, β, isFermi::Bool, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false) return DLRGrid(Euv, β, rtol, isFermi, symmetry; rebuild = rebuild, folder = folder, algorithm = algorithm, verbose = verbose) end @@ -133,29 +175,7 @@ Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω) Base.length(dlrGrid::DLRGrid) = length(dlrGrid.ω) rank(dlrGrid::DLRGrid) = length(dlrGrid.ω) -function _load!(dlrGrid::DLRGrid, folder, filename, algorithm = :functional, verbose = false) - searchdir(path, key) = filter(x -> occursin(key, x), readdir(path)) - - function finddlr(folder, filename) - 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 - - folder = isnothing(folder) ? [] : collect(folder) - push!(folder, string(@__DIR__, "/../basis/")) - - dlrfile = finddlr(folder, filename) - - if isnothing(dlrfile) - _build!(dlrGrid, nothing, nothing, algorithm) - return - end +function _load!(dlrGrid::DLRGrid, dlrfile, verbose = false) grid = readdlm(dlrfile, comments = true, comment_char = '#') # println("reading $filename") @@ -205,6 +225,7 @@ function _build!(dlrGrid::DLRGrid, folder, filename, algorithm, verbose = false) push!(dlrGrid.n, n) push!(dlrGrid.ωn, isFermi ? (2n + 1.0) * π / β : 2n * π / β) end + # println(rank) end diff --git a/src/functional/builder.jl b/src/functional/builder.jl index 9e66a6f..f387cfd 100644 --- a/src/functional/builder.jl +++ b/src/functional/builder.jl @@ -107,13 +107,13 @@ function printCandidate(basis, idx) @printf("%3i : ω=%24.8f ∈ (%24.8f, %24.8f) -> error=%24.16g\n", basis.N, basis.grid[idx], lower, upper, basis.residual[idx]) end -function QR(Λ, rtol, proj, g0; N = nothing) +function QR(Λ, rtol, proj, g0; N = nothing, verbose = false) basis = Basis(Λ, rtol) # println(g0) for g in g0 idx = addBasis!(basis, proj, Float(g)) # @printf("%3i : ω=%24.8f ∈ (%24.8f, %24.8f) -> error=%24.16g\n", 1, g, 0, Λ, basis.residual[idx]) - printCandidate(basis, idx) + verbose && printCandidate(basis, idx) end # @code_warntype Residual(basis, proj, Float(1.0)) @@ -125,7 +125,7 @@ function QR(Λ, rtol, proj, g0; N = nothing) newω = basis.candidate[ωi] idx = addBasis!(basis, proj, newω) - printCandidate(basis, idx) + verbose && printCandidate(basis, idx) # println(length(basis.grid)) # println(idx) # lower = (idx == 1) ? 0 : basis.grid[idx - 1] @@ -136,9 +136,9 @@ function QR(Λ, rtol, proj, g0; N = nothing) # plotResidual(basis, proj, Float(0), Float(100), candidate, residual) maxResidual, ωi = findmax(basis.candidateResidual) end - testOrthgonal(basis) + testOrthgonal(basis, verbose) # @printf("residual = %.10e, Fnorm/F0 = %.10e\n", residual, residualF(freq, Q, Λ)) - @printf("residual = %.10e\n", maximum(basis.candidateResidual)) + verbose && @printf("residual = %.10e\n", maximum(basis.candidateResidual)) # plotResidual(basis, proj, Float(0), Float(100), basis.candidate, basis.candidateResidual) return basis end @@ -215,11 +215,11 @@ function Residual(basis, proj, g::Float) end -function testOrthgonal(basis) - println("testing orthognalization...") +function testOrthgonal(basis, verbose) + # println("testing orthognalization...") II = basis.Q * basis.proj * basis.Q' maxerr = maximum(abs.(II - I)) - println("Max Orthognalization Error: ", maxerr) + verbose && println("Max Orthognalization Error: ", maxerr) end """