d_jacobi Subroutine

subroutine d_jacobi(a,d,v,nrot)
  real(8),dimension(:,:),intent(inout)               :: a
  real(8),dimension(size(a,1)),intent(out)           :: d
  real(8),dimension(size(a,1),size(a,2)),intent(out) :: v
  integer,intent(out)                                :: nrot
  integer                                            :: i,ip,iq,n
  real(8)                                            :: c,g,h,s,sm,t,tau,theta,tresh
  real(8),dimension(size(d))                         :: b,z
  !
  n=size(a,1);if(size(a,2)/=n)stop "Error in Jacobi: size(a)!=n**2 - a not a square matrix"
  !
  v = eye(n)
  forall(i=1:n)b(i)=a(i,i)
  !
  d=b
  z=0d0
  nrot=0
  !
  do i=1,50
     sm=sum(abs(a),mask=upper_triangle(n,n))
     if (sm == 0.0) return
     tresh=merge(0.2d0*sm/n**2, 0d0, i<4 )
     do ip=1,n-1
        do iq=ip+1,n
           g=100d0*abs(a(ip,iq))
           if ((i > 4) .and. (abs(d(ip))+g == abs(d(ip))) &
                .and. (abs(d(iq))+g == abs(d(iq)))) then
              a(ip,iq)=0.0
           else if (abs(a(ip,iq)) > tresh) then
              h=d(iq)-d(ip)
              if (abs(h)+g == abs(h)) then
                 t=a(ip,iq)/h
              else
                 theta=0.5d0*h/a(ip,iq)
                 t=1.0d0/(abs(theta)+sqrt(1.0d0+theta**2))
                 if (theta < 0.0) t=-t
              end if
              c=1d0/sqrt(1+t**2)
              s=t*c
              tau=s/(1d0+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0d0
              call jrotate(a(1:ip-1,ip),a(1:ip-1,iq))
              call jrotate(a(ip,ip+1:iq-1),a(ip+1:iq-1,iq))
              call jrotate(a(ip,iq+1:n),a(iq,iq+1:n))
              call jrotate(v(:,ip),v(:,iq))
              nrot=nrot+1
           end if
        end do
     end do
     b(:)=b(:)+z(:)
     d(:)=b(:)
     z(:)=0d0
  end do
  stop "Too many iteration in Jacobi"
contains
  subroutine jrotate(a1,a2)
    real(8),dimension(:), intent(inout) :: a1,a2
    real(8),dimension(size(a1))         :: wk1
    wk1(:)=a1(:)
    a1(:)=a1(:)-s*(a2(:)+a1(:)*tau)
    a2(:)=a2(:)+s*(wk1(:)-a2(:)*tau)
  end subroutine jrotate
end subroutine d_jacobi