w3ai00.f

Go to the documentation of this file.

C> @file
C> @brief Real array to 16 bit packed format.
C> @author Ralph Jones @date 1985-07-31
 
C> Converts IEEE floating point numbers to 16 bit
C> packed office note 84 format. The floating point number are
C> converted to 16 bit signed scaled integers.
C>
C> Program history log:
C> - Ralph Jones 1989-10-20 Convert cyber 205 version of w3ai00 to cray.
C> - Ralph Jones 1990-03-18 Change to use cray integer*2 packer.
C> - Ralph Jones 1990-10-11 Special version to pack grids larger than
C> 32743 words. Will do old and new version.
C> - Ralph Jones 1991-02-16 Changes so equivalence of pack and real8
C> arrays will work.
C> - Ralph Jones 1993-06-10 Changes for array size (512,512) 262144 words.
C> - Boi Vuong 1998-03-10 Remove the cdir$ integer=64 directive.
C> - Stephen Gilbert 1998-11-18 Changed to pack IEEE values for the IBM SP
C>
C> @param[in] REAL8 Array of cray floating point numbers.
C> @param[in] LABEL Six 8-byte integer words. Must have first 8 of 12 32 bit
C> word office note 84 label. word 6 must have in bits 31-00 the number of
C> real words in array real8 if j is greater than 32743.
C> j in bits 15-0 of the 4th id word is set zero.
C> @param[out] PACK Packed output array of integer words of size 6 + (j+3)/4 ,
C> j = no. points in label (from word 4 bits 15-00). Label will be copied to pack words 1-4.
C> - Pack will contain the following in words 5-6:
C>  - word 5 bits 63-48 Number of bytes in whole record. will not be correct if j > 32743.
C>  - word 5 bits 47-32 Exclusive-or checksum by 16 bit words of whole array pack excluding checksum itself.
C>  - word 5 bits 31-00 Center value a = mean of max and min values. converted to ibm 32 floating point number.
C>  - word 6 bits 63-48 Zero.
C>  - word 6 bits 47-32 16 bit shift value n. the least integer such that abs(x-a)/2**n lt 1 for all x in real8. limited to +-127.
C>  - word 6 bits 31-00 Number of words in real8 if > 32743, right adjusted if <= 32743 set zero.
C>
C> @note Pack and label may be equivalenced.  n, the number of points
C> in a grid is now in 32 bit id word 12.
C>
C> @author Ralph Jones @date 1985-07-31
      SUBROUTINE w3ai00(REAL8,PACK,LABEL)
C
       REAL           REAL8(*)
       REAL           XX(262144)
C
       INTEGER(8)     KK(262144)
       INTEGER(8)     LABEL(6)
       INTEGER(8)     PACK(*)
       INTEGER(8)     TPACK(6)
       INTEGER(8)     MASK16,MASK32,MASKN,IBYTES,IXOR
       INTEGER(8)     IB,N
       REAL(8)        B
       REAL(4)        X,A
       real(4)        rtemp(2)
       integer(8)     irtemp
       equivalence(irtemp,rtemp(1))
C
       SAVE
C
       equivalence(b,ib)
C
       DATA  mask16/z'000000000000FFFF'/
       DATA  mask32/z'00000000FFFFFFFF'/
       DATA  maskn /z'0000FFFF00000000'/
C
C TRANSFER LABEL DATA TO WORDS 1-4.  GET WORD COUNT, COMPUTE BYTES.
C
       DO 10 i = 1,4
         tpack(i) = label(i)
 10    CONTINUE
C
       tpack(5) = 0
       tpack(6) = 0
C
C      GET J, THE NUMBER OF WORDS IN A GRID, IF ZERO GET THE
C      GET J FROM OFFICE NOTE 84 ID WORD 12.
C
       j       = iand(mask16,tpack(4))
       IF (j.EQ.0) THEN
         tpack(6) = label(6)
         j       = iand(mask32,tpack(6))
         IF (j.EQ.0) THEN
           print *,' W3AI00: ERROR, NO. OF WORDS IN GRID = 0'
           RETURN
         ENDIF
         IF (j.GT.262144) THEN
           print *,' W3AI00: ERROR, NO. OF WORDS IN GRID = ',j
           print *,' THERE IS A LIMIT OF 262144 WORDS.'
           RETURN
         ENDIF
       ENDIF
       m       = j + 24
