anchor.diff:    !
     TYPE SpAMM_cubes
        TYPE(SpAMM_cubes), POINTER  :: Next
  +     INTEGER                     :: Lev
        INTEGER                     :: Hash
        REAL(DOUBLE)                :: Bound
        INTEGER,     DIMENSION(3,2) :: Box
  
positive.code:     IF(ASSOCIATED(qA%Quad11)) &
           CALL Print_Quad(qA%Quad11)
  
  111 FORMAT("Cuboid[{",I6,",",I6,":",I6,",",I6,"}, (* ",I6,", ",I6,", ",F12.6,"*)")
  
    END SUBROUTINE Print_Quad
  
    SUBROUTINE NewQuNode(qA,count)
      INTEGER :: count
      TYPE(QuTree), POINTER :: qA
      ALLOCATE(qA);  count=count+1
      qA%Num=count
      qA%Norm=Zero
      NULLIFY(qA%Quad00)
      NULLIFY(qA%Quad01)
      NULLIFY(qA%Quad10)
      NULLIFY(qA%Quad11)
    END SUBROUTINE NewQuNode
  
  
    FUNCTION SpAMM(qA,qB,tolerance) RESULT(qC)
      IMPLICIT NONE
      TYPE(QuTree), POINTER                :: qA,qB
      TYPE(QuTree),  POINTER      :: qC
      REAL(DOUBLE),OPTIONAL       :: tolerance
      !
      IF(PRESENT(tolerance))THEN
         SpAMM_tolerance=tolerance
      ELSE
         SpAMM_tolerance=1D-6
      ENDIF
      !
      SpAMM_quadcount=0
      CALL NewQuNode(qC,SpAMM_quadcount)
      !
      qC%Num=1
      qC%Lev=1
      SpAMM_multiplies=Zero
      qC%Box(:,1)=(/1,1/)
      qC%Box(:,2)=(/SpAMM_tiles,SpAMM_tiles/)*BLOCK_SIZE
      !
      CALL SpAMM_Multiply(qC,qA,qB)
    END FUNCTION SpAMM
  
  
    RECURSIVE SUBROUTINE SpAMM_Multiply(qC,qA,qB)
      IMPLICIT NONE
      TYPE(QuTree), POINTER :: qC,qA,qB
      REAL(DOUBLE), DIMENSION(1:BLOCK_SIZE,1:BLOCK_SIZE) :: CTmp
      LOGICAL :: do00x00,do01x10,do10x00,do11x10, &
                 do00x01,do01x11,do10x01,do11x11
      ! Bounds
      qC%Box(1,1)=qA%Box(2,1)
      qC%Box(1,2)=qA%Box(2,2)
      qC%Box(2,1)=qB%Box(1,1)
  
