anchor.diff:    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))
   
  -!  WRITE(*,*)LOG10(DBLE(N))/LOG10(2D0),SpAMM_tolerance,CDiffNorm/CErrBound,SpAMM_multiplies/DBLE(N**3)
  -!  WRITE(*,*)SpAMM_tolerance,CDiffNorm/CErrBound,SpAMM_multiplies/DBLE(N**3)
  -  WRITE(*,*)SpAMM_tolerance,CDiffNorm,CErrBound,SpAMM_multiplies/DBLE(N**3)
  +  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)
  
positive.code:         ! >>> COMMENTS: If you make "gamma" larger, the accounting for
          ! >>> number of multiplies becomes funny for very small thresholds
          !
          A(I,J)=EXP(-1.0D0*ABS(I-J)) *A(I,J)
          B(I,J)=EXP(-5.0D0*ABS(I-J)) *B(I,J)
       ENDDO
    ENDDO
    !
    SpAMM_tiles=BLOCK_SIZE*CEILING(DBLE(N)/BLOCK_SIZE)
    SpAMM_levels=FLOOR(LOG(DBLE(M))/LOG(2D0))
    !
    qA=>Dense2Quad(A)
  
negative.code:     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
  