quad_emu_mod.F90

Go to the documentation of this file.

#if defined(NECSX) && defined(REALHUGE)

MODULE quad_emu_mod

!        NEC 2007. Emulation of vectorized quadruple precision

USE parkind2  ,ONLY : jprh

IMPLICIT NONE

EXTERNAL qvmulqv, ismulqv, qvaddqv, qvsubqv
EXTERNAL r16vcpqv, qvcpr16v

TYPE :: quads
  REAL(8) :: h, l
END TYPE quads
TYPE :: quadv
  REAL(8), POINTER, DIMENSION(:) :: h, l
END TYPE quadv

CONTAINS

SUBROUTINE qxmy (V1,V2,N,V3)
USE parkind1  ,ONLY : jpim
USE parkind2  ,ONLY : jprh
  ! ***************************
  ! V3(1:N) = V1(1:N) * V2(1:N)
  ! ***************************
  INTEGER(KIND=JPIM), INTENT(IN) :: N
  REAL(KIND=JPRH), DIMENSION(N), INTENT(IN)  :: V1, V2
  REAL(KIND=JPRH), DIMENSION(N), INTENT(OUT) :: V3
  type(quadv),SAVE :: qv1, qv2, qv3
  LOGICAL FIRST
  DATA first / .true. /
  INTEGER, SAVE :: NMAX = -1

  IF (first  .OR.  n>nmax) THEN
    IF(.NOT.first)DEALLOCATE (qv1%H,qv1%L,qv2%H,qv2%L,qv3%H,qv3%L)
    first = .false. 
    nmax = n
    ALLOCATE ( qv1%H(nmax) , qv1%L(nmax) )
    ALLOCATE ( qv2%H(nmax) , qv2%L(nmax) )
    ALLOCATE ( qv3%H(nmax) , qv3%L(nmax) )
  ENDIF

  ! QV1(:) = V1(:)
  CALL r16vcpqv(v1,qv1%H,qv1%L,n)

  ! QV2(:) = V2(:)
  CALL r16vcpqv(v2,qv2%H,qv2%L,n)

  ! QV3(:) = QV1(:) * QV2(:)

  CALL qvmulqv(qv1%H(1),qv1%L(1),qv2%H(1),qv2%L(1),qv3%H(1),qv3%L(1),n)

  ! V3(:) = QV3(:)
  CALL qvcpr16v(qv3%H,qv3%L,v3,n)

RETURN
END SUBROUTINE qxmy

SUBROUTINE qaxpy (ALPHA,V1,V2,N)
USE parkind1  ,ONLY : jpim
USE parkind2  ,ONLY : jprh
  ! ***********************************
  ! V2(1:N) = ALPHA * V1(1:N) + V2(1:N)
  ! ***********************************
  INTEGER(KIND=JPIM), INTENT(IN) :: N
  REAL(KIND=JPRH),               INTENT(IN)    :: ALPHA
  REAL(KIND=JPRH), DIMENSION(N), INTENT(IN)    :: V1
  REAL(KIND=JPRH), DIMENSION(N), INTENT(INOUT) :: V2
  type(quadv),SAVE :: qalpha
  type(quadv),SAVE :: qv1, qv2, qvtmp
  LOGICAL FIRST
  DATA first / .true. /
  INTEGER, SAVE :: NMAX = -1

  IF (first  .OR.  n>nmax) THEN
    IF(.NOT.first)THEN
      DEALLOCATE (qalpha%H,qalpha%L,qvtmp%H,qvtmp%L,qv1%H,qv1%L,qv2%H,qv2%L)
    ENDIF
    first = .false.
    nmax = n
    ALLOCATE ( qalpha%H(nmax) , qalpha%L(nmax) )
    ALLOCATE ( qv1%H(nmax) , qv1%L(nmax) )
    ALLOCATE ( qv2%H(nmax) , qv2%L(nmax) )
    ALLOCATE ( qvtmp%H(nmax) , qvtmp%L(nmax) )
  ENDIF

  ! QALPHA(:) = ALPHA
  CALL r16scpqv(alpha,qalpha%H,qalpha%L,n)

  ! QV1(:) = V1(:)
  CALL r16vcpqv(v1,qv1%H,qv1%L,n)

  ! QV2(:) = V2(:)
  CALL r16vcpqv(v2,qv2%H,qv2%L,n)

  ! QVtmp(:) = ALPHA*QV1(:)
  CALL qvmulqv(qalpha%H,qalpha%L,qv1%H,qv1%L,qvtmp%H,qvtmp%L,n)

  ! QVtmp(:) = QVtmp(:) + QV2(:)
  CALL qvaddqv(qvtmp%H,qvtmp%L,qv2%H,qv2%L,qvtmp%H,qvtmp%L,n)

  ! V2(:) = QVtmp(:)
  CALL qvcpr16v(qvtmp%H,qvtmp%L,v2,n)