C
C      COMPUTE THE NUMBER OF 64 BIT INTEGER CRAY WORDS NEEDED FOR
C      PACKED DATA.
C
       IF (mod(m,4).NE.0) THEN
         iword = (m + 3) / 4
       ELSE
         iword = m / 4
       ENDIF
C
       ibytes = m + m
C
C      STORE NUMBER OF BYTES IN RECORD IN BITS 63-48 OF WORD 5.
C      BITS ARE NUMBERED LEFT TO RIGHT 63 T0 00
C
       tpack(5) = ishft(ibytes,48_8)
C
C FIND MAX, MIN OF DATA, COMPUTE A AND N.
C
         rmax = real8(1)
         rmin = rmax
         DO 20 i = 2,j
           rmax = amax1(rmax,real8(i))
           rmin = amin1(rmin,real8(i))
 20      CONTINUE
C
         a = 0.5 * (rmax + rmin)
         x = rmax - a
         IF (rmax.NE.rmin) THEN
C           CALL USDCTI(X,B,1,1,ISTAT)
           CALL q9e3i6(x,b,1,istat)
           IF (istat.NE.0) print *,' W3AI00-USDCTI OVERFLOW ERROR 1'
           n = iand(ishft(ib,-56_8),127_8)
           n = 4 * (n - 64)
           IF (btest(ib,55_8)) GO TO  30
           n = n - 1
           IF (btest(ib,54_8)) GO TO  30
           n = n - 1
           IF (btest(ib,53_8)) GO TO  30
           n = n - 1
 30      CONTINUE
           n = max(-127_8,min(127_8,n))
         ELSE
C
C        FIELD IS ZERO OR A CONSTANT
C
           n = 0
         ENDIF
C
C     CONVERT AVERAGE VALUE FROM IEEE F.P. TO IBM370 32 BIT
C     STORE IBM370 32 BIT F.P. AVG. VALUE IN BITS 31 - 00 OF WORD 5.
C
C         CALL USSCTI(A,TPACK(5),5,1,ISTAT)
         CALL q9ei32(a,rtemp(2),1,istat)
           IF (istat.NE.0) print *,' W3AI00-USDCTI OVERFLOW ERROR 2'
         tpack(5)=ior(tpack(5),irtemp)
C
C      STORE SCALING VALUE N IN BITS 47 - 32 OF WORD 6.
C
         tpack(6)  = ior(iand(maskn,ishft(n,32_8)),tpack(6))
C
C NOW PACK UP THE DATA, AND SCALE IT TO FIT AN INTEGER*2 WORD
C
         twon = 2.0 ** (15 - n)
         DO 40 i = 1,j
           xx(i) = (real8(i) - a) * twon
           kk(i) = xx(i) + sign(0.5,xx(i))
           IF (kk(i).GE.(-32767)) THEN
             kk(i) = min(32767_8,kk(i))
           ELSE
             kk(i) = -32767
           ENDIF
             kk(i) = iand(kk(i),mask16)
 40      CONTINUE
C
C        SHIFT THE INTEGER*2 DATA TO FIT 4 IN A 64 BIT WORD
C
         lim  = (j / 4 ) * 4
         irem = j - lim
         DO 50 i = 1,lim,4
           kk(i)   = ishft(kk(i),  48_8)
           kk(i+1) = ishft(kk(i+1),32_8)
           kk(i+2) = ishft(kk(i+2),16_8)
 50      CONTINUE
C
C        SHIFT THE REMAINING 1, 2, OR 3 INTEGER*2 WORDS
C
         IF (irem.EQ.1) THEN
           kk(lim+1)   = ishft(kk(lim+1),48_8)
         ENDIF
C
         IF (irem.EQ.2) THEN
           kk(lim+1)   = ishft(kk(lim+1),48_8)
           kk(lim+2)   = ishft(kk(lim+2),32_8)
         ENDIF
C
         IF (irem.EQ.3) THEN
           kk(lim+1)   = ishft(kk(lim+1),48_8)
           kk(lim+2)   = ishft(kk(lim+2),32_8)
           kk(lim+3)   = ishft(kk(lim+3),16_8)
         ENDIF
