c_jacobi Subroutine

subroutine c_jacobi(A,D,U,sweep)
  implicit none
  complex(8),dimension(:,:),intent(inout)               :: a
  real(8),dimension(size(a,1)),intent(out)              :: d
  complex(8),dimension(size(a,1),size(a,2)),intent(out) :: U
  integer                                               :: n
  integer                                               :: p, q, j
  real(8)                                               :: red, off, thresh
  real(8)                                               :: t, delta, invc, s
  complex(8)                                            :: x, y, Apq
  real(8)                                               :: ev(2,size(a,1))
  integer                                               :: sweep
  real(8)                                               :: SYM_EPS=tiny(1d0)
  !
  n=size(a,1);if(size(a,2)/=n)stop "Error in Jacobi: size(a)!=n**2 - a not a square matrix"
  !
  do p = 1, n
     ev(1,p) = 0d0
     ev(2,p) = dble(A(p,p))
     d(p) = ev(2,p)
  enddo
  U = eye(n)
  red = 0.04d0/n**4
  do sweep = 1, 50
     off = 0
     do q = 2, n
        do p = 1, q - 1
           off = off + abs(A(p,q))**2
        enddo
     enddo
     if( .not. off .gt. SYM_EPS ) return !goto 10
     thresh = 0
     if( sweep .lt. 4 ) thresh = off*red
     do q = 2, n
        do p = 1, q - 1
           Apq = A(p,q)
           off = abs(Apq)**2
           if( sweep .gt. 4 .and. off .lt. SYM_EPS*(ev(2,p)**2 + ev(2,q)**2) ) then
              A(p,q) = 0
           else if( off .gt. thresh ) then
              t = .5D0*(ev(2,p) - ev(2,q))
              t = 1/(t + sign(sqrt(t**2 + off), t))
              delta = t*off
              ev(1,p) = ev(1,p) + delta
              ev(2,p) = d(p) + ev(1,p)
              ev(1,q) = ev(1,q) - delta
              ev(2,q) = d(q) + ev(1,q)
              invc = sqrt(delta*t + 1)
              s = t/invc
              t = delta/(invc + 1)
              do j = 1, p - 1
                 x = A(j,p)
                 y = A(j,q)
                 A(j,p) = x + s*(conjg(Apq)*y - t*x)
                 A(j,q) = y - s*(Apq*x + t*y)
              enddo
              do j = p + 1, q - 1
                 x = A(p,j)
                 y = A(j,q)
                 A(p,j) = x + s*(Apq*conjg(y) - t*x)
                 A(j,q) = y - s*(Apq*conjg(x) + t*y)
              enddo
              do j = q + 1, n
                 x = A(p,j)
                 y = A(q,j)
                 A(p,j) = x + s*(Apq*y - t*x)
                 A(q,j) = y - s*(conjg(Apq)*x + t*y)
              enddo
              A(p,q) = 0d0
              do j = 1, n
                 x = U(p,j)
                 y = U(q,j)
                 ! #if UCOLS
                 ! U(p,j) = x + s*(conjg(Apq)*y - t*x)
                 ! U(q,j) = y - s*(Apq*x + t*y)
                 ! #else
                 U(p,j) = x + s*(Apq*y - t*x)
                 U(q,j) = y - s*(conjg(Apq)*x + t*y)
                 ! #endif
              enddo
           endif
        enddo
     enddo
     do p = 1, n
        ev(1,p) = 0
        d(p) = ev(2,p)
     enddo
  enddo
  stop "Too many iteration in Jacobi"
  ! 10  continue
  !     ! if( sort .eq. 0 ) return
  !     do p = 1, n - 1
  !        j = p
  !        t = d(p)
  !        do q = p + 1, n
  !           if( 1*(t - d(q)) .gt. 0 ) then
  !              j = q
  !              t = d(q)
  !           endif
  !        enddo
  !        if( j .ne. p ) then
  !           d(j) = d(p)
  !           d(p) = t
  !           do q = 1, n
  !              x = U(p,q)
  !              U(p,q) = U(j,q)
  !              U(j,q) = x
  !           enddo
  !        endif
  !     enddo
end subroutine C_jacobi