anchor.diff:  
     FUNCTION Quad2Dense(qA) RESULT(A)
       REAL(DOUBLE),DIMENSION(SpAMM_tiles*BLOCK_SIZE,SpAMM_tiles*BLOCK_SIZE) :: A
  -    TYPE(QuTree) :: qA
  +    TYPE(QuTree),POINTER :: qA
       CALL SpAMM_Quad2Mat(A,qA)
     END FUNCTION Quad2Dense
  +
     RECURSIVE SUBROUTINE SpAMM_Quad2Mat(A,qA)
       INTEGER :: I1,I2,J1,J2
       REAL(DOUBLE),DIMENSION(:,:) :: A
  -    TYPE(QuTree) :: qA
  -
  +    TYPE(QuTree),POINTER :: qA
  +
  +    IF(.NOT.ASSOCIATED(qA))RETURN
       IF(qA%Norm==Zero)RETURN
  -
  +    !
       IF(ALLOCATED(qA%Blok))THEN
          I1=qA%Box(1,1)
          I2=qA%Box(1,2)
  
positive.code:        qC%Blok=qC%Blok+CTmp
         ! Bounds
         qC%Norm=qC%Norm+SUM(CTmp(1:BLOCK_SIZE,1:BLOCK_SIZE)**2)
         ! Print
         ! #ifdef SpAMM_PRINT
         ! C_IJ=SUM_K A_IK*B_KJ
         SpAMM_stream%Box(1,:)=qA%Box(2,:) ! I
         SpAMM_stream%Box(2,:)=qA%Box(1,:) ! K
         SpAMM_stream%Box(3,:)=qB%Box(1,:) ! J
         ALLOCATE(SpAMM_stream%Next)
         SpAMM_stream=>SpAMM_stream%Next
         NULLIFY(SpAMM_stream%Next)
  
  !!$
  !       WRITE(*,111)Spamm_Stream%next%Box(1,:),spamm_Stream%next%Box(2,:),spamm_Stream%next%Box(3,:)
  111    FORMAT("Cuboid[{",I6,",",I6,",",I6,"},{",I6,",",I6,",",I6,"}],  ")
  
  
         ! #endif
         RETURN
      ENDIF
      ! New nodes
      IF(.NOT.ASSOCIATED(qC%Quad00)) &
         CALL NewQuNode(qC%Quad00,SpAMM_quadcount)
      IF(.NOT.ASSOCIATED(qC%Quad01)) &
         CALL NewQuNode(qC%Quad01,SpAMM_quadcount)
      IF(.NOT.ASSOCIATED(qC%Quad10)) &
         CALL NewQuNode(qC%Quad10,SpAMM_quadcount)
      IF(.NOT.ASSOCIATED(qC%Quad11)) &
         CALL NewQuNode(qC%Quad11,SpAMM_quadcount)
      ! SpAMM=Sparse Approximate Matrix-Matrix
  
    !  IF(qA%Siz<=2*BLOCK_SIZE)THEN
  
      do00x00=qA%Quad00%Norm*qB%Quad00%Norm>SpAMM_tolerance
      do01x10=qA%Quad01%Norm*qB%Quad10%Norm>SpAMM_tolerance
      do00x01=qA%Quad00%Norm*qB%Quad01%Norm>SpAMM_tolerance
      do01x11=qA%Quad01%Norm*qB%Quad11%Norm>SpAMM_tolerance
      do10x00=qA%Quad10%Norm*qB%Quad00%Norm>SpAMM_tolerance
      do11x10=qA%Quad11%Norm*qB%Quad10%Norm>SpAMM_tolerance
      do10x01=qA%Quad10%Norm*qB%Quad01%Norm>SpAMM_tolerance
      do11x11=qA%Quad11%Norm*qB%Quad11%Norm>SpAMM_tolerance
  
      ! #ifdef SpAMM_PRINT
  !!$    IF(.NOT.do00x00.AND. &
  !!$       .NOT.do01x10.AND. &
  !!$       .NOT.do00x01.AND. &
  !!$       .NOT.do01x11.AND. &
  !!$       .NOT.do10x00.AND. &
  !!$       .NOT.do11x10.AND. &
  !!$       .NOT.do10x01.AND. &
  !!$       .NOT.do11x11)THEN
  !!$
  !!$       WRITE(*,*)' Lev = ',qA%Lev
  
         SpAMM_stream%Box(1,:)=qA%Box(2,:) ! I
         SpAMM_stream%Box(2,:)=qA%Box(1,:) ! K
         SpAMM_stream%Box(3,:)=qB%Box(1,:) ! J
         ALLOCATE(SpAMM_stream%Next)
         SpAMM_stream=>SpAMM_stream%Next
         NULLIFY(SpAMM_stream%Next)
  !!    ENDIF
      ! #endif
  
  
      ! 00=00*00+01*10
      IF(do00x00) &
         CALL SpAMM_Multiply(qC%Quad00,qA%Quad00,qB%Quad00)
      IF(do01x10) &
         CALL SpAMM_Multiply(qC%Quad00,qA%Quad01,qB%Quad10)
      ! 01=00*01+01*11
      IF(do00x01) &
           CALL SpAMM_Multiply(qC%Quad01,qA%Quad00,qB%Quad01)
      IF(do01x11) &
              CALL SpAMM_Multiply(qC%Quad01,qA%Quad01,qB%Quad11)
      ! 10=10*00+11*10
      IF(do10x00) &
           CALL SpAMM_Multiply(qC%Quad10,qA%Quad10,qB%Quad00)
      IF(do11x10) &
           CALL SpAMM_Multiply(qC%Quad10,qA%Quad11,qB%Quad10)
      ! 11=10*01+11*11
      IF(do10x01) &
           CALL SpAMM_Multiply(qC%Quad11,qA%Quad10,qB%Quad01)
      IF(do11x11) &
           CALL SpAMM_Multiply(qC%Quad11,qA%Quad11,qB%Quad11)
      ! Bounds
      qC%Norm=qC%Quad00%Norm+qC%Quad01%Norm+qC%Quad10%Norm+qC%Quad11%Norm
      qC%Quad00%Norm=SQRT(qC%Quad00%Norm)
      qC%Quad10%Norm=SQRT(qC%Quad10%Norm)
      qC%Quad01%Norm=SQRT(qC%Quad01%Norm)
      qC%Quad11%Norm=SQRT(qC%Quad11%Norm)
      IF(qC%Lev==1)THEN
         qC%Norm=SQRT(qC%Norm)
      ENDIF
    END SUBROUTINE SpAMM_Multiply
  
  !!$
  !!$       FUNCTION RANDOM_INT(Limits)
  !!$          INTEGER                :: RANDOM_INT
  !!$          INTEGER, DIMENSION(2)  :: Limits
  !!$          REAL(DOUBLE)           :: Delta
  !!$          REAL(DOUBLE), EXTERNAL :: Random
  !!$!
  !!$          Delta=DBLE(Limits(2)-Limits(1)+1)
  !!$          RANDOM_INT=Limits(1)+INT(Delta*Random())
  !!$!
  !!$       END FUNCTION RANDOM_INT
  !!$
  !!$      FUNCTION RANDOM_DBL(Limits)
  !!$         REAL(DOUBLE)              :: RANDOM_DBL
  !!$         REAL(DOUBLE),DIMENSION(2) :: Limits
  !!$         REAL(DOUBLE)              :: Delta
  !!$         REAL(DOUBLE), EXTERNAL    :: Random
  !!$!
  !!$         Delta=Limits(2)-Limits(1)+0.0D0
  !!$         RANDOM_DBL=Limits(1)+Delta*Random()
  !!$!
  !!$      END FUNCTION RANDOM_DBL
  
  
  
  END MODULE SpAMM_TYPES
  
  
  PROGRAM SpAMM_TEST
    USE SpAMM_TYPES
    IMPLICIT NONE
    INTEGER           ::  L,M,N_OLD,N,I,J,K,II,JJ
    REAL(DOUBLE),ALLOCATABLE, DIMENSION(:,:) ::  A, B, C1,C2,CDiff,A_NOPADDING
    TYPE(QuTree),POINTER  :: qA,qB,qC
    TYPE(SpAMM_cubes),POINTER  :: Stream
    REAL(DOUBLE)          :: CNorm,CDiffNorm,CErrBound,CElementErr
    CHARACTER *100        :: Buffer
  
    CALL GETARG(1,BUFFER)
    READ(BUFFER,*)N
    CALL GETARG(2,BUFFER)
    READ(BUFFER,*)SpAMM_tolerance
  !  CALL GETARG(3,BUFFER)
  !  OPEN(UNIT=66,FILE=TRIM(ADJUSTL(BUFFER)))
  
  !  WRITE(*,*)N,' READING ',TRIM(ADJUSTL(BUFFER))
  
  
  !!$  ALLOCATE(A_NOPADDING(1:N,1:N))
  !!$
  !!$  DO I=1,N
  !!$     DO J=1,N
  !!$        READ(66,*)L,M,A_NOPADDING(L,M)
  !!$     ENDDO
  !!$  ENDDO
  !!$  !
  !!$
  !!$  K=CEILING(LOG10(DBLE(N))/LOG10(2D0))
  !!$  !
  !!$  N_OLD=N
  !!$  N=2**K
  
  
  
  
    !
    ALLOCATE(A(1:N,1:N))
    ALLOCATE(B(1:N,1:N))
    ALLOCATE(C1(1:N,1:N))
    ALLOCATE(C2(1:N,1:N))
    ALLOCATE(CDiff(1:N,1:N))
    !
    A=Zero
    B=Zero
    A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
    B(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
  
  
    CALL RANDOM_SEED
    CALL RANDOM_NUMBER(A)
    CALL RANDOM_NUMBER(B)
  
    DO I=1,N
       DO J=1,N
          !
          ! >>> COMMENTS: If you make "gamma" larger, the accounting for
          ! >>> number of multiplies becomes funny for very small thresholds
          !
  !        IF(I==J)THEN
          A(I,J)=EXP(-5.0D-1*ABS(I-J)) !*A(I,J)
          B(I,J)=EXP(-2.0D-1*ABS(I-J)) !*B(I,J)
  !        ELSE
  !           A(I,J)=0
  !           B(I,J)=0
  !        ENDIF
       ENDDO
    ENDDO
    !
    SpAMM_tiles=CEILING(DBLE(N)/BLOCK_SIZE)!*BLOCK_SIZE
    SpAMM_levels=FLOOR(LOG(DBLE(M))/LOG(2D0))
    !
    qA=>Dense2Quad(A)
    qB=>Dense2Quad(B)
    ! #ifdef SpAMM_print
    ALLOCATE(SpAMM_stream)
    Stream=>SpAMM_stream
    ! #ifdef
  
    qC=>SpAMM(qA,qB,SpAMM_tolerance)
  
  
    WRITE(44,*)'Graphics3D[{Opacity[0.4],'
    DO WHILE(ASSOCIATED(Stream%Next))
       Stream=>Stream%Next
       WRITE(44,111)Stream%Box(1,1),Stream%Box(2,1),Stream%Box(3,1),Stream%Box(1,2),Stream%Box(2,2),Stream%Box(3,2)
  111 FORMAT("Cuboid[{",I6,",",I6,",",I6,"},{",I6,",",I6,",",I6,"}], ")
    ENDDO
    WRITE(44,*)'}]'
    WRITE(*,*)' DONE SPAMM '
  
    STOP
  
    C2=Quad2Dense(qC)
    C1=MATMUL(A,B)
    CDiff=C1-C2
    CNorm=SQRT(SUM(C1(1:N,1:N)**2))
    CDiffNorm=SQRT(SUM(CDiff(1:N,1:N)**2))
    CErrBound=SpAMM_tolerance*SQRT(SUM(A(1:N,1:N)**2))*SQRT(SUM(B(1:N,1:N)**2))
  
    DO I=1,N
       DO J=1,N
          CDiff(I,J)=ABS(CDiff(I,J)) !/C1(I,J))
       ENDDO
    ENDDO
    CElementErr=MAXVAL(CDIFF)
  
    DO I=1,N
       DO J=1,N
          IF(CDiff(I,J)==CElementErr)THEN
             !           WRITE(*,*)I,J
             GOTO 101
          ENDIF
       ENDDO
    ENDDO
  101 CONTINUE
  
    WRITE(*,33)N,SpAMM_tolerance,CDiffNorm/CErrBound,CElementErr,SpAMM_multiplies/DBLE(SpAMM_tiles**3)
  33 FORMAT(I6,', ', 3(D12.6,', '),D12.6)
  
    ! WRITE(*,*)C1(I,J),C2(I,J),(C1(I,J)-C2(I,J))/C1(I,J)
  
    STOP
  
    !
    CALL DeleteQuad(qA)
    CALL DeleteQuad(qB)
    CALL DeleteQuad(qC)
    !
    DEALLOCATE(A)
    DEALLOCATE(B)
    DEALLOCATE(C1)
    DEALLOCATE(C2)
    DEALLOCATE(CDiff)
  
  END PROGRAM SpAMM_TEST
  
  
  
  !!$      ALLOCATE(GlblQT%Next)
  !!$      ALLOCATE(GlblQT%Next%Next)
  !!$      ALLOCATE(GlblQT%Next%Next%Next)
  !!$      ALLOCATE(GlblQT%Next%Next%Next%Next)
  !!$      !
  !!$      GlblQT=>GlblQT%Next
  !!$      GlblQT%Quad=>qA%Quad00
  !!$      GlblQT%Hash=Hash+1*Shift
  !!$      Hash00=GlblQT%Hash
  !!$      !
  !!$      GlblQT=>GlblQT%Next
  !!$      GlblQT%Quad=>qA%Quad01
  !!$      GlblQT%Hash=Hash+2*Shift
  !!$      Hash01=GlblQT%Hash
  !!$      !
  !!$      GlblQT=>GlblQT%Next
  !!$      GlblQT%Quad=>qA%Quad10
  !!$      GlblQT%Hash=Hash+3*Shift
  !!$      Hash10=GlblQT%Hash
  !!$      !
  !!$      GlblQT=>GlblQT%Next
  !!$      GlblQT%Quad=>qA%Quad11
  !!$      GlblQT%Hash=Hash+4*Shift
  !!$      Hash11=GlblQT%Hash
  !
  !!$!--------------------------------------------------------------
  !!$!    F90 wrapper for SFCOrder77, which circumvents the lack
  !!$!    of INTEGER(KIND=8) (INTEGER*8) support for cheazy
  !!$!    F90 compilers (pgf,nag...)
  !!$     SUBROUTINE SFCOrder(N,RD,PointI,SFC_KEY)
  !!$        INTEGER                       :: N
  !!$        INTEGER                       :: SFC_KEY
  !!$
  !!$!        TYPE(DBL_RNK2)                :: R
  !!$!        TYPE(INT_VECT)                :: Point
  !!$!        TYPE(DBL_VECT)                :: RKey
  !!$
  !!$        REAL(DOUBLE),DIMENSION(:),ALLOCATABLE :: RKeyD
  !!$        REAL(DOUBLE),DIMENSION(:,:),ALLOCATABLE :: RD
  !!$        INTEGER, DIMENSION(:),ALLOCATABLE :: PointI
  !!$        INTEGER(INT8),ALLOCATABLE, &
  !!$                         DIMENSION(:) :: IKey
  !!$!
  !!$
  !!$        INTEGER                       :: I
  !!$        INTEGER::N_Neighbor
  !!$        REAL(DOUBLE)::Gamma,Alpha,NLFac
  !!$
  !!$           ALLOCATE(IKey(N))
  !!$           CALL SFCOrder77(N,RD,PointI,IKey,.FALSE.)
  !!$           DEALLOCATE(IKey)
  !!$
  !!$     END SUBROUTINE SFCOrder
  !!$
  !!$     SUBROUTINE Sort_DBL_INT(X,Y,N_O,Ordr_O)
  !!$       REAL(DOUBLE),DIMENSION(:) :: XD
  !!$       INTEGER,DIMENSION(:) :: YI
  !!$
  !!$        TYPE(DBL_VECT)     :: X
  !!$        TYPE(INT_VECT)     :: Y
  !!$        INTEGER,OPTIONAL   :: N_O,Ordr_O
  !!$        INTEGER                       :: N,Ordr
  !!$        Ordr=-2
  !!$        N=MIN(SIZE(XD),SIZE(YI))
  !!$        IF(PRESENT(N_O))THEN
  !!$           IF(N_O>N)CALL Halt(' Dimensioning off in DblSort ')
  !!$           N=N_O
  !!$        ENDIF
  !!$        IF(PRESENT(Ordr_O))Ordr=Ordr_O
  !!$        CALL DblIntSort77(N,XD,YI,Ordr)
  !!$    END SUBROUTINE Sort_DBL_INT
  !!$
  !!$!--------------------------------------------------------------
  !!$!--------------------------------------------------------------
  !!$  SUBROUTINE SFCOrder77(N,R,Point,Key,Hilbert)
  !!$!    Use DerivedTypes
  !!$    IMPLICIT NONE
  !!$    INTEGER,PARAMETER::BitNum=21 ! number of bits to represent an integer
  !!$    INTEGER::IMax,I,J,N,Ix,Iy,Iz,Point(N)
  !!$    INTEGER(INT8):: Key(N)
  !!$    INTEGER(INT8),EXTERNAL:: Interleave,HilbertKey
  !!$    REAL(DOUBLE)::R(3,N),RMin(3),Ratio,MaxDiff,Diff
  !!$    LOGICAL::Hilbert
  !!$
  !!$    RMin(1)=1.D10; RMin(2)=1.D10; RMin(3)=1.D10
  !!$
  !!$    DO J=1,N
  !!$       Point(J)=J
  !!$       RMin(1)=DMIN1(RMin(1),R(1,J))
  !!$       RMin(2)=DMIN1(RMin(2),R(2,J))
  !!$       RMin(3)=DMIN1(RMin(3),R(3,J))
  !!$    ENDDO
  !!$
  !!$    MaxDiff = -1.0D0 ! any negative value will do
  !!$    DO J = 1,N
  !!$      DO I = 1, 3
  !!$        Diff = R(I,J) - RMin(I)
  !!$        IF( Diff > MaxDiff ) THEN
  !!$          MaxDiff = Diff
  !!$        ENDIF
  !!$      ENDDO
  !!$    ENDDO
  !!$
  !!$    IF( MaxDiff .LE. 0) THEN
  !!$      WRITE(*,*) 'MaxDiff = ',MaxDiff
  !!$      STOP 'ERROR: MaxDiff must be positive!'
  !!$    ENDIF
  !!$    Ratio = (2**BitNum-1.0)/MaxDiff
  !!$
  !!$    IMax = -Big_Int
  !!$    IF(Hilbert)THEN
  !!$       DO J=1,N
  !!$          Ix=DNINT((R(1,J)-RMin(1))*Ratio)
  !!$          IF(Ix > IMax) IMax = Ix
  !!$          Iy=DNINT((R(2,J)-RMin(2))*Ratio)
  !!$          IF(Iy > IMax) IMax = Iy
  !!$          Iz=DNINT((R(3,J)-RMin(3))*Ratio)
  !!$          IF(Iz > IMax) IMax = Iz
  !!$          Key(J)=Interleave(BitNum,Ix,Iy,Iz)
  !!$          Key(J)=HilbertKey(BitNum,Key(J))
  !!$       END DO
  !!$
  !!$       IF( IMax /= (2**BitNum-1)) THEN
  !!$!        stop 'ERROR: numerical accuracy problem !!'
  !!$       END IF
  !!$    ELSE ! Peano option
  !!$       DO J=1,N
  !!$          Ix=DNINT((R(1,J)-RMin(1))*Ratio)
  !!$          Iy=DNINT((R(2,J)-RMin(2))*Ratio)
  !!$          Iz=DNINT((R(3,J)-RMin(3))*Ratio)
  !!$          Key(J)=Interleave(BitNum,Ix,Iy,Iz)
  !!$       END DO
  !!$    ENDIF
  !!$    CALL I8Sort(Key,Point,N,2)
  !!$  END SUBROUTINE SFCOrder77
  !!$
  !!$!-----------------------------------------------------------------
  !!$!   Bit shuffle for a triple of integers
  !!$  FUNCTION Interleave(BitNum,Ix,Iy,Iz)
  !!$     USE DerivedTypes
  !!$     IMPLICIT NONE
  !!$     INTEGER I,K,Ix,Iy,Iz,BitNum,EndNum
  !!$     INTEGER(INT8) Interleave
  !!$     K=0
  !!$     Interleave=0
  !!$     EndNum = 3*(BitNum-1)
  !!$     DO I=0,EndNum,3
  !!$        IF(BTEST(Iz,K))Interleave=IBSET(Interleave,I  )
  !!$        IF(BTEST(Iy,K))Interleave=IBSET(Interleave,I+1)
  !!$        IF(BTEST(Ix,K))Interleave=IBSET(Interleave,I+2)
  !!$        K=K+1
  !!$     ENDDO
  !!$  END FUNCTION Interleave
  !!$
  !!$!-----------------------------------------------------------------
  !!$! Hilbert ordering for three-dimensional data, (3*BitNum)-bit
  !!$! implementation (stored in 64-bit integer).
  !!$! Interleaved key based on Algortithm H2 of Faloutsos and Roseman, PODS 89.
  !!$  FUNCTION HilbertKey(BitNum,Key)
  !!$!     USE DerivedTypes
  !!$     IMPLICIT NONE
  !!$     INTEGER I,J,K,L,BitNum,BeginNum
  !!$     INTEGER(INT8) Key
  !!$     INTEGER(INT8) HilbertKey
  !!$     INTEGER SubKey(21)
  !!$     INTEGER ToState(0:7,12)
  !!$     INTEGER ToBinry(0:7,12)
  !!$!--------------------------------------------------------------
  !!$!  State table, based on Fig 3 , T. Bially,
  !!$!  IEEE Trans. on Information Theory, 1969, pp. 658, vol. IT-15
  !!$!                                 0, 1, 2, 3, 4, 5, 6, 7
  !!$     DATA (ToState(I,1 ),I=0,7) / 9,12, 2, 1, 2, 3, 9, 1/
  !!$     DATA (ToBinry(I,1 ),I=0,7) / 0, 1, 3, 2, 7, 6, 4, 5/
  !!$     DATA (ToState(I,2 ),I=0,7) / 5, 3, 2, 2, 3, 5, 1, 4/
  !!$     DATA (ToBinry(I,2 ),I=0,7) / 4, 3, 5, 2, 7, 0, 6, 1/
  !!$     DATA (ToState(I,3 ),I=0,7) / 2, 3, 6, 3, 1, 7, 7, 1/
  !!$     DATA (ToBinry(I,3 ),I=0,7) / 6, 5, 1, 2, 7, 4, 0, 3/
  !!$     DATA (ToState(I,4 ),I=0,7) /10, 9, 4, 2, 5, 2, 4, 9/
  !!$     DATA (ToBinry(I,4 ),I=0,7) / 6, 7, 5, 4, 1, 0, 2, 3/
  !!$     DATA (ToState(I,5 ),I=0,7) / 5, 2, 5, 6, 8, 4, 4, 8/
  !!$     DATA (ToBinry(I,5 ),I=0,7) / 2, 1, 5, 6, 3, 0, 4, 7/
  !!$     DATA (ToState(I,6 ),I=0,7) / 6, 6, 5, 3, 7, 8, 3, 5/
  !!$     DATA (ToBinry(I,6 ),I=0,7) / 2, 5, 3, 4, 1, 6, 0, 7/
  !!$     DATA (ToState(I,7 ),I=0,7) / 6, 7,11,12,11, 7, 6, 3/
  !!$     DATA (ToBinry(I,7 ),I=0,7) / 4, 5, 7, 6, 3, 2, 0, 1/
  !!$     DATA (ToState(I,8 ),I=0,7) / 8, 6,10,11, 8,11, 5, 6/
  !!$     DATA (ToBinry(I,8 ),I=0,7) / 2, 3, 1, 0, 5, 4, 6, 7/
  !!$     DATA (ToState(I,9 ),I=0,7) /12,10, 1, 4,10,12, 9, 9/
  !!$     DATA (ToBinry(I,9 ),I=0,7) / 0, 7, 1, 6, 3, 4, 2, 5/
  !!$     DATA (ToState(I,10),I=0,7) / 8, 4, 4, 8,10, 9,10,11/
  !!$     DATA (ToBinry(I,10),I=0,7) / 4, 7, 3, 0, 5, 6, 2, 1/
  !!$     DATA (ToState(I,11),I=0,7) / 7, 8,12,10,11,11,10,12/
  !!$     DATA (ToBinry(I,11),I=0,7) / 6, 1, 7, 0, 5, 2, 4, 3/
  !!$     DATA (ToState(I,12),I=0,7) / 1, 7, 7, 1, 9,12,11,12/
  !!$     DATA (ToBinry(I,12),I=0,7) / 0, 3, 7, 4, 1, 2, 6, 5/
  !!$
  !!$! Form 3 bit SubKeys from left to right (Step 3 of H2)
  !!$     DO I=1,BitNum
  !!$        SubKey(I)=0
  !!$     ENDDO
  !!$     K=1
  !!$     BeginNum=3*(BitNum-1)
  !!$     DO I=BeginNum,0,-3
  !!$        IF(BTEST(Key,I  ))SubKey(K)=IBSET(SubKey(K),0)
  !!$        IF(BTEST(Key,I+1))SubKey(K)=IBSET(SubKey(K),1)
  !!$        IF(BTEST(Key,I+2))SubKey(K)=IBSET(SubKey(K),2)
  !!$        IF(K.GT.BitNum)STOP ' gt BitNum '
  !!$        IF(I.LT.0)STOP ' lt 0 '
  !!$        K=K+1
  !!$     ENDDO
  !!$! Change each 3 bit SubKey according to the output key
  !!$! from the current state and move to a new state (Step 4 of H2)
  !!$     I=1
  !!$     DO K=1,BitNum
  !!$        L=SubKey(K)
  !!$        SubKey(K)=ToBinry(L,I)
  !!$        I        =ToState(L,I)
  !!$     ENDDO
  !!$
  !!$! Reassemble key from sub keys, with SubKey(21)
  !!$! the rightmost 3 bits (Step 5 of H2)
  !!$     HilbertKey=0
  !!$     I=0
  !!$     DO K=BitNum,1,-1
  !!$        IF(BTEST(SubKey(K),0))HilbertKey=IBSET(HilbertKey,I  )
  !!$        IF(BTEST(SubKey(K),1))HilbertKey=IBSET(HilbertKey,I+1)
  !!$        IF(BTEST(SubKey(K),2))HilbertKey=IBSET(HilbertKey,I+2)
  !!$        I=I+3
  !!$     ENDDO
  !!$  END FUNCTION HilbertKey
  
negative.code: 
  !------------------------------------------------------------------------------
  !    This code is part of the MondoSCF suite of programs for linear scaling
  !    electronic structure theory and ab initio molecular dynamics.
  