anchor.diff:  PROGRAM SpAMM_TEST
     USE SpAMM_TYPES
     IMPLICIT NONE
  -  INTEGER           ::  N,I,J,M,II,JJ
  -  REAL(DOUBLE),ALLOCATABLE, DIMENSION(:,:) ::  A, B, C1,C2,CDiff
  +  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
   
  
positive.code:   INTEGER, PARAMETER :: DOUBLE=KIND(0.D0)          !--Real*8
    REAL(DOUBLE),PARAMETER :: Zero=0D0
    !
    !----------------------------------------------------------------
    REAL(DOUBLE)       :: SpAMM_tolerance,SpAMM_multiplies
    INTEGER            :: SpAMM_tiles,    SpAMM_levels,     SpAMM_quadcount
    INTEGER, PARAMETER :: BLOCK_SIZE=2
    !----------------------------------------------------------------
    TYPE QuTree
       INTEGER                      :: Siz
       INTEGER                      :: Lev     ! Level of this box
  
negative.code: !!$      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
  