Added routine to import configuration in lex format

src/LatticeGPU.jl

@@ -45,5 +45,6 @@ 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
 
 end # module

src/Space/Space.jl

@@ -216,7 +216,7 @@ Given a point `x` with index `p`, this routine returns the index of the points
         ic = mod(div(p[2]-1,lp.rbkS[id]),lp.rbk[id])
         if (ic == lp.rbk[id]-1)
             ru = p[2] - (lp.rbk[id]-1)*lp.rbkS[id]
-        else
+sfqcd_L12_b3.5320_k0.137101170000000_r0_id1n1        else
             ru = p[2] + lp.rbkS[id]
         end
         rd = p[2]
@@ -306,9 +306,16 @@ end
 end
 
 
-@inline function point_index(pt::NTuple{4, Int64}, lp::SpaceParm)
+@inline function point_index(pt::CartesianIndex, lp::SpaceParm)
 
-    
+    b = 1
+    r = 1
+    for i in 1:length(pt)
+        b = b + ((pt[i]-1)%lp.blk[i])*lp.blkS[i]
+        r = r + div((pt[i]-1),lp.blk[i])*lp.rbkS[i]
+    end
+
+    return (b,r)
 end
     
 export up, dw, updw, point_index, point_coord, point_time

src/YM/YM.jl

@@ -145,4 +145,7 @@ export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
 include("YMsf.jl")
 export sfcoupling, bndfield, setbndfield
 
+include("YMio.jl")
+export import_lex64
+
 end

src/YM/YMio.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    YMio.jl
+### 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_lex64(fname, lp::SpaceParm)
+
+    fp = open(fname, "r")
+
+    dtr = [4,1,2,3]
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    Ubnd = Array{SU3{Float64}, 1}(undef, lp.ndim-1)
+    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]
+                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,dtr[id],r] = SU3{Float64}(V[1,id,i4],V[2,id,i4],V[3,id,i4],
+                                                         V[4,id,i4],V[5,id,i4],V[6,id,i4])
+                    end
+                end
+            end
+        end
+    end
+
+    read!(fp, V)
+    for id in 2:lp.ndim
+        Ubnd[dtr[id]] = SU3{Float64}(V[1,id,1],V[2,id,1],V[3,id,1],
+                                     V[4,id,1],V[5,id,1],V[6,id,1])
+    end
+
+    close(fp)
+
+    return CuArray(Ucpu), Ubnd
+end