anchor.diff:        !
         !     If the fsum of off-diagonal terms is zero returns
         !
  -      if ( fsum < 1.d-10 ) return
  -
  +      if ( fsum < 1.d-10 ) then
  +        eigval(1:3) = eigval(1:3)*factor
  +        return
  +      endif
         !
         !     Performs the sweep in three rotations. One per off diagonal term
         !
  
positive.code:   !> Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using
    !> Jacobi iteration adapted from numerical recpies
    subroutine eigen3 (tensor, eigval, princ)
      real(kind=kreal) :: tensor(6)     !< tensor
      real(kind=kreal) :: eigval(3)     !< vector containing the eigvalches
      real(kind=kreal) :: princ(3, 3)   !< matrix containing the three principal column vectors
  
      integer, parameter :: msweep = 50
      integer :: i,j, is, ip, iq, ir
      real(kind=kreal) :: fsum, od, theta, t, c, s, tau, g, h, hd, btens(3,3)
  
      btens(1,1)=tensor(1); btens(2,2)=tensor(2); btens(3,3)=tensor(3)
      btens(1,2)=tensor(4); btens(2,1)=btens(1,2)
  
negative.code:     !
      !     Initialise princ to the identity
      !
      factor = 0.d0
      do i = 1, 3
        do j = 1, 3
          princ (i, j) = 0.d0
        end do
        princ (i, i) = 1.d0
        eigval (i) = btens (i, i)
        factor = factor + dabs(btens(i,i))
      end do
      !     Scaling and iszero/isnan exception
      if( factor == 0.d0 .or. factor /= factor ) then
        return
      else
        eigval(1:3) = eigval(1:3)/factor
        btens(1:3,1:3) = btens(1:3,1:3)/factor
      end if
  
      !
      !     Starts sweeping.
      !
  