RETURN
END SUBROUTINE qaxpy

SUBROUTINE sqrtivdv (I1, I2, N, X)
USE parkind1  ,ONLY : jpim, jpib
USE parkind2  ,ONLY : jprh
  ! *********************************************
  ! X(1:N) = SQRT ( REAL(I1(1:N) / REAL(I2(1:N) )
  ! *********************************************
  INTEGER(KIND=JPIM), INTENT(IN) :: N
  INTEGER(KIND=JPIB), INTENT(IN), DIMENSION(N) :: I1
  INTEGER(KIND=JPIB), INTENT(IN), DIMENSION(N) :: I2
  REAL(KIND=JPRH), INTENT(INOUT), DIMENSION(N) :: X
  type(quadv),SAVE :: qv1, qv2, qv3, qvx
  LOGICAL FIRST
  DATA first / .true. /
  INTEGER, SAVE :: NMAX = -1

  IF (first  .OR.  n>nmax) THEN
    IF(.NOT.first)THEN
      DEALLOCATE (qv1%H,qv1%L,qv2%H,qv2%L,qv3%H,qv3%L)
      DEALLOCATE (qvx%H,qvx%L)
    ENDIF
    first = .false.
    nmax = n
    ALLOCATE ( qv1%H(nmax) , qv1%L(nmax) )
    ALLOCATE ( qv2%H(nmax) , qv2%L(nmax) )
    ALLOCATE ( qv3%H(nmax) , qv3%L(nmax) )
    ALLOCATE ( qvx%H(nmax) , qvx%L(nmax) )
  ENDIF
  ! QV1(:) = I1(:)
  CALL i4vcpqv(i1,qv1%H,qv1%L,n)
  ! QV2(:) = I2(:)
  CALL i4vcpqv(i2,qv2%H,qv2%L,n)
  ! QV3(:) = I1(:) / I2(:)
  CALL qvdivqv(qv1%H,qv1%L,qv2%H,qv2%L,qv3%H,qv3%L,n)
  ! QVx(:) = sqrt ( QV3(:) )
  CALL qvsqrt (qv3%H,qv3%L,qvx%H,qvx%L,n)
  ! X(:) = QVx(:)
  CALL qvcpr16v(qvx%H,qvx%L,x,n)

RETURN
END SUBROUTINE sqrtivdv

END MODULE quad_emu_mod



! New routines ------------------------------------------------
      SUBROUTINE i4vcpqv(I4,AHI,ALO,N)
! For SX vectorized version
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER(KIND=JPIM) :: N
      INTEGER(KIND=JPIB) :: I4(n)
      REAL(KIND=JPRB) :: AHI(n),ALO(n)
      REAL(KIND=JPRH) :: R16(n)

      r16(1:n) = i4(1:n)
      CALL r16vcpqv(r16,ahi,alo,n)

      END SUBROUTINE i4vcpqv

      SUBROUTINE r16vcpqv(R16,AHI,ALO,N)
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
#ifdef SCALAR
! Fortran equivalent code
       INTEGER(KIND=JPIM) :: N
       REAL(KIND=JPRH) :: R16(n)
       REAL(KIND=JPRB) :: AHI(n),ALO(n)
       INTEGER I
       DO i=1,n
         ahi(i)=r16(i)
         alo(i)=r16(i)-ahi(i)
       ENDDO
#else
! For SX vectorized version
      INTEGER(KIND=JPIM) :: N
      REAL(KIND=JPRB) :: R16(2,n)
      REAL(KIND=JPRB) :: AHI(n),ALO(n)

      REAL(KIND=JPRB) :: AA,BB
      INTEGER(KIND=JPIB) :: II,JJ,MSK1,MSK2,MSK3
      equivalence(aa,ii)
      equivalence(bb,jj)
      parameter(msk1=z'7FF0000000000000')
      parameter(msk2=z'FFF0000000000000')
      parameter(msk3=z'0340000000000000')

      INTEGER I

      DO i=1,n
        ahi(i)=r16(1,i)
        bb=r16(1,i)
        jj=iand(jj,msk1)
        IF(jj.LT.msk3) THEN
          alo(i)=0.d0
        ELSE
          aa=r16(2,i)
          ii=iand(ii,msk2)
          alo(i)=r16(2,i)-aa
        ENDIF
      ENDDO
#endif
      END SUBROUTINE r16vcpqv

      SUBROUTINE qvcpr16v(AHI,ALO,R16,N)
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
#ifdef SCALAR
! Fortran equivalent code
       INTEGER N
       REAL(KIND=JPRB) :: AHI(n),ALO(n)
       REAL(KIND=JPRH) :: R16(n)
       INTEGER I
       DO i=1,n
         r16(i)=ahi(i)
         r16(i)=r16(i)+alo(i)
       ENDDO
#else
! For SX vectorized version
      REAL(KIND=JPRB) :: AHI(n),ALO(n)
      REAL(KIND=JPRB) :: R16(2,n)
      INTEGER N
!
      REAL(KIND=JPRB) :: A1,A2,AA,ZHI,ZLO,AGOMI,R
      INTEGER(KIND=JPIB) :: II,MSK
      equivalence(aa,ii)
      parameter(msk=z'FFF0000000000000')
      INTEGER I
!
      DO i=1,n
         a1=ahi(i)
         a2=alo(i)
         aa=a1
         ii=iand(ii,msk)
         agomi=aa*2d0**(-52)
         zhi=a1
         zlo=a2+agomi
         r=a1*a2
         IF (r.LT.0) zhi=a1-agomi
         r16(1,i)=zhi
         r16(2,i)=(a1-zhi)+zlo
      ENDDO
#endif
      END SUBROUTINE qvcpr16v
! New routines ------------------------------------------------

      SUBROUTINE r16scpqv(R16,AHI,ALO,N)
! ... Quadruple precision native to Quadruple emulation pair
!     QV(1:N) = R16
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: R16(2)
      REAL(KIND=JPRB) :: AHI(n),ALO(n)

      REAL(KIND=JPRB) :: AA
      INTEGER(KIND=JPIB) :: II,MSK
      equivalence(aa,ii)
      parameter(msk=z'FFF0000000000000')

      INTEGER I

      DO i=1,n
        ahi(i)=r16(1)
        aa=r16(2)
        ii=iand(ii,msk)
        alo(i)=r16(2)-aa
      ENDDO

      END SUBROUTINE r16scpqv

      SUBROUTINE qvmulqv(AHI,ALO,CHI,CLO,ACHI,ACLO,N)
! ... multiplication of Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),CHI(n),CLO(n),ACHI(n),ACLO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: ZHI,ZLO,ZZ
      REAL(KIND=JPRB) :: A1,A2,C1,C2,T
      INTEGER I

!     longmul(a, c) RESULT(ac)
!     z = exactmul2(a%hi, c%hi)
!     zz = ((a%hi + a%lo) * c%lo + a%lo * c%hi) + z%lo
!     ac%hi = z%hi + zz
!     ac%lo = (z%hi - ac%hi) + zz

!     exactmul2(a, c) RESULT(ac)
!     t = constant * a
!     a1 = (a - t) + t
!     a2 = a - a1
!     t = constant * c
!     c1 = (c - t) + t
!     c2 = c - c1
!     ac%hi = a * c
!     ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2

      DO i=1,n
        t = constant * ahi(i)
        a1 = (ahi(i) - t) + t
        a2 = ahi(i) - a1
        t = constant * chi(i)
        c1 = (chi(i) - t) + t
        c2 = chi(i) - c1
        zhi = ahi(i) * chi(i)
        zlo = (((a1 * c1 - zhi) + a1 * c2) + c1 * a2) + c2 * a2
        zz = ((ahi(i) + alo(i)) * clo(i) + alo(i) * chi(i)) + zlo
        achi(i) = zhi + zz
        aclo(i) = (zhi - achi(i)) + zz
      ENDDO

      END SUBROUTINE qvmulqv
      SUBROUTINE qvdivqv(AHI,ALO,CHI,CLO,ACHI,ACLO,N)
! ... division of Quadruple emulation pair 
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),CHI(n),CLO(n),ACHI(n),ACLO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: Z,ZZ,QHI,QLO
      REAL(KIND=JPRB) :: A1,A2,C1,C2,T
      INTEGER I

