From ffc6818ba5786f758aaffa97fcfd7ac20450648e Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Sat, 13 Nov 2021 23:13:19 +0100
Subject: [PATCH] Added routine to import configuration in lex format

---
 src/LatticeGPU.jl  |  1 +
 src/Space/Space.jl | 13 +++++++++---
 src/YM/YM.jl       |  3 +++
 src/YM/YMio.jl     | 50 ++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 64 insertions(+), 3 deletions(-)
 create mode 100644 src/YM/YMio.jl

diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index 17d6277..b852c52 100644
--- a/src/LatticeGPU.jl
+++ b/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
diff --git a/src/Space/Space.jl b/src/Space/Space.jl
index 5ff9001..8fa67a7 100644
--- a/src/Space/Space.jl
+++ b/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
diff --git a/src/YM/YM.jl b/src/YM/YM.jl
index 3b2d10a..547e0dc 100644
--- a/src/YM/YM.jl
+++ b/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
diff --git a/src/YM/YMio.jl b/src/YM/YMio.jl
new file mode 100644
index 0000000..923115f
--- /dev/null
+++ b/src/YM/YMio.jl
@@ -0,0 +1,50 @@
+###
+### "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
-- 
GitLab