Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Contains routines for data compression: packing integers of Nbits bits each into a
10 : !> stream of 64-bit words, and unpacking them. Per number of bits, fypp generates a
11 : !> routine that assembles each packed word in a register from constant shifts and stores
12 : !> it once, for every group of 64 integers (which fills exactly Nbits words). Moreover,
13 : !> packing at a width per group of integers (the width the group needs), in 8 streams
14 : !> whose shifts are uniform over the lanes.
15 : !> \par History
16 : !> 04.2008 created [Manuel Guidon]
17 : !> 10.2026 generated by fypp instead of a generated source file [Hans Pabst]
18 : !> \author Manuel Guidon
19 : ! **************************************************************************************************
20 : MODULE hfx_compression_core_methods
21 :
22 : USE kinds, ONLY: int_8
23 : #include "../base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 : PRIVATE
27 :
28 : PUBLIC :: bits2ints_specific, groups2ints, ints2bits_specific, ints2groups
29 :
30 : CONTAINS
31 :
32 : ! **************************************************************************************************
33 : !> \brief Packs integers into a stream of 64-bit words, Nbits bits per integer.
34 : !> \param Nbits number of bits per integer (1 to 64); bits beyond them are ignored
35 : !> \param Ndata number of integers
36 : !> \param packed_data the stream of bits, CEILING(Ndata*Nbits/64) words
37 : !> \param full_data the integers
38 : ! **************************************************************************************************
39 0 : SUBROUTINE ints2bits_specific(Nbits, Ndata, packed_data, full_data)
40 : INTEGER, INTENT(IN) :: Nbits, Ndata
41 : INTEGER(KIND=int_8), INTENT(OUT) :: packed_data(*)
42 : INTEGER(KIND=int_8), INTENT(IN) :: full_data(*)
43 :
44 0 : SELECT CASE (Nbits)
45 : #:for nb in range(1, 65)
46 : CASE (${nb}$)
47 0 : CALL ints2bits_${nb}$ (Ndata, packed_data, full_data)
48 : #:endfor
49 : CASE DEFAULT
50 0 : CPABORT("Invalid number of bits")
51 : END SELECT
52 0 : END SUBROUTINE ints2bits_specific
53 :
54 : ! **************************************************************************************************
55 : !> \brief Unpacks integers from a stream of 64-bit words, Nbits bits per integer.
56 : !> \param Nbits number of bits per integer (1 to 64)
57 : !> \param Ndata number of integers
58 : !> \param packed_data the stream of bits, CEILING(Ndata*Nbits/64) words
59 : !> \param full_data the integers
60 : ! **************************************************************************************************
61 0 : SUBROUTINE bits2ints_specific(Nbits, Ndata, packed_data, full_data)
62 : INTEGER, INTENT(IN) :: Nbits, Ndata
63 : INTEGER(KIND=int_8), INTENT(IN) :: packed_data(*)
64 : INTEGER(KIND=int_8), INTENT(OUT) :: full_data(*)
65 :
66 0 : SELECT CASE (Nbits)
67 : #:for nb in range(1, 65)
68 : CASE (${nb}$)
69 0 : CALL bits2ints_${nb}$ (Ndata, packed_data, full_data)
70 : #:endfor
71 : CASE DEFAULT
72 0 : CPABORT("Invalid number of bits")
73 : END SELECT
74 0 : END SUBROUTINE bits2ints_specific
75 :
76 : #:for nb in range(1, 65)
77 : #:set mask = (1 << nb) - 1 if nb < 64 else -1
78 : ! **************************************************************************************************
79 : !> \brief Packs integers of ${nb}$ bits each (see ints2bits_specific).
80 : !> \param Ndata ...
81 : !> \param packed_data ...
82 : !> \param full_data ...
83 : ! **************************************************************************************************
84 0 : SUBROUTINE ints2bits_${nb}$ (Ndata, packed_data, full_data)
85 : INTEGER, INTENT(IN) :: Ndata
86 : INTEGER(KIND=int_8), INTENT(OUT) :: packed_data(*)
87 : INTEGER(KIND=int_8), INTENT(IN) :: full_data(*)
88 :
89 : INTEGER(KIND=int_8), PARAMETER :: mask = ${mask}$_int_8
90 :
91 : INTEGER :: igroup, ipack, idata
92 : INTEGER(KIND=int_8) :: pack_tmp
93 :
94 0 : DO igroup = 0, Ndata/64 - 1
95 0 : idata = 64*igroup
96 0 : ipack = ${nb}$*igroup
97 : #:for w in range(nb)
98 : #:set js = [j for j in range(64) if j*nb < 64*w + 64 and j*nb + nb > 64*w]
99 0 : pack_tmp = ISHFT(IAND(full_data(idata + ${js[0] + 1}$), mask), ${js[0]*nb - 64*w}$)
100 : #:for j in js[1:]
101 0 : pack_tmp = IOR(pack_tmp, ISHFT(IAND(full_data(idata + ${j + 1}$), mask), ${j*nb - 64*w}$))
102 : #:endfor
103 0 : packed_data(ipack + ${w + 1}$) = pack_tmp
104 : #:endfor
105 : END DO
106 0 : CALL ints2bits_tail(${nb}$, Ndata, packed_data, full_data)
107 0 : END SUBROUTINE ints2bits_${nb}$
108 :
109 : ! **************************************************************************************************
110 : !> \brief Unpacks integers of ${nb}$ bits each (see bits2ints_specific).
111 : !> \param Ndata ...
112 : !> \param packed_data ...
113 : !> \param full_data ...
114 : ! **************************************************************************************************
115 0 : SUBROUTINE bits2ints_${nb}$ (Ndata, packed_data, full_data)
116 : INTEGER, INTENT(IN) :: Ndata
117 : INTEGER(KIND=int_8), INTENT(IN) :: packed_data(*)
118 : INTEGER(KIND=int_8), INTENT(OUT) :: full_data(*)
119 :
120 : INTEGER(KIND=int_8), PARAMETER :: mask = ${mask}$_int_8
121 :
122 : INTEGER :: igroup, ipack, idata
123 :
124 0 : DO igroup = 0, Ndata/64 - 1
125 0 : idata = 64*igroup
126 0 : ipack = ${nb}$*igroup
127 : #:for j in range(64)
128 : #:set w = (j*nb)//64
129 : #:set off = (j*nb) % 64
130 : #:if off + nb > 64
131 : full_data(idata + ${j + 1}$) = IAND(IOR(ISHFT(packed_data(ipack + ${w + 1}$), ${-off}$), &
132 0 : ISHFT(packed_data(ipack + ${w + 2}$), ${64 - off}$)), mask)
133 : #:else
134 0 : full_data(idata + ${j + 1}$) = IAND(ISHFT(packed_data(ipack + ${w + 1}$), ${-off}$), mask)
135 : #:endif
136 : #:endfor
137 : END DO
138 0 : CALL bits2ints_tail(${nb}$, Ndata, packed_data, full_data)
139 0 : END SUBROUTINE bits2ints_${nb}$
140 :
141 : #:endfor
142 : ! **************************************************************************************************
143 : !> \brief Packs the integers after the last whole group of 64 (fewer than 64 integers).
144 : !> \param Nbits number of bits per integer
145 : !> \param Ndata number of integers including the whole groups
146 : !> \param packed_data the stream of bits
147 : !> \param full_data the integers
148 : ! **************************************************************************************************
149 0 : SUBROUTINE ints2bits_tail(Nbits, Ndata, packed_data, full_data)
150 : INTEGER, INTENT(IN) :: Nbits, Ndata
151 : INTEGER(KIND=int_8), INTENT(INOUT) :: packed_data(*)
152 : INTEGER(KIND=int_8), INTENT(IN) :: full_data(*)
153 :
154 : INTEGER :: ibit, idata, ipack, ioff, ndone
155 : INTEGER(KIND=int_8) :: mask, data_tmp
156 :
157 0 : ndone = (Ndata/64)*64
158 0 : IF (ndone < Ndata) THEN
159 0 : mask = ISHFT(1_int_8, Nbits) - 1 ! all bits if Nbits is 64
160 0 : ipack = (ndone/64)*Nbits
161 0 : packed_data(ipack + 1:ipack + ((Ndata - ndone)*Nbits + 63)/64) = 0
162 0 : ibit = 0
163 0 : DO idata = ndone + 1, Ndata
164 0 : data_tmp = IAND(full_data(idata), mask)
165 0 : ioff = MOD(ibit, 64)
166 0 : packed_data(ipack + ibit/64 + 1) = IOR(packed_data(ipack + ibit/64 + 1), ISHFT(data_tmp, ioff))
167 0 : IF (ioff + Nbits > 64) packed_data(ipack + ibit/64 + 2) = ISHFT(data_tmp, ioff - 64)
168 0 : ibit = ibit + Nbits
169 : END DO
170 : END IF
171 0 : END SUBROUTINE ints2bits_tail
172 :
173 : ! **************************************************************************************************
174 : !> \brief Unpacks the integers after the last whole group of 64 (fewer than 64 integers).
175 : !> \param Nbits number of bits per integer
176 : !> \param Ndata number of integers including the whole groups
177 : !> \param packed_data the stream of bits
178 : !> \param full_data the integers
179 : ! **************************************************************************************************
180 0 : SUBROUTINE bits2ints_tail(Nbits, Ndata, packed_data, full_data)
181 : INTEGER, INTENT(IN) :: Nbits, Ndata
182 : INTEGER(KIND=int_8), INTENT(IN) :: packed_data(*)
183 : INTEGER(KIND=int_8), INTENT(INOUT) :: full_data(*)
184 :
185 : INTEGER :: ibit, idata, ipack, ioff, ndone
186 : INTEGER(KIND=int_8) :: mask, data_tmp
187 :
188 0 : ndone = (Ndata/64)*64
189 0 : IF (ndone < Ndata) THEN
190 0 : mask = ISHFT(1_int_8, Nbits) - 1 ! all bits if Nbits is 64
191 0 : ipack = (ndone/64)*Nbits
192 0 : ibit = 0
193 0 : DO idata = ndone + 1, Ndata
194 0 : ioff = MOD(ibit, 64)
195 0 : data_tmp = ISHFT(packed_data(ipack + ibit/64 + 1), -ioff)
196 0 : IF (ioff + Nbits > 64) data_tmp = IOR(data_tmp, ISHFT(packed_data(ipack + ibit/64 + 2), 64 - ioff))
197 0 : full_data(idata) = IAND(data_tmp, mask)
198 0 : ibit = ibit + Nbits
199 : END DO
200 : END IF
201 0 : END SUBROUTINE bits2ints_tail
202 :
203 : ! **************************************************************************************************
204 : !> \brief Packs integers at a width per group, the width the group's largest magnitude needs:
205 : !> header bytes (the group size, then a width per group), then 8 interleaved streams
206 : !> (integer 8*c+l in stream l, all at the same bit position), so that every shift is
207 : !> uniform over the 8 lanes; a group is Ngroup/8 such columns at its width. The integers
208 : !> are signed, of Nbits bits (sign folded).
209 : !> \param Nbits number of bits per integer (1 to 64); bits beyond them are ignored
210 : !> \param Ngroup number of integers per group, a multiple of 8 (8 to 248)
211 : !> \param packed_data the stream: words for 1+CEILING(n/Ngroup) header bytes, then at most
212 : !> 8*CEILING(((n+7)/8)*Nbits/64) words (n integers)
213 : !> \param full_data the integers
214 : !> \return number of words written
215 : ! **************************************************************************************************
216 0 : FUNCTION ints2groups(Nbits, Ngroup, packed_data, full_data) RESULT(nwords)
217 : INTEGER, INTENT(IN) :: Nbits, Ngroup
218 : INTEGER(KIND=int_8), CONTIGUOUS, INTENT(OUT) :: packed_data(:)
219 : INTEGER(KIND=int_8), CONTIGUOUS, INTENT(IN) :: full_data(:)
220 : INTEGER :: nwords
221 :
222 : INTEGER :: b, c, i, j, lo, m, ncols, Ndata, &
223 : ngroups, nhead, nmax, off, s, w
224 : INTEGER(KIND=int_8) :: z(256)
225 :
226 0 : CPASSERT(0 < Ngroup .AND. Ngroup <= 248 .AND. 0 == MOD(Ngroup, 8))
227 0 : Ndata = SIZE(full_data)
228 0 : ngroups = (Ndata + Ngroup - 1)/Ngroup
229 0 : nhead = (1 + ngroups + 7)/8
230 0 : nmax = nhead + 8*((((Ndata + 7)/8)*Nbits + 63)/64) ! words at most
231 0 : CPASSERT(nmax <= SIZE(packed_data))
232 0 : packed_data(1:nhead) = 0
233 0 : packed_data(1) = INT(Ngroup, KIND=int_8) ! group size
234 0 : s = 64 - Nbits
235 0 : j = 0
236 0 : off = 0
237 0 : DO i = 0, ngroups - 1
238 0 : m = MIN(Ngroup, Ndata - Ngroup*i)
239 0 : ncols = (m + 7)/8
240 0 : z(m + 1:8*ncols) = 0
241 : ! sign-extended from Nbits, then the sign folded (zig-zag)
242 0 : z(1:m) = SHIFTA(ISHFT(full_data(Ngroup*i + 1:Ngroup*i + m), s), s)
243 0 : z(1:m) = IEOR(ISHFT(z(1:m), 1), SHIFTA(z(1:m), 63))
244 0 : w = 64 - LEADZ(IANY(z(1:m)))
245 0 : b = i + 1 ! header byte of the width
246 0 : packed_data(b/8 + 1) = IOR(packed_data(b/8 + 1), ISHFT(INT(w, KIND=int_8), 8*MOD(b, 8)))
247 0 : IF (w > 0) THEN
248 0 : DO c = 0, ncols - 1
249 0 : lo = nhead + 8*j
250 0 : IF (off /= 0) THEN
251 0 : packed_data(lo + 1:lo + 8) = IOR(packed_data(lo + 1:lo + 8), ISHFT(z(8*c + 1:8*c + 8), off))
252 : ELSE
253 0 : packed_data(lo + 1:lo + 8) = z(8*c + 1:8*c + 8)
254 : END IF
255 0 : IF (off + w > 64) THEN ! spills into the next word
256 0 : packed_data(lo + 9:lo + 16) = ISHFT(z(8*c + 1:8*c + 8), off - 64)
257 : END IF
258 0 : off = off + w
259 0 : IF (off >= 64) THEN
260 0 : off = off - 64
261 0 : j = j + 1
262 : END IF
263 : END DO
264 : END IF
265 : END DO
266 0 : nwords = nhead + 8*j
267 0 : IF (off /= 0) nwords = nwords + 8 ! partial last word
268 0 : END FUNCTION ints2groups
269 :
270 : ! **************************************************************************************************
271 : !> \brief Unpacks integers that ints2groups packed (signed, the sign unfolded); the group size
272 : !> is read from the stream.
273 : !> \param packed_data the stream
274 : !> \param full_data the integers
275 : !> \return number of words read
276 : ! **************************************************************************************************
277 0 : FUNCTION groups2ints(packed_data, full_data) RESULT(nwords)
278 : INTEGER(KIND=int_8), CONTIGUOUS, INTENT(IN) :: packed_data(:)
279 : INTEGER(KIND=int_8), CONTIGUOUS, INTENT(OUT) :: full_data(:)
280 : INTEGER :: nwords
281 :
282 : INTEGER :: b, base, c, i, j, lo, m, mc, Ndata, &
283 : ngroup, ngroups, nhead, off, w
284 : INTEGER(KIND=int_8) :: mask, z(8)
285 :
286 0 : ngroup = INT(IAND(packed_data(1), 255_int_8))
287 0 : CPASSERT(0 < ngroup .AND. 0 == MOD(ngroup, 8))
288 0 : Ndata = SIZE(full_data)
289 0 : ngroups = (Ndata + ngroup - 1)/ngroup
290 0 : nhead = (1 + ngroups + 7)/8
291 0 : j = 0
292 0 : off = 0
293 0 : DO i = 0, ngroups - 1
294 0 : m = MIN(ngroup, Ndata - ngroup*i)
295 0 : b = i + 1 ! header byte of the width
296 0 : w = INT(IAND(ISHFT(packed_data(b/8 + 1), -8*MOD(b, 8)), 255_int_8))
297 0 : mask = ISHFT(-1_int_8, w - 64)
298 0 : DO c = 0, (m + 7)/8 - 1
299 0 : base = ngroup*i + 8*c
300 0 : mc = MIN(8, m - 8*c)
301 0 : IF (w > 0) THEN
302 0 : lo = nhead + 8*j
303 0 : IF (off + w > 64) THEN ! spills into the next word
304 0 : z = IAND(IOR(ISHFT(packed_data(lo + 1:lo + 8), -off), ISHFT(packed_data(lo + 9:lo + 16), 64 - off)), mask)
305 : ELSE
306 0 : z = IAND(ISHFT(packed_data(lo + 1:lo + 8), -off), mask)
307 : END IF
308 0 : off = off + w
309 0 : IF (off >= 64) THEN
310 0 : off = off - 64
311 0 : j = j + 1
312 : END IF
313 : ELSE
314 0 : z = 0
315 : END IF
316 0 : full_data(base + 1:base + mc) = IEOR(ISHFT(z(1:mc), -1), -IAND(z(1:mc), 1_int_8)) ! the sign unfolded
317 : END DO
318 : END DO
319 0 : nwords = nhead + 8*j
320 0 : IF (off /= 0) nwords = nwords + 8 ! partial last word
321 0 : END FUNCTION groups2ints
322 :
323 : END MODULE hfx_compression_core_methods
|