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