Merge branch 'master' into 'master'

Master

See merge request alramos/latticegpu.jl!4

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)

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,7 +175,7 @@ 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)

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})

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)
+    niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return nothing
+    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)
+    niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return nothing
+    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)
+    niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return nothing
+    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)
+    niter = CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
-    return nothing
+    return niter
 end
 
 

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)
 

test/dirac/test_fp_fa.jl

@@ -2,17 +2,15 @@ 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
 
             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);
 
@@ -56,9 +54,9 @@ function fP_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1
         end
         return sum(abs.(res-res_th))
 
-end
+    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)
+    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
 
@@ -87,46 +85,6 @@ function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1
 
             end end
 
-    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]
-                    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
-    
-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)
 
@@ -145,21 +103,19 @@ end
             end
 
         end
-
         return sum(abs.(res-res_th))
-
+    end
-end
 
 
-difA = fA_test();
+    difA = fA_test();
-difP = fP_test();
+    difP = fP_test();
 
-if difA > 1.0e-15 
+        if difA > 1.0e-15
             error("fA test failed with error ", difA)
-elseif difP > 1.0e-15
+        elseif difP > 1.0e-15
             error("fP test failed with error ", difP)
-else
+        else
             print("fA & fP tests passed with errors: ", difA," and ",difP,"!\n")
-end
+        end
 
 end