C
C        PACK THE DATA BY USE OF IOR FOUR TO A WORD
C
         ii = 7
         DO 60 i = 1,lim,4
           pack(ii) = ior(ior(ior(kk(i),kk(i+1)),kk(i+2)),kk(i+3))
           ii       = ii + 1
 60      CONTINUE
C
C        PACK THE LAST 1, 2, OR 3 INTEGER*2 WORDS
C
         IF (irem.EQ.1) THEN
           pack(iword) = kk(lim+1)
         ENDIF
C
         IF (irem.EQ.2) THEN
           pack(iword) = ior(kk(i),kk(i+1))
         ENDIF
C
         IF (irem.EQ.3) THEN
           pack(iword) = ior(ior(kk(i),kk(i+1)),kk(i+2))
         ENDIF
C
C      MOVE LABEL FROM TEMPORARY ARRAY TO PACK
C
       DO 70 i = 1,6
         pack(i) = tpack(i)
 70    CONTINUE
C
C      COMPUTE CHECKSUM AND STORE
C
         ixor = 0
C
C      COMPUTES A 64 BIT CHECKSUM 1ST
C
         DO 80 i = 1,iword
           ixor = ieor(ixor,pack(i))
 80      CONTINUE
C
C      COMPUTES A 32 BIT CHECKSUM 2ND
C
         ixor = ieor(ishft(ixor,-32_8),iand(ixor,mask32))
C
C      COMPUTES A 16 BIT CHECKSUM 3RD
C
         ixor = ieor(ishft(ixor,-16_8),iand(ixor,mask16))
C
C      STORE 16 BIT CHECK SUM OF RECORD IN BITS 47-32 OF WORD 5.
C
       pack(5) = ior(ishft(ixor,32_8),pack(5))
C
       RETURN
       END
 
 
C> Convert IEEE 32 bit task 754 floating point numbers
C> to IBM370 32 bit floating point numbers.
C>
C> Program history log:
C> - Ralph Jones 1990-06-04 Convert to sun fortran 1.3.
C> - Ralph Jones 1990-07-14 Change ishft to lshift or lrshft.
C> - Ralph Jones 1991-03-28 Change to silicongraphics 3.3 fortran 77.
C> - Ralph Jones 1992-07-20 Change to ibm aix xl fortran.
C> - Ralph Jones 1995-11-15 Add save statement.
C> - Stepen Gilbert 1998-11-18 Specified 4-byte Integer values.
C>
C> @param[in] A - Real*4 array of IEEE 32 bit floating point numbers.
C> @param[in] N - Number of words to convert to IBM370 32 bit F.P.
C> @param[out] B - Real*4 array of IBM370 32 bit floating point numbers.
C> @param[out] ISTAT:
C> - 0: All numbers converted.
C> - -1: N is less than one.
C> - +K: K infinity or nan numbers were found.
C>
C> @note See IEEE task 754 standard floating point arithmetic for
C> more information about IEEE F.P.
C>
C> @author Ralph Jones @date 1990-06-04
      SUBROUTINE q9ei32(A,B,N,ISTAT)
C
       INTEGER(4)      A(*)
       INTEGER(4)      B(*)
       INTEGER(4)      SIGN,MASKFR,IBIT8,MASKSN,ITEMP,IBMEXP,IBX7
       INTEGER(4)      ISIGN
C
       SAVE
C
       DATA  maskfr/z'00FFFFFF'/
       DATA  ibit8 /z'00800000'/
       DATA  masksn/z'7FFFFFFF'/
       DATA  sign  /z'80000000'/
C
           IF (n.LT.1) THEN
             istat = -1
             RETURN
           ENDIF
C
           istat = 0
C
         DO 30 i = 1,n
C
C          SIGN BIT OFF
C
           isign = 0
           itemp = a(i)
C
C          TEST SIGN BIT
C
           IF (itemp.EQ.0) GO TO 20
C
           IF (itemp.LT.0)  THEN
C
C            SIGN BIT ON
C
             isign = sign
C
C            TURN SIGN BIT OFF
C
             itemp = iand(itemp,masksn)
C
           END IF
C
           ibmexp = ishft(itemp,-23_4)
C
C          TEST FOR INDIFINITE OR NAN NUMBER
C
           IF (ibmexp.EQ.255) GO TO 10
