anchor.diff:      TYPE(QuTree),POINTER :: qA
       CALL SpAMM_DeleteQuad(qA)
       DEALLOCATE(qA)
  +    NULLIFY(qA)
     END SUBROUTINE DeleteQuad
   
     RECURSIVE SUBROUTINE SpAMM_DeleteQuad(qA)
  
positive.code:     IF(ALLOCATED(qA%Blok))THEN
         DEALLOCATE(qA%Blok,STAT=Status)
         NULLIFY(qA%Blok)
         RETURN
      END IF
      IF(ASSOCIATED(qA%Quad00))THEN
         CALL SpAMM_DeleteQuad(qA%Quad00)
         DEALLOCATE(qA%Quad00)
         NULLIFY(qA%Quad00)
      ENDIF
      IF(ASSOCIATED(qA%Quad01))THEN
         CALL SpAMM_DeleteQuad(qA%Quad01)
         DEALLOCATE(qA%Quad01)
         NULLIFY(qA%Quad01)
      ENDIF
      IF(ASSOCIATED(qA%Quad10))THEN
         CALL SpAMM_DeleteQuad(qA%Quad10)
         DEALLOCATE(qA%Quad10)
         NULLIFY(qA%Quad10)
      ENDIF
      IF(ASSOCIATED(qA%Quad11))THEN
         CALL SpAMM_DeleteQuad(qA%Quad11)
         DEALLOCATE(qA%Quad11)
         NULLIFY(qA%Quad11)
      ENDIF
    END SUBROUTINE SpAMM_DeleteQuad
  
    RECURSIVE SUBROUTINE Print_Quad(qA)
      TYPE(QuTree),POINTER :: qA
      IF(ASSOCIATED(qA))THEN
         IF(.NOT.ASSOCIATED(qA%Quad00))RETURN
         WRITE(*,111)qA%Box(:,1),qA%Box(:,2),qA%Num,qA%Siz,qA%Norm
        ENDIF
      IF(ASSOCIATED(qA%Quad00))&
           WRITE(*,111)qA%Quad00%Box(:,1),qA%Quad00%Box(:,2),qA%Quad00%Num,qA%Quad00%Siz,qA%Quad00%Norm
      IF(ASSOCIATED(qA%Quad10))&
  