!     longdiv(a, c) RESULT(ac)
!     z = a%hi / c%hi
!     q = exactmul2(c%hi, z)
!     zz = ((((a%hi - q%hi) - q%lo) + a%lo) - z*c%lo) / (c%hi + c%lo)
!     ac%hi = z + zz
!     ac%lo = (z - ac%hi) + zz

!     exactmul2(a, c) RESULT(ac)
!     t = constant * a
!     a1 = (a - t) + t
!     a2 = a - a1
!     t = constant * c
!     c1 = (c - t) + t
!     c2 = c - c1
!     ac%hi = a * c
!     ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2

      DO i=1,n
        z = ahi(i)/chi(i)
        t = constant * chi(i)
        a1 = (chi(i) - t) + t
        a2 = chi(i) - a1
        t = constant * z
        c1 = (z - t) + t
        c2 = z - c1
        qhi = chi(i) * z
        qlo = (((a1 * c1 - qhi) + a1 * c2) + c1 * a2) + c2 * a2
        zz = ((((ahi(i) - qhi) - qlo) + alo(i)) - z*clo(i)) / (chi(i) + clo(i))
        achi(i) = z + zz
        aclo(i) = (z - achi(i)) + zz
      ENDDO

      END SUBROUTINE qvdivqv
      SUBROUTINE qvaddqv(AHI,ALO,CHI,CLO,ACHI,ACLO,N)
