anchor.diff:      IF (I.EQ.1) THEN
          TrP=Trace(P)
          TrdP=Trace(dP)
  +       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  +       !       TrP=TrdP
          IF(ABS(TrP-TrdP)>SpAMM_tolerance*1D1)THEN
             WRITE(*,*)' 1TrP = ',TrP,TrdP
             STOP
          ENDIF
  -    ENDIF
  -
  -    WRITE(*,*)' SpAMM SpAMM SpAMM SpAMM SpAMM SpAMM '
  +       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  +    ENDIF
  +    !
       CALL Multiply(P,P,P2)             ! The only multiplication is a square
  -    WRITE(*,*)' SpAMM SpAMM SpAMM SpAMM SpAMM SpAMM '
  -
       dP2=MATMUL(dP,dP)
  +
  +
  +!    P2=>Dense2Quad(Quad2Dense(P2))
  +    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       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)
  +
  +             sum=Zero
  +             DO KK=1,NN
  +                IF(ABS(dP(II,KK)*dP(KK,J))>1D-1*SpAMM_tolerance)THEN
  +                   WRITE(*,*)dP(II,KK)*dP(KK,J),sum
  +                   sum=sum+dP(II,KK)*dP(KK,J)
  +                ENDIF
  +             ENDDO
  +             WRITE(*,*)' SpAMM_tolerance = ',SpAMM_tolerance
  +            WRITE(*,*)'ARAW',II,J,dP2(II,J),dP2b(II,J),dP2(II,J)-dP2b(II,J)
  +            WRITE(*,*)'ASUM',II,J,sum,dP2b(II,J),dP2(II,J)-sum
  +
  +!             CALL Print_quad(P2)
  + !            P=>Dense2Quad(dP2)
                WRITE(*,*)'++++++++++++++++++++++++++++++++++++++++++++'
                WRITE(*,*)' dense P*P'
  -             CALL Print_quad(P)
  +
  +
  + !            CALL Print_quad(P)
                STOP
             ENDIF
          ENDDO
       ENDDO
  -
  +    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       TrP2=Trace(P2)
       TrdP2=Trace(dP2)
  -
  -    IF(ABS(TrP2-TrdP2)>SpAMM_tolerance*1D1)THEN
  +    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  +    !    TrP2=TrdP2
  +    IF(ABS(TrP2-TrdP2)>SpAMM_tolerance)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
  
positive.code:   INTEGER, PARAMETER :: INT8=SELECTED_INT_KIND(18) !--Integer*8
    !
    INTEGER, PARAMETER :: SINGLE=KIND(0.0)           !--Real*4
    INTEGER, PARAMETER :: DOUBLE=KIND(0.0D0)           !--Real*8
    REAL(DOUBLE),PARAMETER :: Zero=0D0,Half=5D-1,One=1D0,Two=2D0
    !
    TYPE QuTree
       INTEGER                      :: Siz
  
negative.code:   NULLIFY(qTmp2)
    CALL NewQuNode(qTmp1,init=.TRUE.)
    CALL NewQuNode(qTmp2,init=.TRUE.)
    !--------------------------------------------------
    SpAMM_average_multiplies=Zero
    IF(.NOT.DoFilter)THEN
       SpAMM_tolerance=SpAMM_tolerance*SQRT(Two/DBLE(Nel))
    ENDIF
    !--------------------------------------------------
    Occ0 = 0.D0;Occ1 = 0.D0;Occ2 = 0.D0;Occ3 = 0.D0
    DO I=1,TC2_cycles
       IF(DoFilter)THEN
          CALL TC2_filter(qP,qTmp1,qTmp2,Half*FLOAT(NEl),Occ0,I)
       ELSE
          CALL TC2(qP,qTmp1,qTmp2,Half*FLOAT(NEl),Occ0,I)
       ENDIF
       Occ3 = Occ2;Occ2 = Occ1;Occ1 = Occ0
       SpAMM_average_multiplies=SpAMM_average_multiplies+SpAMM_multiplies
    ENDDO
    !
    IF(DoFilter)THEN
       CALL Filter_qutree(qF,SpAMM_tolerance)
       CALL Multiply_qutree_tms_qutree(qP,qF,qTmp1,tolerance=1D-100)
    ELSE
       CALL Multiply(qP,qF,qTmp1)
    ENDIF
    TrE=Trace(qTmp1)
    RelativeErrorE=ABS(TrE-TargetTrE)/ABS(TargetTrE)
    RelativeErrorN=ABS(Occ0-Half*DBLE(Nel))/(Half*DBLE(Nel))
    !
    SpAMM_average_multiplies=SpAMM_average_multiplies+SpAMM_multiplies
    SpAMM_average_multiplies=SpAMM_average_multiplies/DBLE(TC2_cycles+1)
    !
    WRITE(*,66)N_old,N,SpAMM_BLOCK_SIZE,SpAMM_tolerance,RelativeErrorE,RelativeErrorN, &
               SpAMM_average_multiplies
  66 FORMAT(I6,", ",I6,", ",I3,", ",E8.2,7(", ",E18.8))
    !
  END PROGRAM SpAMM_TEST
  