negative.code:        qC%Blok=qC%Blok+CTmp
         ! Bounds
         qC%Norm=qC%Norm+SUM(CTmp(1:BLOCK_SIZE,1:BLOCK_SIZE)**2)
      ELSE
         qC%Norm=Zero
         ! SpAMM=Sparse Approximate Matrix-Matrix
         do00x00=ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qB%Quad00)
         do01x10=ASSOCIATED(qA%Quad01).AND.ASSOCIATED(qB%Quad10)
         do00x01=ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qB%Quad01)
         do01x11=ASSOCIATED(qA%Quad01).AND.ASSOCIATED(qB%Quad11)
         do10x00=ASSOCIATED(qA%Quad10).AND.ASSOCIATED(qB%Quad00)
         do11x10=ASSOCIATED(qA%Quad11).AND.ASSOCIATED(qB%Quad10)
         do10x01=ASSOCIATED(qA%Quad10).AND.ASSOCIATED(qB%Quad01)
         do11x11=ASSOCIATED(qA%Quad11).AND.ASSOCIATED(qB%Quad11)
         ! SpAMM, SpAMM and SpAMM
         IF(do00x00)do00x00=do00x00.AND.(qA%Quad00%Norm*qB%Quad00%Norm>SpAMM_tolerance)
         IF(do01x10)do01x10=do01x10.AND.(qA%Quad01%Norm*qB%Quad10%Norm>SpAMM_tolerance)
         IF(do00x01)do00x01=do00x01.AND.(qA%Quad00%Norm*qB%Quad01%Norm>SpAMM_tolerance)
         IF(do01x11)do01x11=do01x11.AND.(qA%Quad01%Norm*qB%Quad11%Norm>SpAMM_tolerance)
         IF(do10x00)do10x00=do10x00.AND.(qA%Quad10%Norm*qB%Quad00%Norm>SpAMM_tolerance)
         IF(do11x10)do11x10=do11x10.AND.(qA%Quad11%Norm*qB%Quad10%Norm>SpAMM_tolerance)
         IF(do10x01)do10x01=do10x01.AND.(qA%Quad10%Norm*qB%Quad01%Norm>SpAMM_tolerance)
         IF(do11x11)do11x11=do11x11.AND.(qA%Quad11%Norm*qB%Quad11%Norm>SpAMM_tolerance)
         ! 00=00*00+01*10
         IF(do00x00.OR.do01x10)THEN
            IF(.NOT.ASSOCIATED(qC%Quad00))CALL NewQuNode(qC%Quad00,SpAMM_quadcount)
            IF(do00x00)CALL SpAMM_Multiply(qC%Quad00,qA%Quad00,qB%Quad00)
            IF(do01x10)CALL SpAMM_Multiply(qC%Quad00,qA%Quad01,qB%Quad10)
            qC%Norm=qC%Norm+qC%Quad00%Norm
            qC%Quad00%Norm=SQRT(qC%Quad00%Norm)
         ELSEIF(ASSOCIATED(qC%Quad00))THEN
            CALL DeleteQuad(qC%Quad00)
         ENDIF
         ! 01=00*01+01*11
         IF(do00x01.OR.do01x11)THEN
            IF(.NOT.ASSOCIATED(qC%Quad01))CALL NewQuNode(qC%Quad01,SpAMM_quadcount)
            IF(do00x01)CALL SpAMM_Multiply(qC%Quad01,qA%Quad00,qB%Quad01)
            IF(do01x11)CALL SpAMM_Multiply(qC%Quad01,qA%Quad01,qB%Quad11)
            qC%Norm=qC%Norm+qC%Quad01%Norm
            qC%Quad01%Norm=SQRT(qC%Quad01%Norm)
         ELSEIF(ASSOCIATED(qC%Quad01))THEN
            CALL DeleteQuad(qC%Quad01)
         ENDIF
         ! 10=10*00+11*10
         IF(do10x00.OR.do11x10)THEN
            IF(.NOT.ASSOCIATED(qC%Quad10))CALL NewQuNode(qC%Quad10,SpAMM_quadcount)
            IF(do10x00)CALL SpAMM_Multiply(qC%Quad10,qA%Quad10,qB%Quad00)
            IF(do11x10)CALL SpAMM_Multiply(qC%Quad10,qA%Quad11,qB%Quad10)
            qC%Norm=qC%Norm+qC%Quad10%Norm
            qC%Quad10%Norm=SQRT(qC%Quad10%Norm)
         ELSEIF(ASSOCIATED(qC%Quad10))THEN
            CALL DeleteQuad(qC%Quad10)
         ENDIF
         ! 11=10*01+11*11
         IF(do10x01.OR.do11x11)THEN
            IF(.NOT.ASSOCIATED(qC%Quad11))CALL NewQuNode(qC%Quad11,SpAMM_quadcount)
            IF(do10x01)CALL SpAMM_Multiply(qC%Quad11,qA%Quad10,qB%Quad01)
            IF(do11x11)CALL SpAMM_Multiply(qC%Quad11,qA%Quad11,qB%Quad11)
            qC%Norm=qC%Norm+qC%Quad11%Norm
            qC%Quad11%Norm=SQRT(qC%Quad11%Norm)
         ELSEIF(ASSOCIATED(qC%Quad11))THEN
            CALL DeleteQuad(qC%Quad11)
         ENDIF
         IF(qC%Lev==1)THEN
            qC%Norm=SQRT(qC%Norm)
         ENDIF
      ENDIF
  
    END SUBROUTINE SpAMM_Multiply
  
       SUBROUTINE TC2(P,P2,Tmp1,Norm,TrP,I)
         TYPE(QuTree),POINTER :: P,P2,Tmp1
         REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2
         INTEGER      :: I
         !-------------------------------------------------------------------------------
         IF (I.EQ.1) THEN
            TrP=Trace(P)
         ENDIF
         WRITE(*,*)'1'
         CALL Multiply(P,P,P2)             ! The only multiplication is a square
         WRITE(*,*)'2'
         TrP2=Trace(P2)
         WRITE(*,*)'3'
         CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
         CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria
         WRITE(*,33)I,NORM,TrP,CR1,CR2
  33     FORMAT(I4,", ",F10.5,", ",F16.8,", ",F10.5,", ",F10.5)
         IF (CR1 < CR2) THEN               ! Too many states
            P=>P2
         ELSE
            CALL Multiply(P,Two)
         WRITE(*,*)'4'
            CALL Multiply(P2,-One)
         WRITE(*,*)'5'
            CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
         WRITE(*,*)'6'
            P=>Tmp1
         ENDIF
         WRITE(*,*)'7'
         TrP=Trace(P)
         WRITE(*,*)'8'
       END SUBROUTINE TC2
  
     END MODULE SpAMM_TYPES
  
  PROGRAM SpAMM_TEST
    USE SpAMM_TYPES
    IMPLICIT NONE
    INTEGER           ::  L,M,N_OLD,N,Nel,I,J,K,II,JJ
    REAL(DOUBLE),ALLOCATABLE, DIMENSION(:,:) ::  A, B, C,C2,CDiff,A_NOPADDING
    TYPE(QuTree),POINTER  :: qC,qP,qF,qTmp1,qTmp2
    TYPE(SpAMM_cubes),POINTER  :: Stream
    REAL(DOUBLE)          :: CNorm,CDiffNorm,CErrBound,CElementErr,Opacity,Occ0,Occ1,Occ2,Occ3
    CHARACTER *100        :: Name,Buffer
  
    CALL GETARG(1,NAME)
    NAME=TRIM(ADJUSTL(NAME))
    CALL GETARG(2,BUFFER)
    READ(BUFFER,*)N
    CALL GETARG(3,BUFFER)
    READ(BUFFER,*)Nel
    CALL GETARG(4,BUFFER)
    READ(BUFFER,*)SpAMM_tolerance
    !
    K=CEILING(LOG10(DBLE(N))/LOG10(2D0))
    N_OLD=N
    N=2**K
    SpAMM_tiles=CEILING(DBLE(N)/BLOCK_SIZE)!*BLOCK_SIZE
    SpAMM_levels=CEILING(LOG(DBLE(SpAMM_tiles))/LOG(2D0))+1
    !
    ALLOCATE(A_NOPADDING(1:N_OLD,1:N_OLD))
    !--------------------------------------------------
    OPEN(UNIT=66,FILE=TRIM(NAME)//".f_ortho")
    DO I=1,N_OLD
       DO J=1,N_OLD
          READ(66,*)L,M,A_NOPADDING(L,M)
       ENDDO
    ENDDO
    CLOSE(66)
    ALLOCATE(A(1:N,1:N))
    A=Zero
    A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
    qF=>Dense2Quad(A)
    CALL Print_Quad(qF)
    !--------------------------------------------------
    OPEN(UNIT=66,FILE=TRIM(NAME)//".p_guess_ortho")
    DO I=1,N_OLD
       DO J=1,N_OLD
          READ(66,*)L,M,A_NOPADDING(L,M)
       ENDDO
    ENDDO
    CLOSE(66)
    A=Zero
    A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
    qP=>Dense2Quad(A)
    CALL Print_Quad(qP)
    !--------------------------------------------------
  !  DEALLOCATE(A)
    DEALLOCATE(A_NOPADDING)
    !
    ALLOCATE(qTmp1)
    ALLOCATE(qTmp2)
  !!$
  !!$
  !!$  ALLOCATE(C(1:N,1:N))
  !!$  ALLOCATE(C2(1:N,1:N))
  !!$
  !!$  CALL Multiply(qP,qP,1D-30)
  !!$  C=Quad2Dense(qC)
  !!$  !
  !!$  C2=MATMUL(A,A)
  !!$
  !!$  DO I=1,4
  !!$     DO J=1,4
  !!$        WRITE(*,*)C(I,J),C2(I,J)
  !!$     ENDDO
  !!$  ENDDO
  !!$
  !!$!  STOP
    !
    Occ0 = 0.D0
    Occ1 = 0.D0
    Occ2 = 0.D0
    Occ3 = 0.D0
    DO I=1,10
      CALL TC2(qP,qTmp1,qTmp2,Half*DBLE(NEl),Occ0,I)
      Occ3 = Occ2
      Occ2 = Occ1
      Occ1 = Occ0
   ENDDO
  
  
  !!$#ifdef SpAMM_print
  !!$  ALLOCATE(SpAMM_stream)
  !!$  Stream=>SpAMM_stream
  !!$#ifdef
  !!$  qC=>SpAMM(qA,qB,SpAMM_tolerance)
  !!$  WRITE(44,*)'Graphics3D[{Yellow,'
  !!$  DO WHILE(ASSOCIATED(Stream%Next))
  !!$     Stream=>Stream%Next
  !!$     IF(Stream%Lev==SpAMM_levels)THEN
  !!$        Opacity=1D0
  !!$     ELSE
  !!$        Opacity=0.4D0*DBLE(Stream%Lev)/DBLE(SpAMM_levels)
  !!$     ENDIF
  !!$     WRITE(44,111)Opacity,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("Opacity[",F10.5,"],Cuboid[{",I6,",",I6,",",I6,"},{",I6,",",I6,",",I6,"}], ")
  !!$  ENDDO
  !!$
  !!$  WRITE(44,*)'Boxed->False,ViewAngle -> Automatic, ViewCenter -> {0.5, 0.5, 0.5}, '
  !!$  WRITE(44,*)'ViewMatrix -> Automatic, ViewPoint -> {-1.17551, -1.22153, -2.92849},'
  !!$  WRITE(44,*)'ViewRange -> All, ViewVector -> Automatic,'
  !!$  WRITE(44,*)'ViewVertical -> {-0.961419, -0.126619, -0.244213}]'
  !!$  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
  