! ... addition of Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),CHI(n),CLO(n),ACHI(n),ACLO(n)

      REAL(KIND=JPRB) :: Z,ZZ,Q
      INTEGER I

!     longadd(a, c) RESULT(ac)
!     z = a%hi + c%hi
!     q = a%hi - z
!     zz = (((q + c%hi) + (a%hi - (q + z))) + a%lo) + c%lo
!     ac%hi = z + zz
!     ac%lo = (z - ac%hi) + zz

      DO i=1,n
        z = ahi(i) + chi(i)
        q = ahi(i) - z
        zz = (((q + chi(i)) + (ahi(i) - (q + z))) + alo(i)) + clo(i)
        achi(i) = z + zz
        aclo(i) = (z - achi(i)) + zz
      ENDDO

      END SUBROUTINE qvaddqv
      SUBROUTINE qvsubqv(AHI,ALO,CHI,CLO,ACHI,ACLO,N)
! ... subtraction of Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),CHI(n),CLO(n),ACHI(n),ACLO(n)

      REAL(KIND=JPRB) :: Z,ZZ,Q
      INTEGER I

!     longsub(a, c) RESULT(ac)
!     z = a%hi - c%hi
!     q = a%hi - z
!     zz = (((q - c%hi) + (a%hi - (q + z))) + a%lo) - c%lo
!     ac%hi = z + zz
!     ac%lo = (z - ac%hi) + zz

      DO i=1,n
        z = ahi(i) - chi(i)
        q = ahi(i) - z
        zz = (((q - chi(i)) + (ahi(i) - (q + z))) + alo(i)) - clo(i)
        achi(i) = z + zz
        aclo(i) = (z - achi(i)) + zz
      ENDDO

      END SUBROUTINE qvsubqv
      SUBROUTINE dvmuldv(A,C,ACHI,ACLO,N)
! ... multiplication of double precision to Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: A(n),C(n),ACHI(n),ACLO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: A1,A2,C1,C2,T
      INTEGER I

!     exactmul2(a, c) RESULT(ac)
!     t = constant * a
!     a1 = (a - t) + t
!     a2 = a - a1
!     t = constant * c
!     c1 = (c - t) + t
!     c2 = c - c1
!     ac%hi = a * c
!     ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2

      DO i=1,n
        t = constant * a(i)
        a1 = (a(i) - t) + t
        a2 = a(i) - a1
        t = constant * c(i)
        c1 = (c(i) - t) + t
        c2 = c(i) - c1
        achi(i) = a(i) * c(i)
        aclo(i) = (((a1 * c1 - achi(i)) + a1 * c2) + c1 * a2) + c2 * a2
      ENDDO

      END SUBROUTINE dvmuldv
      SUBROUTINE ismulqv(IS,QHI,QLO,QIHI,QILO,N)
