anchor.diff:    CALL NewQuNode(qTmp1,init=.TRUE.)
     CALL NewQuNode(qTmp2,init=.TRUE.)
     !
  +
  +!!$ ./a.out bc  bc 256  296 1D-7
  +!!$BLOCK_SIZE=1
  +!!$1D-4, Tr(P)=148.00410675, O(N)/N^3 =  0.03703, TrE =   -453.64326697675500
  +!!$1D-6, Tr(P)=148.00007063, O(N)/N^3 =  0.24725, TrE =   -453.49036394380886
  +!!$1D-7, Tr(P)=148.00003758, O(N)/N^3 =  0.34967, TrE =   -453.48904088853374
  +!!$1D-8  Tr(P)=148.00000463, O(N)/N^3 =  0.52270, TrE =   -453.48864186234346
  +!!$1D-10 Tr(P)=148.00000006, O(N)/N^3 =  0.81716, TrE =   -453.48859793980563
  +
  +!!$
  +!!$BLOCK_SIZE=2 (DOUBLE)
  +!!$1D-4, Tr(P)=147.99964372, O(N)/N^3 =  0.11211, TrE =   -453.53013549760510
  +!!$1D-5, Tr(P)=148.00000991, O(N)/N^3 =  0.24183, TrE =   -453.49417884535342
  +!!$1D-6, Tr(P)=148.00000677, O(N)/N^3 =  0.42200, TrE =   -453.48909106857525
  +!!$1D-7, Tr(P)=148.00000486, O(N)/N^3 =  0.61572, TrE =   -453.48864613262896
  +!!$1D-8, Tr(P)=148.00000090, O(N)/N^3 =  0.77632, TrE =   -453.48860167043756
  +!!$1D-10,Tr(P)=148.00000000, O(N)/N^3 =  0.95261, TrE =   -453.48859749814483
  +!!$BLOCK_SIZE=2 (SINGLE)
  +!!$1D-4, Tr(P)=147.99964905, O(N)/N^3 =  0.11212, TrE =   -453.53015
  +!!$1D-5, Tr(P)=148.00001526, O(N)/N^3 =  0.24183, TrE =   -453.49417
  +!!$1D-6, Tr(P)=148.00000000, O(N)/N^3 =  0.42211, TrE =   -453.48907
  +!!$1D-7, Tr(P)=148.00001526, O(N)/N^3 =  0.61580, TrE =   -453.48868
  +!!$1D-8, Tr(P)=148.00001526, O(N)/N^3 =  0.77637, TrE =   -453.48865
  +!!$1D-10,Tr(P)=148.00000000, O(N)/N^3 =  0.95261, TrE =   -453.48865
  +
  +!!$BLOCK_SIZE=256
  +!!$      Tr(P)=148.00000000, O(N)/N^3 =  1.00000, TrE =   -453.48859746697627
  +
  +
     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)
  +  DO I=1,40
  +    CALL TC2(qP,qTmp1,qTmp2,Half*FLOAT(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)
  + CALL Multiply(qP,qF,qTmp1)
  + TrE=Trace(qTmp1)
  + WRITE(*,*)' TrE = ',TrE
   
   END PROGRAM SpAMM_TEST
  
positive.code:     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
  
negative.code: 
  !    WRITE(*,*)'TRACE , Lev =',qA%Lev,' Num = ',qA%Num,ASSOCIATED(qA%Quad00),ASSOCIATED(qA%Quad11)
  
      IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
         trc=Zero
         IF(.NOT.ALLOCATED(qA%Blok))RETURN
         DO I=1,SpAMM_BLOCK_SIZE
            trc=trc+qA%Blok(I,I)
         ENDDO
  !       write(*,*)'Num = ',qA%Num,'Bx = ',qA%Box(1,1),' trc = ',trc
      ELSEIF(.NOT.(ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qA%Quad11)))THEN
         Trc=Zero
         WRITE(*,*)' Z '
     ELSEIF(.NOT.ASSOCIATED(qA%Quad11))THEN
  !       WRITE(*,*)' 00 ',ASSOCIATED(qA%Quad00)
  !       WRITE(*,*)' 00 ',ASSOCIATED(qA%Quad00),ASSOCIATED(qA%Quad11)
         trc=Trace_qutree(qA%Quad00)
  !       WRITE(*,*)' 00 ',trc
      ELSEIF(.NOT.ASSOCIATED(qA%Quad00))THEN
         trc=Trace_qutree(qA%Quad11)
  !       WRITE(*,*)' 11 ',trc
      ELSE
         trc=Trace_qutree(qA%Quad00)+Trace_qutree(qA%Quad11)
  !       WRITE(*,*)' 00+11 ',trc
      ENDIF
    END FUNCTION Trace_qutree
    !=================================================================
  