anchor.diff:    A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
     qP=>Dense2Quad(A)
     !--------------------------------------------------
  -!  DEALLOCATE(A)
  +  DEALLOCATE(A)
     DEALLOCATE(A_NOPADDING)
     !
     NULLIFY(qTmp1)
  
positive.code:        ALLOCATE(qA)
         qA%Num=1
         qA%Lev=1
         qA%Siz=SpAMM_tiles*BLOCK_SIZE
         qA%Box(:,1)=(/1,1/)
         qA%Box(:,2)=(/SpAMM_tiles,SpAMM_tiles/)*BLOCK_SIZE
      ELSE
         IF(.NOT.ASSOCIATED(qA))THEN
            ALLOCATE(qA)
  
negative.code:     IMPLICIT NONE
      TYPE(QuTree), POINTER       :: qA,qB,qC
      REAL(DOUBLE),OPTIONAL       :: tolerance
      REAL(DOUBLE)                :: Saved_tolerance,norm
      INTEGER                     :: SpAMM_count
      IF(PRESENT(tolerance))THEN
         Saved_tolerance=SpAMM_tolerance
         SpAMM_tolerance=tolerance
      ENDIF
      CALL DeleteQuad(qC)
      CALL NewQuNode(qC,init=.TRUE.)
      SpAMM_count=1
      SpAMM_multiplies=0
      CALL SpAMM_Multiply(qC,qA,qB,count=SpAMM_count)
      qC%Norm=SQRT(SpAMM_NormReduce(qC))
      IF(PRESENT(tolerance))THEN
         SpAMM_tolerance=Saved_tolerance
      ENDIF
    END SUBROUTINE Multiply_qutree_tms_qutree
    !
    RECURSIVE FUNCTION SpAMM_NormReduce(qC) RESULT(norm)
      IMPLICIT NONE
      TYPE(QuTree), POINTER :: qC
      REAL(DOUBLE) :: norm
      ! Associated
      IF(.NOT.ASSOCIATED(qC))THEN
         norm=Zero
         RETURN
      ELSEIF(qC%Siz==SpAMM_BLOCK_SIZE.AND.ALLOCATED(qC%Blok))THEN
         norm=SUM(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE)**2)
         qC%Norm=SQRT(norm)
      ELSE
         norm=SpAMM_NormReduce(qC%Quad00) + &
              SpAMM_NormReduce(qC%Quad01) + &
              SpAMM_NormReduce(qC%Quad10) + &
              SpAMM_NormReduce(qC%Quad11)
         qC%Norm=SQRT(norm)
      ENDIF
    END FUNCTION SpAMM_NormReduce
  
    RECURSIVE SUBROUTINE SpAMM_Multiply(qC,qA,qB,count)
      IMPLICIT NONE
      TYPE(QuTree), POINTER :: qC,qA,qB
      INTEGER :: count
      REAL(DOUBLE) :: norm
      REAL(DOUBLE), DIMENSION(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE) :: CTmp
      ! Associated
      IF(ASSOCIATED(qA).AND.ASSOCIATED(qB))THEN
      ! Bounds
      IF(qA%Norm*qB%Norm<SpAMM_tolerance)RETURN
      ! Allocate
      IF(.NOT.ASSOCIATED(qC))CALL NewQuNode(qC,count=count)
      ! Boxing
       qC%Siz=qA%Siz
       qC%Lev=qA%Lev
       qC%Box(1,:)=qA%Box(1,:)
       qC%Box(2,:)=qB%Box(2,:)
         ! Blocks
         IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
            ! Allocate
            IF(.NOT.ALLOCATED(qC%Blok))THEN
               ALLOCATE(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE))
               qC%Blok=Zero
            END IF
            ! Count
            SpAMM_multiplies=SpAMM_multiplies+1
            ! Accumulate
            qC%Blok=qC%Blok+MATMUL(qA%Blok,qB%Blok)
         ELSE
            ! 00=00*00+01*10
            CALL SpAMM_Multiply(qC%Quad00,qA%Quad00,qB%Quad00,count=count)
            CALL SpAMM_Multiply(qC%Quad00,qA%Quad01,qB%Quad10,count=count)
            ! 01=00*01+01*11
            CALL SpAMM_Multiply(qC%Quad01,qA%Quad00,qB%Quad01,count=count)
            CALL SpAMM_Multiply(qC%Quad01,qA%Quad01,qB%Quad11,count=count)
            ! 10=10*00+11*10
            CALL SpAMM_Multiply(qC%Quad10,qA%Quad10,qB%Quad00,count=count)
            CALL SpAMM_Multiply(qC%Quad10,qA%Quad11,qB%Quad10,count=count)
            ! 11=10*01+11*11
            CALL SpAMM_Multiply(qC%Quad11,qA%Quad10,qB%Quad01,count=count)
            CALL SpAMM_Multiply(qC%Quad11,qA%Quad11,qB%Quad11,count=count)
         ENDIF
      ENDIF
    END SUBROUTINE SpAMM_Multiply
  
    FUNCTION Trace_dns(A) RESULT(trc)
      REAL(DOUBLE) :: trc
      REAL(DOUBLE),DIMENSION(:,:) :: A
      INTEGER I
      trc=Zero
      DO I=1,SIZE(A,1)
         trc=trc+A(I,I)
         !       WRITE(*,*)I,' TRACE ',A(I,I)
      ENDDO
    END FUNCTION Trace_dns
  
    SUBROUTINE TC2_test(P,P2,Tmp1,Norm,TrP,I)
      TYPE(QuTree),POINTER :: P,P2,Tmp1 , Tmp2
      REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2,    TrdP,TrdP2
      INTEGER      :: I,  II,J
  
      INTEGER, PARAMETER:: NN=8
      REAL(DOUBLE),DIMENSION(NN,NN) :: dP,dP2,dTmp1, dTmp2,dP2b
      !-------------------------------------------------------------------------------
      dP=Quad2Dense(P)
  
      IF (I.EQ.1) THEN
         TrP=Trace(P)
         TrdP=Trace(dP)
         IF(ABS(TrP-TrdP)>SpAMM_tolerance*1D1)THEN
            WRITE(*,*)' 1TrP = ',TrP,TrdP
            STOP
         ENDIF
      ENDIF
  
      CALL Multiply(P,P,P2)             ! The only multiplication is a square
  
      dP2=MATMUL(dP,dP)
      dP2b=Quad2Dense(P2)
  
      DO II=1,NN
         DO J=1,NN
            IF(ABS(dP2(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
              WRITE(*,*)'A',II,J,dP2(II,J),dP2b(II,J),dP2(II,J)-dP2b(II,J)
               WRITE(*,*)'++++++++++++++++++++++++++++++++++++++++++++'
               WRITE(*,*)' sparse P*P'
               CALL Print_quad(P2)
               P=>Dense2Quad(dP2)
               WRITE(*,*)'++++++++++++++++++++++++++++++++++++++++++++'
               WRITE(*,*)' dense P*P'
               CALL Print_quad(P)
               STOP
            ENDIF
         ENDDO
      ENDDO
  
      TrP2=Trace(P2)
      TrdP2=Trace(dP2)
  
      IF(ABS(TrP2-TrdP2)>SpAMM_tolerance*1D1)THEN
         WRITE(*,*)'TrP = ',TrP2,TrdP2
         STOP "P2 STOP "
      ENDIF
  
      CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
      CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria
  
      WRITE(*,33)I,NORM,TrP,SpAMM_multiplies/SpAMM_dimension**3
  33  FORMAT(I4,", N =",F8.5,", Tr(P)=",F8.5,", O(N)/N^3 = ",F8.5)
  
      IF (CR1 < CR2) THEN               ! Too many states
  
        P=>Copy(P2)
        dP2b=Quad2Dense(P)
  
        dP=dP2
         DO II=1,NN
            DO J=1,NN
               IF(ABS(dP(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
                  WRITE(*,*)' B ',II,J,dP(II,J)-dP2b(II,J)
                  STOP
               ENDIF
            ENDDO
         ENDDO
      ELSE
         CALL Multiply(P,Two)
         CALL Multiply(P2,-One)
         CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
         dP2b=Quad2Dense(Tmp1)
  
         P=>Copy(Tmp1)
         dP=Two*dP-dP2
  
         DO II=1,NN
            DO J=1,NN
               IF(ABS(dP(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
                  WRITE(*,*)II,J,dP(II,J)-dP2b(II,J)
                  STOP
               ENDIF
            ENDDO
         ENDDO
  
      ENDIF
      TrP=Trace(P)
      TrdP=Trace(dP)
  
      IF(ABS(TrP-TrdP)>SpAMM_tolerance*1D1)THEN
         WRITE(*,*)'TrP = ',TrP,TrdP
         STOP
      ENDIF
  
  
    END SUBROUTINE TC2_TEST
  
    SUBROUTINE TC2(P,P2,Tmp1,Norm,TrP,I)
      TYPE(QuTree),POINTER :: P,P2,Tmp1 , Tmp2
      REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2,    TrdP,TrdP2
      INTEGER      :: I,  II,J
      !-------------------------------------------------------------------------------
      IF (I.EQ.1)TrP=Trace(P)
      CALL Multiply(P,P,P2)             ! The only multiplication is a square
      TrP2=Trace(P2)
      CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
      CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria
      WRITE(*,33)I,NORM,TrP,SpAMM_multiplies/SpAMM_tiles**3
  33  FORMAT(I4,", N =",F12.8,", Tr(P)=",F12.8,", O(N)/N^3 = ",F8.5)
      IF (CR1 < CR2) THEN               ! Too many states
        P=>Copy(P2)
      ELSE
         CALL Multiply(P,Two)
         CALL Multiply(P2,-One)
         CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
         P=>Copy(Tmp1)
      ENDIF
      TrP=Trace(P)
    END SUBROUTINE TC2
  
  END MODULE SpAMM_TYPES
  
  PROGRAM SpAMM_TEST
    USE SpAMM_TYPES
  