! ... integer scale of Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER IS,N
      REAL(KIND=JPRB) :: QHI(n),QLO(n),QIHI(n),QILO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: A1,A2,Q1,Q2,T,ZHI,ZLO,ZZ
      INTEGER I

      t = constant * is
      a1 = (dble(is) - t) + t
      a2 = dble(is) - a1
      DO i=1,n
        t = constant * qhi(i)
        q1 = (qhi(i) - t) + t
        q2 = qhi(i) - q1
        zhi = is * qhi(i)
        zlo = (((a1 * q1 - zhi) + a1 * q2) + q1 * a2) + q2 * a2
        zz = is * qlo(i) + zlo
        qihi(i) = zhi + zz
        qilo(i) = (zhi - qihi(i)) + zz
      ENDDO

      END SUBROUTINE ismulqv
      SUBROUTINE qvdivis(AHI,ALO,IS,AIHI,AILO,N)
! ... integer division of Quadruple emulation pair
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER IS,N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),AIHI(n),AILO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: Z,ZZ,QHI,QLO
      REAL(KIND=JPRB) :: A1,A2,C1,C2,T
      INTEGER I

!     div_quad_int(a, b) RESULT(c)
!     Divide a quadruple-precision number (a) by an integer (b)
!     z = a%hi / b
!     q = exactmul2(DBLE(b), z)
!     zz = (((a%hi - q%hi) - q%lo) + a%lo) / b
!     c%hi = z + zz
!     c%lo = (z - c%hi) + zz

!     exactmul2(a, c) RESULT(ac)
!     t = constant * a
!     a1 = (a - t) + t
!     a2 = a - a1
!     t = constant * c
!     c1 = (c - t) + t
!     c2 = c - c1
!     ac%hi = a * c
!     ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2

      t = constant * is
      a1 = (dble(is) - t) + t
      a2 = dble(is) - a1
      DO i=1,n
        z = ahi(i)/is
        t = constant * z
        c1 = (z - t) + t
        c2 = z - c1
        qhi = dble(is) * z
        qlo = (((a1 * c1 - qhi) + a1 * c2) + c1 * a2) + c2 * a2
        zz = (((ahi(i) - qhi) - qlo) + alo(i)) / is
        aihi(i) = z + zz
        ailo(i) = (z - aihi(i)) + zz
      ENDDO

      END SUBROUTINE qvdivis

      SUBROUTINE qvsqrt(AHI,ALO,BHI,BLO,N)
      USE parkind1  ,ONLY : jpim, jpib, jprb
      USE parkind2  ,ONLY : jprh
      INTEGER N
      REAL(KIND=JPRB) :: AHI(n),ALO(n),BHI(n),BLO(n)

      INTEGER, PARAMETER :: NBITS = digits(1.0d0)
      REAL(KIND=JPRB), PARAMETER :: CONSTANT = 2.d0**(nbits-nbits/2)+1.d0

      REAL(KIND=JPRB) :: A1,A2,C1,C2,T
      REAL(KIND=JPRB) :: TTHI,TTLO,TT,RES
      INTEGER I

!     longsqrt(a) RESULT(b)
!     t = SQRT(a%hi)
!     tt = exactmul2(t, t)
!     res = (((a%hi - tt%hi) - tt%lo) + a%lo) * 0.5_dp / t
!     b%hi = t + res
!     b%lo = (t - b%hi) + res

!     exactmul2(a, c) RESULT(ac)
!     t = constant * a
!     a1 = (a - t) + t
!     a2 = a - a1
!     t = constant * c
!     c1 = (c - t) + t
!     c2 = c - c1
!     ac%hi = a * c
!     ac%lo = (((a1 * c1 - ac%hi) + a1 * c2) + c1 * a2) + c2 * a2

      DO i=1,n
        tt = sqrt(ahi(i))
        t = constant * tt
        a1 = (tt - t) + t
        a2 = tt - a1
        t = constant * tt
        c1 = (tt - t) + t
        c2 = tt - c1
        tthi = tt * tt
        ttlo = (((a1 * c1 - tthi) + a1 * c2) + c1 * a2) + c2 * a2
        res = (((ahi(i) - tthi) - ttlo) + alo(i)) * 0.5d0 / tt
        bhi(i) = tt + res
        blo(i) = (tt - bhi(i)) + res
      ENDDO

      END SUBROUTINE qvsqrt
#else
MODULE quad_emu_mod
END MODULE quad_emu_mod
#endif