So I understand that for the write there’s nothing to do.
regarding my struggles with the vector notation (I’m not sure that I used the term “broadcasting” in the correct way isn’t it the . notation? like in .= .+= abs.( ) ?) in the function below you can see a lot of short loops on ii that were made in order to avoid the allocations I was getting with vector notation, I also tried to use StaticArrays that worked in some cases. The whole experience has been, to me “a nightmare” and that’s makes me angry as i really like a lot julia.
(Btw that way I was able (mostly by clever algorithm in another function) to reach 2.4x speed over C)
function read_xtc_atoms(file::XtcFile, frame::Integer, coords::AbstractMatrix{T}) where {T <: Real}
local smallnum::Int
local smaller::Int
minint = MVector{3, Int32}(undef)
maxint = MVector{3, Int32}(undef)
thiscoord = MVector{3, Int}(undef)
prevcoord = MVector{3, Int}(undef)
# magicints = file.magicints
FIRSTINDEX = 9
seek(file.file, file.offsets[frame] + 36) # we don't read box here
size = ntoh(read(file.file, Int32))
if size <= 9
for j in 1:3
for i in 1:size
coords[i, j] = ntoh(read(file.file, Float32))
end
end
else
precision = ntoh(read(file.file, Float32))
read!(file.file, minint)
minint .= ntoh.(minint)
read!(file.file, maxint)
maxint .= ntoh.(maxint)
small_idx = ntoh(read(file.file, Int32))
smallidx = Int(small_idx)
smaller = magicints[max(FIRSTINDEX, smallidx - 1) + 1] Ă· 2
smallnum = magicints[smallidx + 1] Ă· 2
sizesmall = SA[magicints[smallidx + 1], magicints[smallidx + 1], magicints[smallidx + 1]]
nbytes = ntoh(read(file.file, Int32))
buffer = BitBuffer(nbytes)
readbytes!(file.file, buffer.bits, nbytes)
sizeint = SA[maxint[1] - minint[1] + 1,
maxint[2] - minint[2] + 1,
maxint[3] - minint[3] + 1]
if any( sizeint .> 2^24-1)
bitsizeint = sizeofint.(sizeint)
bitsize = 0
else
bitsize = sizeofints(sizeint)
end
i = 1 # atom index
run = 0
while i <= file.natoms
if bitsize == 0
for ii in 1:3
thiscoord[ii] = receivebits!(buffer, bitsizeint[ii])
end
else
receiveints!(buffer, bitsize, sizeint, thiscoord)
end
for ii in 1:3
thiscoord[ii] += minint[ii]
end
flag = receivebits!(buffer, 1)
is_smaller = 0
if flag == 1
run = receivebits!(buffer, 5)
is_smaller = run % 3
run -= is_smaller
is_smaller -= 1
end
if run == 0
for ii in 1:3
coords[ii, i] = thiscoord[ii] / precision
end
i += 1
else
for ii in 1:3
prevcoord[ii] = thiscoord[ii]
end
for k in 1:3:run
receiveints!(buffer, smallidx, sizesmall, thiscoord)
for ii in 1:3
thiscoord[ii] += prevcoord[ii] - smallnum # smallnum ??? WTF!!
end
if k == 1 # exchange first with second atom WTF!!
thiscoord, prevcoord = prevcoord, thiscoord
for ii in 1:3
coords[ii,i] = prevcoord[ii] / precision
end
i += 1
else
for ii in 1:3
prevcoord[ii] = thiscoord[ii]
end
end
for ii in 1:3
coords[ii,i] = thiscoord[ii] / precision
end
i += 1
end
end
smallidx += is_smaller;
if is_smaller < 0
smallnum = smaller
if smallidx > FIRSTINDEX
smaller = magicints[smallidx] >>> 1
else
smaller = 0
end
elseif is_smaller > 0
smaller = smallnum
smallnum = magicints[smallidx + 1] >>> 1
end
sizesmall = SA[magicints[smallidx + 1], magicints[smallidx + 1], magicints[smallidx + 1]]
end
end
return nothing
end