negative.code:     IF(ASSOCIATED(qA%Quad11)) &
           CALL Print_Quad(qA%Quad11)
  
  111 FORMAT("Cuboid[{",I8,",",I8,":",I8,",",I8,"}, (* ",I6,", ",I6,", ",F12.6,"*)")
  
    END SUBROUTINE Print_Quad
  
    !=================================================================
    ! QuTree trace routines: a=trace[A]
    !=================================================================
    RECURSIVE FUNCTION Trace_qutree(qA) RESULT(Trc)
      IMPLICIT NONE
      TYPE(QuTree), POINTER  :: qA
      REAL(DOUBLE) :: trc
      INTEGER :: I
      WRITE(*,*)ASSOCIATED(qA%Quad00),ASSOCIATED(qA%Quad11)
  
      IF(.NOT.(ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qA%Quad11)))THEN
         Trc=Zero
         RETURN
      ELSEIF(qA%Siz==BLOCK_SIZE)THEN
         IF(.NOT.ALLOCATED(qA%Blok))THEN
            WRITE(*,*)'ASSOCIATED ',ASSOCIATED(qA)
            WRITE(*,*)' QA NORM = ',qA%Norm
            STOP 'hosed'
         ENDIF
         trc=Zero
         DO I=1,BLOCK_SIZE
            trc=trc+qA%Blok(I,I)
         ENDDO
      ELSEIF(.NOT.ASSOCIATED(qA%Quad11))THEN
         trc=Trace_qutree(qA%Quad00)
      ELSEIF(.NOT.ASSOCIATED(qA%Quad00))THEN
         trc=Trace_qutree(qA%Quad11)
      ELSE
         trc=Trace_qutree(qA%Quad00)+Trace_qutree(qA%Quad11)
      ENDIF
    END FUNCTION Trace_qutree
    !=================================================================
    ! QuTree filter routines: B=filter[A,tau]
    !=================================================================
    RECURSIVE SUBROUTINE Filter_qutree(qA,tau)
      IMPLICIT NONE
      TYPE(QuTree), POINTER  :: qA
      REAL(DOUBLE)           :: tau
      IF(qA%Norm<tau)THEN
         CALL DeleteQuad(qA)
      ELSE
         CALL Filter_qutree(qA%Quad00,tau)
         CALL Filter_qutree(qA%Quad01,tau)
         CALL Filter_qutree(qA%Quad10,tau)
         CALL Filter_qutree(qA%Quad11,tau)
      ENDIF
    END SUBROUTINE Filter_qutree
    !=================================================================
    ! QuTree add routines: C=A+B
    !=================================================================
    RECURSIVE SUBROUTINE Add_qutree_pls_qutree(qA,qB,qC)
      IMPLICIT NONE
      TYPE(QuTree), POINTER    :: qA,qB,qC
      !
      IF(.NOT.ASSOCIATED(qA))THEN
         CALL Copy_qutree_eq_qutree(qB,qC)
      ELSEIF(.NOT.ASSOCIATED(qB))THEN
         CALL Copy_qutree_eq_qutree(qA,qC)
      ELSE
         IF(qA%Siz==BLOCK_SIZE)THEN
            IF(.NOT.ALLOCATED(qC%Blok))THEN
               ALLOCATE(qC%Blok(1:BLOCK_SIZE,1:BLOCK_SIZE))
               qC%Norm=0D0
               qC%Blok=0D0
               qC%Lev=qA%Lev
               qC%Siz=BLOCK_SIZE
            ENDIF
            ! Add
            qC%Blok=qA%Blok+qB%Blok
            ! Bounds
            qC%Norm=SUM(qC%Blok(1:BLOCK_SIZE,1:BLOCK_SIZE)**2)
         ELSE
            CALL NewQuNode(qC%Quad00)
            CALL NewQuNode(qC%Quad01)
            CALL NewQuNode(qC%Quad10)
            CALL NewQuNode(qC%Quad11)
            CALL Add_qutree_pls_qutree(qA%Quad00,qB%Quad00,qC%Quad00)
            CALL Add_qutree_pls_qutree(qA%Quad01,qB%Quad01,qC%Quad01)
            CALL Add_qutree_pls_qutree(qA%Quad10,qB%Quad10,qC%Quad10)
            CALL Add_qutree_pls_qutree(qA%Quad11,qB%Quad11,qC%Quad11)
            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
         ENDIF
      ENDIF
    END SUBROUTINE Add_qutree_pls_qutree
    !=================================================================
    ! QuTree copy routines: C=A
    !=================================================================
    RECURSIVE SUBROUTINE Copy_qutree_eq_qutree(qA,qC)
      IMPLICIT NONE
      TYPE(QuTree), POINTER    :: qA,qC
      !
      IF(.NOT.ASSOCIATED(qA))RETURN
      IF(.NOT.ASSOCIATED(qC))ALLOCATE(qC)
      IF(qA%Siz==BLOCK_SIZE)THEN
         IF(.NOT.ALLOCATED(qC%Blok))THEN
            ALLOCATE(qC%Blok(1:BLOCK_SIZE,1:BLOCK_SIZE))
            qC%Norm=qA%Norm
            qC%Lev=qA%Lev
            qC%Siz=BLOCK_SIZE
            qC%Blok=qA%Blok
         ELSE
            ALLOCATE(qC%Quad00)
            ALLOCATE(qC%Quad01)
            ALLOCATE(qC%Quad10)
            ALLOCATE(qC%Quad11)
            CALL Copy_qutree_eq_qutree(qA%Quad00,qC%Quad00)
            CALL Copy_qutree_eq_qutree(qA%Quad01,qC%Quad01)
            CALL Copy_qutree_eq_qutree(qA%Quad10,qC%Quad10)
            CALL Copy_qutree_eq_qutree(qA%Quad11,qC%Quad11)
         ENDIF
      ENDIF
    END SUBROUTINE Copy_qutree_eq_qutree
    !=================================================================
    ! QuTree multiply routines
    !=================================================================
    ! A=A*B
    !-----------------------------------------------------------------
    RECURSIVE SUBROUTINE Multiply_qutree_tms_scalar(qA,B)
      IMPLICIT NONE
      TYPE(QuTree), POINTER    :: qA
      REAL(DOUBLE)             :: B
      IF(.NOT.ASSOCIATED(qA))RETURN
      IF(qA%Siz==BLOCK_SIZE)THEN
         qA%Norm=qA%Norm*B
         qA%Blok=qA%Blok*B
      ELSE
         CALL Multiply_qutree_tms_scalar(qA%Quad00,B)
         CALL Multiply_qutree_tms_scalar(qA%Quad01,B)
         CALL Multiply_qutree_tms_scalar(qA%Quad10,B)
         CALL Multiply_qutree_tms_scalar(qA%Quad11,B)
         qA%Norm=qA%Norm*B
      ENDIF
    END SUBROUTINE Multiply_qutree_tms_scalar
    !-----------------------------------------------------------------
    ! C=A*B
    !-----------------------------------------------------------------
    SUBROUTINE Multiply_qutree_tms_qutree(qA,qB,qC,tolerance)
      IMPLICIT NONE
      TYPE(QuTree), POINTER       :: qA,qB,qC
      REAL(DOUBLE),OPTIONAL       :: tolerance
      REAL(DOUBLE)                :: Saved_tolerance
      IF(PRESENT(tolerance))THEN
         Saved_tolerance=SpAMM_tolerance
         SpAMM_tolerance=tolerance
      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)
      IF(PRESENT(tolerance))THEN
         SpAMM_tolerance=Saved_tolerance
      ENDIF
    END SUBROUTINE Multiply_qutree_tms_qutree
  !
    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%Siz=qA%Siz
      qC%Box(1,1)=qA%Box(2,1)
      qC%Box(1,2)=qA%Box(2,2)
      qC%Box(2,1)=qB%Box(1,1)
  