Propagators file added

src/LatticeGPU.jl

@@ -63,5 +63,6 @@ export Dw!, DwdagDw!, pfrandomize!, g5Dw!
 include("Solvers/Solvers.jl")
 using .Solvers
 export CG!
+export propagator, bndpropagator
 
 end # module

src/Solvers/Propagators.jl

+
+
+
+
+
+
+"""
+    function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64)
+
+Returns the fermionic progapator  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`.
+
+"""
+function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T}
+
+    function krnlg5!(src)
+        b=Int64(CUDA.threadIdx().x)
+        r=Int64(CUDA.blockIdx().x)
+        src[b,r] = dmul(Gamma{5},src[b,r])
+        return nothing
+    end
+  
+    pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+    src = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
+    CUDA.@allowscalar src[point_index(CartesianIndex{lp.ndim}(y),lp)...] = Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4))
+       
+    CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
+    end
+       
+    g5Dw!(pro,U,src,dpar,lp)
+      
+    CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
+    return pro
+end
+
+
+function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T}
+
+    function krnlg5!(src)
+        b=Int64(CUDA.threadIdx().x)
+        r=Int64(CUDA.blockIdx().x)
+        src[b,r] = dmul(Gamma{5},src[b,r])
+        return nothing
+    end
+  
+    pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+    src = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+    pfrandomize!(src,lp,dpar,time)
+    
+    CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
+    end
+       
+    g5Dw!(pro,U,src,dpar,lp)
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
+    end
+      
+    CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
+    return pro
+end
+
+"""
+
+    function bndpropagator(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=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'.
+
+
+"""
+
+
+function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
+
+    function krnlg5!(src)
+        b=Int64(CUDA.threadIdx().x)
+        r=Int64(CUDA.blockIdx().x)
+        src[b,r] = dmul(Gamma{5},src[b,r])
+        return nothing
+    end
+
+    function krnl_assign_bndsrc!(src,U,lp::SpaceParm, c::Int64, s::Int64)
+        b=Int64(CUDA.threadIdx().x)
+        r=Int64(CUDA.blockIdx().x)
+
+        if (point_time((b,r),lp) == 2)
+        bd4, rd4 = dw((b,r), 4, lp) 
+        src[b,r] = gdagpmul(Pgamma{4,1},U[bd4,4,rd4],Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4)))/2
+        end
+
+        return nothing
+    end
+  
+    pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+    src = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
+    
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(src, U, lp, c, s)
+    end
+
+    CUDA.@sync begin
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
+    end
+       
+    g5Dw!(pro,U,src,dpar,lp)
+    
+    CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
+    return pro
+end
+

src/Solvers/Solvers.jl

@@ -22,7 +22,8 @@ using ..Dirac
 include("CG.jl")
 export CG!
 
-
+include("Propagators.jl")
+export propagator, bndpropagator
 
 
 end