1!
2!  Dalton, a molecular electronic structure program
3!  Copyright (C) by the authors of Dalton.
4!
5!  This program is free software; you can redistribute it and/or
6!  modify it under the terms of the GNU Lesser General Public
7!  License version 2.1 as published by the Free Software Foundation.
8!
9!  This program is distributed in the hope that it will be useful,
10!  but WITHOUT ANY WARRANTY; without even the implied warranty of
11!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12!  Lesser General Public License for more details.
13!
14!  If a copy of the GNU LGPL v2.1 was not distributed with this
15!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
16!
17!
18#if (defined (SYS_AIX) && !defined (VAR_INT64))
19#define BIGENDIAN 1
20#else
21#define LITTLEENDIAN 1
22#endif
23C
24C  /* Deck pckr8 */
25      SUBROUTINE PCKR8(BUFFER,NDATA,BUFPCK,NBYTE,IPCKTABLE,LPACK)
26*=====================================================================*
27C
28C     Purpose: compress BUFFER with NDATA real*8 numbers
29C
30C        BUFFER -- array with unpacked real*8 numbers
31C        BUFPCK -- buffer with the compressed numbers
32C        NDATA  -- length of unpacked array BUFFER in R*8 words
33C        NBYTE  -- length of packed array BUFPCK in bytes
34C        IPCKTABLE -- packing table generated by INITPCKR8
35C        LPACK  -- flag, if set to .false. PCKR8 does only a DCOPY
36C
37C
38C     Christof Haettig, August 1998
39C
40*---------------------------------------------------------------------*
41      IMPLICIT NONE
42C
43#include "priunit.h"
44C
45      LOGICAL LOCDBG
46      PARAMETER (LOCDBG = .FALSE.)
47C
48      LOGICAL LPACK
49      INTEGER NBYTE, NDATA
50      INTEGER IPCKTABLE(0:255)
51C
52      REAL*8 BUFFER(NDATA)
53      REAL*8 R8WORD, ZERO
54C
55      PARAMETER (ZERO = 0.0D0)
56C
57      LOGICAL*1 BUFPCK(*)
58C
59      LOGICAL*1 L1BYTE(8), LTAB(4)
60      INTEGER*4 I4WORD(2), ITAB
61C
62      EQUIVALENCE(R8WORD,I4WORD(1))
63      EQUIVALENCE(R8WORD,L1BYTE(1))
64      EQUIVALENCE(ITAB,LTAB(1))
65C
66      INTEGER IB1, IW1, IW2, IT1, L, I, OFF, MBYTE
67      INTEGER*4 M1, MASK1, MASK2, HBITS, ISIGN, ITMP
68C
69#ifdef LITTLEENDIAN
70      PARAMETER (IT1 = 1)
71      PARAMETER (IB1 = 8)
72      PARAMETER (IW1 = 2, IW2 = 1)
73#elif BIGENDIAN
74      PARAMETER (IT1 = 4)
75      PARAMETER (IB1 = 1)
76      PARAMETER (IW1 = 1, IW2 = 2)
77#endif
78C
79#if defined (VAR_GFORTRAN)
80      ! special code to avoid gfortran out of range error
81      PARAMETER (M1 = z'40000000')
82      PARAMETER (MASK1 = -M1)
83#else
84      PARAMETER (MASK1 = z'C0000000')
85#endif
86      PARAMETER (MASK2 = z'3FFFFFFF')
87
88C
89C---------------------------------------------------------------
90C     if packing disabled, copy input data to output array
91C----------------------------------------------------------------
92C
93      IF (.NOT.LPACK) THEN
94        CALL DCOPY(NDATA,BUFFER,1,BUFPCK,1)
95        NBYTE = 8 * NDATA
96        RETURN
97      END IF
98C
99C----------------------------------------------------------------
100C     Pack the input data and copy to output array
101C----------------------------------------------------------------
102C
103#if !(defined (BIGENDIAN) || defined(LITTLEENDIAN))
104      CALL QUIT('PCKR8 not implemented for this operating system.')
105#endif
106      NBYTE = 0
107      ITAB  = 0
108C
109      DO I = 1, NDATA
110C
111         IF (LOCDBG) THEN
112           WRITE (LUPRI,'(//A,I5,A,E30.20)') 'Element #',I,' = ',
113     &           BUFFER(I)
114           CALL R8BITREP(BUFFER(I))
115         END IF
116C
117         R8WORD      = BUFFER(I)
118C
119         ISIGN       = IAND( I4WORD(IW1),MASK1)
120         I4WORD(IW1) = ISHFT(I4WORD(IW1),4)
121         I4WORD(IW1) = IAND( I4WORD(IW1),MASK2)
122         I4WORD(IW1) = IOR(  I4WORD(IW1),ISIGN)
123C
124         HBITS       = ISHFT(I4WORD(IW2),-28)
125         I4WORD(IW2) = ISHFT(I4WORD(IW2),4)
126         I4WORD(IW1) = IOR(  I4WORD(IW1),HBITS)
127C
128         LTAB(IT1)   = L1BYTE(IB1)
129C
130         IF (LOCDBG) CALL I4BITREP(ITAB)
131C
132         IF (LOCDBG) CALL R8BITREP(R8WORD)
133C
134         MBYTE = IPCKTABLE(ITAB)
135C
136#ifdef LITTLEENDIAN
137         OFF = 8 - MBYTE
138         DO L = 1, MBYTE
139           BUFPCK(NBYTE+L) = L1BYTE(OFF+L)
140         END DO
141#elif defined(BIGENDIAN)
142         OFF = MBYTE + 1
143         DO L = 1, MBYTE
144           BUFPCK(NBYTE+L) = L1BYTE(OFF-L)
145         END DO
146#endif
147C
148         NBYTE = NBYTE + MBYTE
149C
150         IF (LOCDBG) THEN
151           WRITE (LUPRI,'(3I5)') ITAB, NBYTE, MBYTE
152         END IF
153C
154      END DO
155C
156C     RETURN
157      END
158*=====================================================================*
159C  /* Deck unpckr8 */
160      SUBROUTINE UNPCKR8(BUFFER,NDATA,BUFPCK,NBYTE,IPCKTABLE,LPACK)
161*=====================================================================*
162C
163C     Purpose: unpack BUFPCK with NBYTE compressed data
164C              containing NDATA real*8 numbers
165C
166C        BUFFER -- array with unpacked real*8 numbers
167C        BUFPCK -- buffer with the compressed numbers
168C        NDATA  -- length of unpacked array BUFUPK in R*8 words
169C        NBYTE  -- length of packed array BUFPCK in bytes
170C        IPCKTABLE -- packing table generated by INITPCKR8
171C        LPACK  -- flag, if set to .false. PCKR8 does only a DCOPY
172C
173C     note that both dimensions NBYTE and NDATA are input!
174C
175C     Christof Haettig, August 1998
176C
177*---------------------------------------------------------------------*
178      IMPLICIT NONE
179C
180#include "priunit.h"
181C
182      LOGICAL LOCDBG
183      PARAMETER (LOCDBG = .FALSE.)
184C
185      LOGICAL LPACK
186      INTEGER NBYTE, NDATA
187      INTEGER IPCKTABLE(0:255)
188C
189#if defined (SYS_CRAY)
190      REAL BUFFER(NDATA)
191      REAL ZERO
192#else
193      DOUBLE PRECISION BUFFER(NDATA)
194      DOUBLE PRECISION R8WORD, ZERO
195#endif
196C
197      PARAMETER (ZERO = 0.0D0)
198C
199      LOGICAL*1 BUFPCK(NBYTE)
200C
201      LOGICAL*1 L1BYTE(8), LTAB(4)
202      INTEGER*4 I4WORD(2), ITAB
203C
204      EQUIVALENCE(R8WORD,I4WORD(1))
205      EQUIVALENCE(R8WORD,L1BYTE(1))
206      EQUIVALENCE(ITAB,LTAB(1))
207C
208      INTEGER IB1, IB4, IW1, IW2, IT1, L, OFF, KBYTE, MBYTE, MDATA
209      INTEGER*4 M1, M4, MASK1, MASK3, MASK4, MASK5, HBITS, ISIGN, ITMP
210C
211#ifdef LITTLEENDIAN
212      PARAMETER (IT1 = 1)
213      PARAMETER (IB1 = 8)
214      PARAMETER (IW1 = 2, IW2 = 1)
215#elif defined(BIGENDIAN)
216      PARAMETER (IT1 = 4)
217      PARAMETER (IB1 = 1)
218      PARAMETER (IW1 = 1, IW2 = 2)
219#endif
220C
221#if defined (VAR_GFORTRAN)
222      ! special code to avoid gfortran out of range error
223      PARAMETER (M1    = z'40000000')
224      PARAMETER (MASK1 = -M1)
225      PARAMETER (MASK3 = z'3C000000')
226      PARAMETER (M4    = z'3C000001')
227      PARAMETER (MASK4 = -M4)
228      PARAMETER (MASK5 = z'0000000F')
229#else
230      PARAMETER (MASK1 = z'C0000000')
231      PARAMETER (MASK3 = z'3C000000')
232      PARAMETER (MASK4 = z'C3FFFFFF')
233      PARAMETER (MASK5 = z'0000000F')
234#endif
235C
236C---------------------------------------------------------------
237C     if packing disabled, copy input data to output array
238C----------------------------------------------------------------
239C
240      IF (.NOT.LPACK) THEN
241        CALL DCOPY(NDATA,BUFPCK,1,BUFFER,1)
242        RETURN
243      END IF
244C
245C---------------------------------------------------------------
246C     Unpack the input data and copy to output array
247C----------------------------------------------------------------
248C
249#if !(defined (BIGENDIAN) || defined(LITTLEENDIAN))
250      CALL QUIT('PCKR8 not implemented for this operating system.')
251#endif
252      MBYTE  = NBYTE
253      MDATA  = NDATA
254      ITAB   = 0
255      R8WORD = ZERO
256C
257      DO WHILE ( MDATA .GT. 0 )
258C
259         LTAB(IT1)   = BUFPCK(MBYTE)
260C
261         IF (LOCDBG) CALL I4BITREP(ITAB)
262C
263         KBYTE = IPCKTABLE(ITAB)
264C
265#ifdef LITTLEENDIAN
266         MBYTE = MBYTE - KBYTE
267         OFF = 8 - KBYTE
268         DO L = 1, KBYTE
269           L1BYTE(OFF+L) = BUFPCK(MBYTE+L)
270         END DO
271#elif defined(BIGENDIAN)
272         DO L = 1, KBYTE
273           L1BYTE(L) = BUFPCK(MBYTE+1-L)
274         END DO
275         MBYTE = MBYTE - KBYTE
276#endif
277C
278         IF (LOCDBG) CALL R8BITREP(R8WORD)
279C
280         HBITS = IAND(I4WORD(IW1),MASK5)
281         HBITS = ISHFT(HBITS,28)
282         I4WORD(IW2) = ISHFT(I4WORD(IW2),-4)
283         I4WORD(IW2) = IOR(I4WORD(IW2),HBITS)
284C
285         ISIGN       = IAND(I4WORD(IW1),MASK1)
286         I4WORD(IW1) = ISHFT(I4WORD(IW1),-4)
287         I4WORD(IW1) = IOR( I4WORD(IW1),ISIGN)
288C
289         IF ( BTEST(ITAB,6) ) THEN
290           I4WORD(IW1) = IAND(I4WORD(IW1),MASK4)
291         ELSE
292           I4WORD(IW1) = IOR( I4WORD(IW1),MASK3)
293         END IF
294C
295         BUFFER(MDATA)   = R8WORD
296C
297         IF (LOCDBG) THEN
298           WRITE (LUPRI,'(//A,I5,A,E30.20)') 'Element #',MDATA,' = ',
299     &           R8WORD
300           WRITE (LUPRI,'(3I5)') ITAB, MBYTE+KBYTE, KBYTE
301         END IF
302C
303         MDATA = MDATA - 1
304C
305      END DO
306C
307C     RETURN
308      END
309*=====================================================================*
310C  /* Deck initpckr8 */
311      SUBROUTINE INITPCKR8(THRESHOLD,IPCKTABLE,ENABLED)
312*=====================================================================*
313C
314C     Purpose: initialize exponent table for packing/unpacking
315C              of arrays with real*8 (double precision) numbers
316C
317C     Christof Haettig, August 1998
318C
319*---------------------------------------------------------------------*
320      IMPLICIT NONE
321C
322#include "priunit.h"
323C
324      LOGICAL LOCDBG
325      PARAMETER (LOCDBG = .FALSE.)
326C
327      INTEGER IPCKTABLE(0:255)
328      LOGICAL ENABLED
329C
330#if defined (SYS_CRAY)
331      REAL THRESHOLD, ZERO, ERRMAX(8)
332      REAL R8WORD0, R8WORD1, PCKTHRS
333#else
334      DOUBLE PRECISION THRESHOLD, ZERO, ERRMAX(8)
335      DOUBLE PRECISION R8WORD0, R8WORD1, PCKTHRS
336#endif
337C
338      PARAMETER (ZERO = 0.0D0)
339C
340      INTEGER IEXP, IWORD
341C
342      LOGICAL*1 L1BYTE0(8), L1BYTE1(8)
343      INTEGER*4 I4WORD0(2), I4WORD1(2)
344      INTEGER*4 IMASK(2,7), JMASK
345C
346      EQUIVALENCE(R8WORD0,I4WORD0(1))
347      EQUIVALENCE(R8WORD0,L1BYTE0(1))
348      EQUIVALENCE(R8WORD1,I4WORD1(1))
349      EQUIVALENCE(R8WORD1,L1BYTE1(1))
350C
351      INTEGER IW1, IW2, L
352C
353#undef NOPACKING
354#if LITTLEENDIAN
355      PARAMETER (IW1 = 2, IW2 = 1)
356#elif defined(BIGENDIAN)
357      PARAMETER (IW1 = 1, IW2 = 2)
358#else
359#define NOPACKING
360#endif
361C
362#if defined (LITTLEENDIAN) || defined(BIGENDIAN)
363      PARAMETER (JMASK = z'3C000000')
364C
365#if defined (VAR_GFORTRAN)
366      ! special code to avoid gfortran out of range error
367      INTEGER*4 Ftmp, F0000000
368      PARAMETER ( Ftmp     = z'10000000' )
369      PARAMETER ( F0000000 = -Ftmp)
370      DATA IMASK / z'00000000', z'000FF000',
371     &             z'00000000', z'00000FF0',
372     &               F0000000 , z'0000000F',
373     &             z'0FF00000', z'00000000',
374     &             z'000FF000', z'00000000',
375     &             z'00000FF0', z'00000000',
376     &             z'0000000F', z'00000000'/
377#else
378      DATA IMASK / z'00000000', z'000FF000',
379     &             z'00000000', z'00000FF0',
380     &             z'F0000000', z'0000000F',
381     &             z'0FF00000', z'00000000',
382     &             z'000FF000', z'00000000',
383     &             z'00000FF0', z'00000000',
384     &             z'0000000F', z'00000000'/
385#endif
386#else
387# define NOPACKING
388#endif
389C
390C---------------------------------------------------------------
391C     check if packing can be used on this archicture:
392C---------------------------------------------------------------
393C
394#if defined NOPACKING
395      ENABLED = .FALSE.
396      RETURN
397#else
398      ENABLED = .TRUE.
399#endif
400C
401C---------------------------------------------------------------
402C     save packing accuracy on common block
403C---------------------------------------------------------------
404C
405      PCKTHRS = THRESHOLD
406C
407C---------------------------------------------------------------
408C     loop over all exponent values and find the maximum errors
409C     introduced by truncating at a certain byte length
410C----------------------------------------------------------------
411C
412      DO IEXP = 0, 255
413C
414         IF (LOCDBG) WRITE (LUPRI,'(//A,I5)') 'IEXP = ', IEXP
415C
416         IPCKTABLE(IEXP) = 8
417C
418         R8WORD0 = ZERO
419C
420         IF (IEXP.NE.0) THEN
421C
422           I4WORD0(IW1) = ISHFT(IEXP,24)
423           IF (LOCDBG) CALL R8BITREP(R8WORD0)
424C
425           I4WORD0(IW1) = ISHFTC(I4WORD0(IW1),-4,30)
426           IF (LOCDBG) CALL R8BITREP(R8WORD0)
427C
428           IF (.NOT. BTEST(I4WORD0(IW1),30) ) THEN
429             I4WORD0(IW1) = I4WORD0(IW1) + JMASK
430           END IF
431C
432         END IF
433C
434         IF (LOCDBG) THEN
435           CALL R8BITREP(R8WORD0)
436           WRITE (LUPRI,'(A,E20.10,/)') 'R8WORD0:',R8WORD0
437         END IF
438C
439         R8WORD1 = R8WORD0
440         IF (IEXP.EQ.0) THEN
441           I4WORD1(IW1) = I4WORD1(IW1) + JMASK
442         END IF
443
444         DO L = 7, 1, -1
445           I4WORD1(IW1) = I4WORD1(IW1) + IMASK(2,L)
446           I4WORD1(IW2) = I4WORD1(IW2) + IMASK(1,L)
447           ERRMAX(L) = DABS(R8WORD0-R8WORD1)
448           IF ( ERRMAX(L) .LT. PCKTHRS ) IPCKTABLE(IEXP) = L
449         END DO
450         ERRMAX(8) = ERRMAX(7) / 256.0D0
451C
452         IF (LOCDBG) THEN
453           WRITE (LUPRI,'(/A)') ' #Bytes / Truncation error:'
454           WRITE (LUPRI,'(8I10)') (L,L=1,8)
455           WRITE (LUPRI,'(8E10.2)') (ERRMAX(L),L=1,8)
456           WRITE (LUPRI,'(A,I5,A,8E10.2)') ' selected length:',
457     &          IPCKTABLE(IEXP),
458     &                    ' -- accepted error:',ERRMAX(IPCKTABLE(IEXP))
459         END IF
460C
461      END DO
462C
463C     RETURN
464      END
465*=====================================================================*
466      SUBROUTINE R8BITREP(R8WORD)
467*---------------------------------------------------------------------*
468C
469C     Purpose: print bit representation of a real*8 number
470C
471C     Christof Haettig, August 1998
472C
473*---------------------------------------------------------------------*
474      IMPLICIT NONE
475C
476#include "priunit.h"
477C
478      REAL*8    R8WORD, R8COPY
479      INTEGER*8 K, J, I8COPY
480      CHARACTER*64 STRING
481
482      EQUIVALENCE (R8COPY,I8COPY)
483
484      R8COPY = R8WORD
485
486      DO J = 1, 64
487         IF ( BTEST(I8COPY,J-1) )  THEN
488            STRING(J:J) = '1'
489         ELSE
490            STRING(J:J) = '0'
491         END IF
492      END DO
493
494      WRITE (LUPRI,'(4(1X,A16))') STRING(01:16), STRING(17:32),
495     &                     STRING(33:48), STRING(49:64)
496
497      RETURN
498      END
499*=====================================================================*
500      SUBROUTINE I4BITREP(I4WORD)
501*---------------------------------------------------------------------*
502C
503C     Purpose: print bit representation of a i*4 number
504C
505C     Christof Haettig, August 1998
506C
507*---------------------------------------------------------------------*
508      IMPLICIT NONE
509C
510#include "priunit.h"
511C
512      INTEGER*4 I4WORD, I4COPY
513      INTEGER J
514      CHARACTER*32 STRING
515
516      I4COPY = I4WORD
517
518      DO J = 1, 32
519        IF ( BTEST(I4COPY,J-1) )  THEN
520           STRING(J:J) = '1'
521        ELSE
522           STRING(J:J) = '0'
523        END IF
524      END DO
525
526      WRITE (LUPRI,'(2(1X,A16))') STRING(1:16), STRING(17:32)
527      WRITE (LUPRI,*) I4COPY
528
529      RETURN
530      END
531*=====================================================================*
532