C
C          TEST FOR ZERO EXPONENT AND FRACTION (UNDERFLOW)
C
           IF (ibmexp.EQ.0) GO TO 20
           ibmexp = ibmexp + 133
           ibx7   = iand(3_4,ibmexp)
           ibmexp = ieor(ibmexp,ibx7)
           ibx7   = ieor(3_4,ibx7)
           itemp  = ior(itemp,ibit8)
           itemp  = ior(ishft(ibmexp,22_4),ishft(iand(itemp,maskfr),
     &     -ibx7))
           b(i)   = ior(itemp,isign)
           GO TO  30
C
 10      CONTINUE
C
C          ADD 1 TO ISTAT FOR INDEFINITE OR NAN NUMBER
C
           istat = istat + 1
C
 20      CONTINUE
           b(i) = 0
C
 30      CONTINUE
C
         RETURN
       END
 
C> Convert ieee 32 bit task 754 floating point numbers
C> to ibm370 64 bit floating point numbers.
C>
C> Program history log:
C> - Ralph Jones 1992-08-02
C> - Ralph Jones 1995-11-15 Add save statement.
C>
C> @param[in] A Real*4 array of IEEE 32 bit floating point numbers.
C> @param[in] N Number of words to convert to IBM370 64 bit F.P.
C> @param[out] B Real*8 array of IBM370 64 bit floating point numbers.
C> @param[out] ISTAT
C> - 0 All numbers converted.
C> - -1 N is less than one.
C> - +K K infinity or nan numbers were found.
C>
C> @note See IEEE task 754 standard floating point arithmetic for
C> more information about IEEE F.P.
C>
C> @author Ralph Jones @date 1992-08-02
      SUBROUTINE q9e3i6(A,B,N,ISTAT)
 
C
       INTEGER(4)      A(N)
       INTEGER(4)      B(2,N)
       INTEGER(4)      SIGN,MASKFR,IBIT8,MASKSN,ITEMP,IEEEXP
       INTEGER(4)      IBMEXP,IBX7,JTEMP,ISIGN
C
       SAVE
C
       DATA  maskfr/z'00FFFFFF'/
       DATA  ibit8 /z'00800000'/
       DATA  masksn/z'7FFFFFFF'/
       DATA  sign  /z'80000000'/
C
           IF (n.LT.1) THEN
             istat = -1
             RETURN
           ENDIF
C
           istat = 0
C
         DO 30 i = 1,n
           isign = 0
           itemp = a(i)
C
C          TEST SIGN BIT
C
           IF (itemp.EQ.0) GO TO 20
C
C          TEST FOR NEGATIVE NUMBERS
C
           IF (itemp.LT.0) THEN
C
C            SIGN BIT ON
C
             isign = sign
C
C            TURN SIGN BIT OFF
C
             itemp = iand(itemp,masksn)
C
           END IF
C
C          GET IEEE EXPONENT
C
           ieeexp = ishft(itemp,-23_4)
C
C          TEST FOR INDIFINITE OR NAN NUMBER
C
           IF (ieeexp.EQ.255) GO TO 10
C
C          TEST FOR ZERO EXPONENT AND FRACTION (UNDERFLOW)
C          CONVERT IEEE EXPONENT (BASE 2) TO IBM EXPONENT
C          (BASE 16)
C
           IF (ieeexp.EQ.0)   GO TO 20
           ibmexp = ieeexp + 133
           ibx7   = iand(3_4,ibmexp)
           ibmexp = ieor(ibmexp,ibx7)
           ibx7   = ieor(3_4,ibx7)
           itemp  = ior(itemp,ibit8)
           jtemp  = ior(ishft(ibmexp,22_4),ishft(iand(itemp,maskfr),
     &     -ibx7))
           b(1,i) = ior(jtemp,isign)
           b(2,i) = 0
           IF (ibx7.GT.0) b(2,i) = ishft(itemp,32_4-ibx7)
           GO TO  30
C
 10      CONTINUE
C          ADD 1 TO ISTAT FOR INDEFINITE OR NAN NUMBER
C
           istat = istat + 1
C
 20      CONTINUE
           b(1,i) = 0
           b(2,i) = 0
C
 30      CONTINUE
C
         RETURN
       END