From c274ecfb94e2d6ed90ed10de21013ea619c0591d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20P=C3=A9rez=20Panadero?= Date: Tue, 19 Dec 2023 12:31:30 +0100 Subject: [PATCH 1/5] DiracIO fixed for twisted mass --- docs/src/solvers.md | 2 +- src/Dirac/DiracIO.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/solvers.md b/docs/src/solvers.md index b6fc606..e6425da 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -76,7 +76,7 @@ Where $$P_\pm = (1 \pm \gamma_0)/2$$. The boundary-to-boundary propagator ```math - - \frac{c_t}{\sqrt{V}} \sum_\vec{x} U_0 ^\dagger (T-1,\vec{x}) P_+ S(T-1,\vec{x}) +\frac{-c_t}{\sqrt{V}} \sum_{\vec{x}} U_0 ^\dagger (T-1,\vec{x}) P_+ S(T-1,\vec{x}) ``` is computed by the function [`bndtobnd`](@ref) diff --git a/src/Dirac/DiracIO.jl b/src/Dirac/DiracIO.jl index 0c5da46..1342f70 100644 --- a/src/Dirac/DiracIO.jl +++ b/src/Dirac/DiracIO.jl @@ -41,7 +41,7 @@ function read_prop(fname::String) footh = Vector{Float64}(undef, 4) lp = SpaceParm{ndim}(iL, (4,4,4,4), ibc, ntw) - dpar = DiracParam{Float64}(SU3fund,foopars[1],foopars[2],ntuple(i -> footh[i], 4),foopars[3]) + dpar = DiracParam{Float64}(SU3fund,foopars[1],foopars[2],ntuple(i -> footh[i], 4),foopars[3],foopars[4]) dtr = (2,3,4,1) @@ -100,7 +100,7 @@ function save_prop(fname::String, psi, lp::SpaceParm{4,M,B,D}, dpar::DiracParam; BDIO_write!(fb, [convert(Int32, B)]) BDIO_write!(fb, [convert(Int32, lp.iL[i]) for i in 1:4]) BDIO_write!(fb, [convert(Int32, lp.ntw[i]) for i in 1:M]) - BDIO_write!(fb, [dpar.m0, dpar.csw, dpar.ct]) + BDIO_write!(fb, [dpar.m0, dpar.csw, dpar.tm, dpar.ct]) BDIO_write!(fb, [dpar.th[i] for i in 1:4]) end BDIO_write_hash!(fb) @@ -175,9 +175,9 @@ function read_dpar(fname::String) footh = Vector{Float64}(undef, 4) lp = SpaceParm{ndim}(iL, (4,4,4,4), ibc, ntw) - dpar = DiracParam{Float64}(SU3fund,foopars[1],foopars[2],ntuple(i -> footh[i], 4),foopars[3]) + dpar = DiracParam{Float64}(SU3fund,foopars[1],foopars[2],ntuple(i -> footh[i], 4),foopars[3],foopars[4]) BDIO_close!(fb) return dpar, lp -end \ No newline at end of file +end From b129e77878fbc632fac205138288cd23f38a4898 Mon Sep 17 00:00:00 2001 From: Fernando Desktop1 Date: Fri, 9 Feb 2024 15:11:13 +0100 Subject: [PATCH 2/5] Typo in dot() --- src/Groups/FundamentalSU3.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Groups/FundamentalSU3.jl b/src/Groups/FundamentalSU3.jl index d429d5d..d4d203b 100644 --- a/src/Groups/FundamentalSU3.jl +++ b/src/Groups/FundamentalSU3.jl @@ -38,7 +38,7 @@ norm2(a::SU3fund{T}) where T <: AbstractFloat = (abs2(a.t1) + abs2 Returns the scalar product of two fundamental elements. The convention is for the product to the linear in the second argument, and anti-linear in the first argument. """ -dot(g1::SU3fund{T},g2::SU3fund{T}) where T <: AbstractFloat = conj(g1.t1)*g2.t1+g1.t2*conj(g2.t2)+g1.t3*conj(g2.t3) +dot(g1::SU3fund{T},g2::SU3fund{T}) where T <: AbstractFloat = conj(g1.t1)*g2.t1+conj(g1.t2)*g2.t2+conj(g1.t3)*g2.t3 """ *(g::SU3{T},b::SU3fund{T}) From 506e8b53cb8542a7d7d000d0d464dbef33fad905 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20P=C3=A9rez=20Panadero?= Date: Wed, 14 Feb 2024 12:03:07 +0100 Subject: [PATCH 3/5] Correction in fA/fP test --- test/dirac/test_fp_fa.jl | 210 ++++++++++++++++----------------------- 1 file changed, 83 insertions(+), 127 deletions(-) diff --git a/test/dirac/test_fp_fa.jl b/test/dirac/test_fp_fa.jl index ad93afb..5551da5 100644 --- a/test/dirac/test_fp_fa.jl +++ b/test/dirac/test_fp_fa.jl @@ -2,164 +2,120 @@ using LatticeGPU using CUDA using TimerOutputs - @timeit "fA_fP test" begin -function fP_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) + function fP_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) - @timeit "fP inversion (x12)" begin + @timeit "fP inversion (x12)" begin - lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); - exptheta = exp.(im.*theta./lp.iL); + lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); + exptheta = exp.(im.*theta./lp.iL); + dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); + dws = DiracWorkspace(SU3fund,Float64,lp); - dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); - dws = DiracWorkspace(SU3fund,Float64,lp); + U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); + psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); - U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); - psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); + res = zeros(lp.iL[4]) - res = zeros(lp.iL[4]) + for s in 1:4 for c in 1:3 + bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); - for s in 1:4 for c in 1:3 - bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); - - for t in 1:lp.iL[4] - #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] + for t in 1:lp.iL[4] + #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; CUDA.@allowscalar (res[t] += norm2(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...])/2) - #end end end - #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + #end end end + #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + + end + + end end end - end end - - end - - @timeit "fP analitical solution" begin + @timeit "fP analitical solution" begin #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.33) - res_th = zeros(lp.iL[4]) + res_th = zeros(lp.iL[4]) - pp3 = ntuple(i -> theta[i]/lp.iL[i],3) - omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) - pp = (-im*omega,pp3...) - Mpp = m + 2* sum((sin.(pp./2)).^2) - Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) - - for i in 2:lp.iL[4] - res_th[i] = (2*3*sinh(omega)/(Rpp^2)) * ( (Mpp + sinh(omega))*exp(-2*omega*(i-1)) - (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))) ) - end - - end - return sum(abs.(res-res_th)) - -end - -function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) - - @timeit "fA inversion (x12)" begin - - lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); - exptheta = exp.(im.*theta./lp.iL); - - dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); - dws = DiracWorkspace(SU3fund,Float64,lp); - - U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); - psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); - - res = im*zeros(lp.iL[4]) - - for s in 1:4 for c in 1:3 - bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); - - for t in 1:lp.iL[4] - #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] - i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; - CUDA.@allowscalar (res[t] += -dot(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...],dmul(Gamma{4},psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...]))/2) - #end end end - #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + pp3 = ntuple(i -> theta[i]/lp.iL[i],3) + omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) + pp = (-im*omega,pp3...) + Mpp = m + 2* sum((sin.(pp./2)).^2) + Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) + for i in 2:lp.iL[4] + res_th[i] = (2*3*sinh(omega)/(Rpp^2)) * ( (Mpp + sinh(omega))*exp(-2*omega*(i-1)) - (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))) ) end - end end + end + return sum(abs.(res-res_th)) - lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); - exptheta = exp.(im.*theta./lp.iL); - - dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0); - dws = DiracWorkspace(SU3fund,Float64,lp); - - U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); - psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); - - res = im*zeros(lp.iL[4]) - - for s in 1:4 for c in 1:3 - bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); - - for t in 1:lp.iL[4] - #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] + end + + function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) + + @timeit "fA inversion (x12)" begin + + lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); + exptheta = exp.(im.*theta./lp.iL); + + dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); + dws = DiracWorkspace(SU3fund,Float64,lp); + + U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); + psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); + + res = im*zeros(lp.iL[4]) + + for s in 1:4 for c in 1:3 + bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); + + for t in 1:lp.iL[4] + #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; CUDA.@allowscalar (res[t] += -dot(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...],dmul(Gamma{4},psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...]))/2) - #end end end - #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) - + #end end end + #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + + end + + end end + end - - end end - -end - #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.32) - -@timeit "fA analitical solution" begin - res_th = zeros(lp.iL[4]) - - pp3 = ntuple(i -> theta[i]/lp.iL[i],3) - omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) - pp = (-im*omega,pp3...) - Mpp = m + 2* sum((sin.(pp./2)).^2) - Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) - - for i in 2:lp.iL[4] - res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) - - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))))) - end #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.32) - @timeit "fA analitical solution" begin - res_th = zeros(lp.iL[4]) + @timeit "fA analitical solution" begin + res_th = zeros(lp.iL[4]) - pp3 = ntuple(i -> theta[i]/lp.iL[i],3) - omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) - pp = (-im*omega,pp3...) - Mpp = m + 2* sum((sin.(pp./2)).^2) - Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) + pp3 = ntuple(i -> theta[i]/lp.iL[i],3) + omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) + pp = (-im*omega,pp3...) + Mpp = m + 2* sum((sin.(pp./2)).^2) + Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) + + for i in 2:lp.iL[4] + res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) + - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))))) + end - for i in 2:lp.iL[4] - res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) - - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))))) end - + return sum(abs.(res-res_th)) end - return sum(abs.(res-res_th)) - -end - - -difA = fA_test(); -difP = fP_test(); - -if difA > 1.0e-15 - error("fA test failed with error ", difA) -elseif difP > 1.0e-15 - error("fP test failed with error ", difP) -else - print("fA & fP tests passed with errors: ", difA," and ",difP,"!\n") -end + + difA = fA_test(); + difP = fP_test(); + + if difA > 1.0e-15 + error("fA test failed with error ", difA) + elseif difP > 1.0e-15 + error("fP test failed with error ", difP) + else + print("fA & fP tests passed with errors: ", difA," and ",difP,"!\n") + end end From 6214577ed1d5d7368d9eadee529b39fd79ea5cdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20P=C3=A9rez=20Panadero?= Date: Mon, 19 Feb 2024 10:15:25 +0100 Subject: [PATCH 4/5] Propagator functions return the number of iterations --- src/Solvers/Propagators.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/Solvers/Propagators.jl b/src/Solvers/Propagators.jl index 58008b2..87bd4df 100644 --- a/src/Solvers/Propagators.jl +++ b/src/Solvers/Propagators.jl @@ -5,7 +5,7 @@ function propagator!(pro,U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) Saves the fermionic progapator in pro for a source at point `y` with color `c` and spin `s`. If the last three arguments are replaced by `time::Int64`, the source is replaced -by a random source in spin and color at t = `time`. +by a random source in spin and color at t = `time`. Returns the number of iterations. """ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T} @@ -27,8 +27,8 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp) - CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return nothing + niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) + return niter end function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T} @@ -48,8 +48,8 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp) - CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return nothing + niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) + return niter end """ @@ -57,7 +57,7 @@ end function bndpropagator!(pro,U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) Saves the propagator from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's' in 'pro'. The factor c_t is included while the factor 1/sqrt(V) is not. -For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) +For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64). Returns the number of iterations. """ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D} @@ -94,8 +94,8 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp) - CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return nothing + niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) + return niter end """ @@ -103,7 +103,7 @@ end function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) Returns the propagator from the t=T boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not. -For the propagator from t=0 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) +For the propagator from t=0 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64). Returns the number of iterations. """ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D} @@ -139,8 +139,8 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp) - CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return nothing + niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) + return niter end From b92f9c92e0b825f6510359e954c6bc4f0a13688e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20P=C3=A9rez=20Panadero?= Date: Mon, 19 Feb 2024 10:23:37 +0100 Subject: [PATCH 5/5] Typo in read_cnfg for SF bc --- src/YM/YMio.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/YM/YMio.jl b/src/YM/YMio.jl index f492497..b8a15d8 100644 --- a/src/YM/YMio.jl +++ b/src/YM/YMio.jl @@ -75,7 +75,7 @@ function read_cnfg(fname::String) end if ibc == BC_SF_AFWB || ibc == BC_SF_ORBI - BDIO_read(fb, V) + BDIO_read(fb, vec(V)) Ubnd = ntuple(i->assign(i, V, 1), 3) BDIO_close!(fb)