Some updates

Manifest.toml

@@ -1,5 +1,10 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ANSIColoredPrinters]]
+git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
+uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
+version = "0.0.1"
+
 [[AbstractFFTs]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
@@ -18,6 +23,14 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
+[[BDIO]]
+deps = ["Documenter", "Nettle"]
+git-tree-sha1 = "31aa5298af4936ad33dd39f43ef0c3f7c45fb197"
+repo-rev = "master"
+repo-url = "https://gitlab.ift.uam-csic.es/alberto/bdio.jl.git"
+uuid = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
+version = "0.1.0"
+
 [[BFloat16s]]
 deps = ["LinearAlgebra", "Test"]
 git-tree-sha1 = "4af69e205efc343068dc8722b8dfec1ade89254a"
@@ -78,6 +91,12 @@ git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.8.5"
 
+[[Documenter]]
+deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233"
+uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+version = "0.27.10"
+
 [[Downloads]]
 deps = ["ArgTools", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -87,6 +106,10 @@ git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 version = "0.1.6"
 
+[[GMP_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
+
 [[GPUArrays]]
 deps = ["AbstractFFTs", "Adapt", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"]
 git-tree-sha1 = "ececbf05f8904c92814bdbd0aafd5540b0bf2e9a"
@@ -99,6 +122,12 @@ git-tree-sha1 = "e8a09182a4440489e2e3dedff5ad3f6bbe555396"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
 version = "0.12.5"
 
+[[IOCapture]]
+deps = ["Logging", "Random"]
+git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
+uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
+version = "0.2.2"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -109,6 +138,12 @@ git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
 version = "1.3.0"
 
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.2"
+
 [[LLVM]]
 deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"]
 git-tree-sha1 = "1b7ba36ea7aa6fa2278118951bad114fbb8359f2"
@@ -171,6 +206,18 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 [[MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 
+[[Nettle]]
+deps = ["Libdl", "Nettle_jll"]
+git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
+uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
+version = "0.5.1"
+
+[[Nettle_jll]]
+deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
+uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
+version = "3.7.2+0"
+
 [[NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
@@ -185,6 +232,12 @@ git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.4.1"
 
+[[Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.1.3"
+
 [[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"

Project.toml

@@ -4,6 +4,7 @@ authors = ["Alberto Ramos <alberto.ramos@ific.uv.es>"]
 version = "0.1.0"
 
 [deps]
+BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"

src/Dirac/Dirac.jl

@@ -19,6 +19,23 @@ using ..Fields
 using ..YM
 using ..Spinors
 
+struct DiracWorkspace{T}
+    sr
+    sp
+    sAp
+    st
+    
+    function DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D}) where {G,T <: AbstractFloat, B,D}
+
+        sr  = scalar_field(Spinor{4,G}, lp)
+        sp  = scalar_field(Spinor{4,G}, lp)
+        sAp = scalar_field(Spinor{4,G}, lp)
+        st  = scalar_field(Spinor{4,G}, lp)
+        return new{T}(sr,sp,sAp,st)
+    end
+end
+export DiracWorkspace
+
 function Dw!(so, U, si, m0, lp::SpaceParm)
 
     @timeit "Dw" begin
@@ -30,7 +47,7 @@ function Dw!(so, U, si, m0, lp::SpaceParm)
     return nothing
 end
 
-function DwdagDw!(so, U, si, m0, lp::SpaceParm)
+function DwdagDw!(so, U, si, m0, st, lp::SpaceParm)
 
     @timeit "DwdagDw" begin
         CUDA.@sync begin
@@ -49,9 +66,10 @@ function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
-<    # For SF:
+    # For SF:
     #  - cttilde affects mass term at x0 = a, T-a
     #  - Spinor can be periodic as long as 0 at x_0=0
+    @inbounds begin 
         so[b,r] = (4+m0)*si[b,r]
         for id in 1:4
             bu, ru = up((b,r), id, lp)
@@ -60,6 +78,7 @@ function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
             so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
                 conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
         end
+    end
 
     return nothing
 end
@@ -68,6 +87,7 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
+    @inbounds begin 
         so[b,r] = (4+m0)*si[b,r]
         for id in 1:4
             bu, ru = up((b,r), id, lp)
@@ -76,12 +96,12 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
             so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
                 conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
         end
-    so[b,r].s[3] = -so[b,r].s[3] 
+        so[b,r] = dmul(Gamma{5}, so[b,r])
-    so[b,r].s[4] = -so[b,r].s[4] 
+    end
         
     return nothing
 end
 
-export Dw
+export Dw!, DwdagDw!
 
 end

src/LatticeGPU.jl

@@ -45,7 +45,7 @@ export Eoft_clover, Eoft_plaq, Qtop
 export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
 export flw, flw_adapt
 export sfcoupling, bndfield, setbndfield
-export import_lex64, import_cern64
+export import_lex64, import_cern64, import_bsfqcd
 
 include("Spinors/Spinors.jl")
 
@@ -54,4 +54,14 @@ export Spinor, Pgamma
 export imm, mimm
 export pmul, gpmul, gdagpmul, dmul
 
+include("Dirac/Dirac.jl")
+using .Dirac
+export DiracWorkspace
+export Dw!, DwdagDw!
+
+
+include("Solvers/Solvers.jl")
+using .Solvers
+export CG!
+
 end # module

src/Spinors/Spinors.jl

@@ -22,7 +22,7 @@ Base.zero(::Type{Spinor{NS,G}}) where {NS,G} = Spinor{NS,G}(ntuple(i->zero(G), N
 """
     norm2(a::Spinor)
 
-Returns the norm of a fundamental element. Same result as sqrt(dot(a,a)).
+Returns the norm squared of a fundamental element. Same result as `dot(a,a)`.
 """
 @generated function norm2(a::Spinor{NS,G}) where {NS,G}
 
@@ -100,52 +100,52 @@ end
 
 Returns ``(1+s\\gamma_N)a``.
 """
-function pmul(::Type{Pgamma{4,1}},  a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{4,1}},  a::Spinor{4,G}) where {NS,G}
     
     r1 = a.s[1]+a.s[3]
     r2 = a.s[2]+a.s[4]
     return Spinor{4,G}((r1,r2,r1,r2))
 end
-function pmul(::Type{Pgamma{4,-1}}, a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{4,-1}}, a::Spinor{4,G}) where {NS,G}
     
     r1 = a.s[1]-a.s[3]
     r2 = a.s[2]-a.s[4] 
     return Spinor{4,G}((r1,r2,-r1,-r2))
 end
     
-function pmul(::Type{Pgamma{1,1}},  a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{1,1}},  a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]+imm(a.s[4])
     r2 = a.s[2]+imm(a.s[3])
     return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
 end
-function pmul(::Type{Pgamma{1,-1}}, a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{1,-1}}, a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]-imm(a.s[4])
     r2 = a.s[2]-imm(a.s[3])
     return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
 end
 
-function pmul(::Type{Pgamma{2,1}},  a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{2,1}},  a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]+a.s[4]
     r2 = a.s[2]-a.s[3]
     return Spinor{4,G}((r1,r2,-r2,r1))
 end
-function pmul(::Type{Pgamma{2,-1}}, a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{2,-1}}, a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]-a.s[4]
     r2 = a.s[2]+a.s[3]
     return Spinor{4,G}((r1,r2,r2,-r1))
 end
 
-function pmul(::Type{Pgamma{3,1}},  a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{3,1}},  a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]+imm(a.s[3])
     r2 = a.s[2]-imm(a.s[4])
     return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
 end
-function pmul(::Type{Pgamma{3,-1}},  a::Spinor{4,G}) where {NS,G}
+@inline function pmul(::Type{Pgamma{3,-1}},  a::Spinor{4,G}) where {NS,G}
 
     r1 = a.s[1]-imm(a.s[3])
     r2 = a.s[2]+imm(a.s[4])
@@ -158,52 +158,52 @@ end
 
 Returns ``g(1+s\\gamma_N)a``
 """
-function gpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
     
     r1 = g*(a.s[1]+a.s[3])
     r2 = g*(a.s[2]+a.s[4])
     return Spinor{4,G}((r1,r2,r1,r2))
 end
-function gpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
     
     r1 = g*(a.s[1]-a.s[3])
     r2 = g*(a.s[2]-a.s[4]) 
     return Spinor{4,G}((r1,r2,-r1,-r2))
 end
     
-function gpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]+imm(a.s[4]))
     r2 = g*(a.s[2]+imm(a.s[3]))
     return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
 end
-function gpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]-imm(a.s[4]))
     r2 = g*(a.s[2]-imm(a.s[3]))
     return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
 end
 
-function gpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]+a.s[4])
     r2 = g*(a.s[2]-a.s[3])
     return Spinor{4,G}((r1,r2,-r2,r1))
 end
-function gpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]-a.s[4])
     r2 = g*(a.s[2]+a.s[3])
     return Spinor{4,G}((r1,r2,r2,-r1))
 end
 
-function gpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]+imm(a.s[3]))
     r2 = g*(a.s[2]-imm(a.s[4]))
     return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
 end
-function gpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g*(a.s[1]-imm(a.s[3]))
     r2 = g*(a.s[2]+imm(a.s[4]))
@@ -215,52 +215,52 @@ end
 
 Returns ``g^+ (1+s\\gamma_N)a``
 """
-function gdagpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{4,1}},  g, a::Spinor{4,G}) where {NS,G}
     
     r1 = g\(a.s[1]+a.s[3])
     r2 = g\(a.s[2]+a.s[4])
     return Spinor{4,G}((r1,r2,r1,r2))
 end
-function gdagpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{4,-1}}, g, a::Spinor{4,G}) where {NS,G}
     
     r1 = g\(a.s[1]-a.s[3])
     r2 = g\(a.s[2]-a.s[4]) 
     return Spinor{4,G}((r1,r2,-r1,-r2))
 end
     
-function gdagpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{1,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]+imm(a.s[4]))
     r2 = g\(a.s[2]+imm(a.s[3]))
     return Spinor{4,G}((r1,r2,mimm(r2),mimm(r1)))
 end
-function gdagpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{1,-1}}, g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]-imm(a.s[4]))
     r2 = g\(a.s[2]-imm(a.s[3]))
     return Spinor{4,G}((r1,r2,imm(r2),imm(r1)))
 end
 
-function gdagpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{2,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]+a.s[4])
     r2 = g\(a.s[2]-a.s[3])
     return Spinor{4,G}((r1,r2,-r2,r1))
 end
-function gdagpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{2,-1}}, g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]-a.s[4])
     r2 = g\(a.s[2]+a.s[3])
     return Spinor{4,G}((r1,r2,r2,-r1))
 end
 
-function gdagpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{3,1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]+imm(a.s[3]))
     r2 = g\(a.s[2]-imm(a.s[4]))
     return Spinor{4,G}((r1,r2,mimm(r1),imm(r2)))
 end
-function gdagpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
+@inline function gdagpmul(::Type{Pgamma{3,-1}},  g, a::Spinor{4,G}) where {NS,G}
 
     r1 = g\(a.s[1]-imm(a.s[3]))
     r2 = g\(a.s[2]+imm(a.s[4]))
@@ -298,22 +298,22 @@ indexing for Dirac basis ``\\Gamma_n``:
 16  identity
 
 """
-dmul(::Type{Gamma{1}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[4]), mimm(a.s[3]), imm(a.s[2]), imm(a.s[1]))
+@inline dmul(::Type{Gamma{1}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[4]), mimm(a.s[3]), imm(a.s[2]), imm(a.s[1]))
-dmul(::Type{Gamma{2}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[4], a.s[3], a.s[2], -a.s[1])
+@inline dmul(::Type{Gamma{2}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[4], a.s[3], a.s[2], -a.s[1])
-dmul(::Type{Gamma{3}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[3]), imm(a.s[4]), imm(a.s[1]), mimm(a.s[2]))
+@inline dmul(::Type{Gamma{3}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[3]), imm(a.s[4]), imm(a.s[1]), mimm(a.s[2]))
-dmul(::Type{Gamma{4}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[3], -a.s[4], -a.s[1], -a.s[2])
+@inline dmul(::Type{Gamma{4}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[3], -a.s[4], -a.s[1], -a.s[2])
-dmul(::Type{Gamma{5}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1],  a.s[2],  -a.s[3], -a.s[4])
+@inline dmul(::Type{Gamma{5}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1],  a.s[2],  -a.s[3], -a.s[4])
-dmul(::Type{Gamma{6}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[4]),  imm(a.s[3]),   imm(a.s[2]),  imm(a.s[1]))
+@inline dmul(::Type{Gamma{6}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[4]),  imm(a.s[3]),   imm(a.s[2]),  imm(a.s[1]))
-dmul(::Type{Gamma{7}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[4], -a.s[3],   a.s[2], -a.s[1])
+@inline dmul(::Type{Gamma{7}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[4], -a.s[3],   a.s[2], -a.s[1])
-dmul(::Type{Gamma{8}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[3]), mimm(a.s[4]),   imm(a.s[1]), mimm(a.s[2]))
+@inline dmul(::Type{Gamma{8}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[3]), mimm(a.s[4]),   imm(a.s[1]), mimm(a.s[2]))
-dmul(::Type{Gamma{9}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[3],  a.s[4],  -a.s[1], -a.s[2])
+@inline dmul(::Type{Gamma{9}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[3],  a.s[4],  -a.s[1], -a.s[2])
-dmul(::Type{Gamma{10}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[2],  a.s[1],  -a.s[4], -a.s[3])
+@inline dmul(::Type{Gamma{10}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[2],  a.s[1],  -a.s[4], -a.s[3])
-dmul(::Type{Gamma{11}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[2]),  imm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
+@inline dmul(::Type{Gamma{11}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(mimm(a.s[2]),  imm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
-dmul(::Type{Gamma{12}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1], -a.s[2],  -a.s[3],  a.s[4])
+@inline dmul(::Type{Gamma{12}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( a.s[1], -a.s[2],  -a.s[3],  a.s[4])
-dmul(::Type{Gamma{13}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[1],  a.s[2],  -a.s[3],  a.s[4])
+@inline dmul(::Type{Gamma{13}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[1],  a.s[2],  -a.s[3],  a.s[4])
-dmul(::Type{Gamma{14}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[2], -a.s[1],  -a.s[4], -a.s[3])
+@inline dmul(::Type{Gamma{14}}, a::Spinor{4,G}) where {G} = Spinor{4,G}(-a.s[2], -a.s[1],  -a.s[4], -a.s[3])
-dmul(::Type{Gamma{15}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[2]), mimm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
+@inline dmul(::Type{Gamma{15}}, a::Spinor{4,G}) where {G} = Spinor{4,G}( imm(a.s[2]), mimm(a.s[1]),   imm(a.s[4]), mimm(a.s[3]))
-dmul(::Type{Gamma{16}}, a::Spinor{4,G}) where {G} = a
+@inline dmul(::Type{Gamma{16}}, a::Spinor{4,G}) where {G} = a
 
 export Spinor, Pgamma
 export norm, norm2, dot, imm, mimm

src/YM/YM.jl

@@ -12,7 +12,7 @@
 
 module YM
 
-using CUDA, Random, TimerOutputs
+using CUDA, Random, TimerOutputs, BDIO
 using ..Space
 using ..Groups
 using ..Fields
@@ -147,6 +147,6 @@ include("YMsf.jl")
 export sfcoupling, bndfield, setbndfield
 
 include("YMio.jl")
-export import_lex64, import_cern64
+export import_lex64, import_cern64, import_bsfqcd
 
 end

src/YM/YMio.jl

@@ -9,6 +9,49 @@
 ### created: Wed Nov 10 12:58:27 2021
 ###                               
 
+"""
+    function import_lex64(fname::String, lp::SpaceParm)
+
+import a double precision configuration in lexicographic format. SF boundary conditions are assummed. 
+"""
+function import_bsfqcd(fname, lp::SpaceParm)
+
+    fp = BDIO_open(fname, "r")
+    while BDIO_seek!(fp)
+        if (BDIO_get_uinfo(fb) == 2)
+            break
+        end
+    end
+
+    
+    dtr = [2,3,4,1]
+
+    assign(id, V, i3) = SU3{Float64}(V[1,dtr[id],i3],V[2,dtr[id],i3],V[3,dtr[id],i3],
+                                     V[4,dtr[id],i3],V[5,dtr[id],i3],V[6,dtr[id],i3])
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 3}(undef, 9, lp.ndim, lp.iL[3])
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                BDIO_read(fp, V)
+                for i3 in 1:lp.iL[3]
+                    b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                    for id in 1:lp.ndim
+                        Ucpu[b,id,r] = assign(id, V, i3)
+                    end
+                end
+            end
+        end
+    end
+
+    BDIO_read(fp, V)
+    Ubnd = ntuple(i->assign(i, V, 1), 3)
+    close(fp)
+
+    return CuArray(Ucpu), Ubnd
+end
+
 """
     function import_lex64(fname::String, lp::SpaceParm)