module jacobi
implicit none
contains

function jacobi_evd(a,e,v) result(sweeps)
real(8) a(:,:),v(:,:),e(:)
real(8) theta,g,h,app,aqq,apq,a1pp,a1qq,c,s,xp(size(a,1)),xq(size(a,1))
integer i,n,changed,sweeps,p,q
n=size(a,1)
v(:,:)=0; do i=1,n; v(i,i)=1; e(i)=a(i,i); end do
changed=1
sweeps=0
do while (changed.eq.1)
changed=0
sweeps=sweeps+1
do p=1,n; do q=p+1,n
	theta=0.5*atan2(2*a(p,q),e(q)-e(p))
	c=cos(theta); s=sin(theta)
	a1pp=c*c*e(p)-2*s*c*a(p,q)+s*s*e(q)
	a1qq=s*s*e(p)+2*s*c*a(p,q)+c*c*e(q)
	if(a1pp.ne.e(p) .or. a1qq.ne.e(q))then
		changed=1
		e(p)=a1pp
		e(q)=a1qq
		a(p,q)=0
		do i=1,p-1
			g=a(i,p)
			h=a(i,q)
			a(i,p)=c*g-s*h
			a(i,q)=s*g+c*h
		end do
		do i=p+1,q-1
			g=a(p,i)
			h=a(i,q)
			a(p,i)=c*g-s*h
			a(i,q)=s*g+c*h
		end do
		do i=q+1,n
			g=a(p,i)
			h=a(q,i)
			a(p,i)=c*g-s*h
			a(q,i)=s*g+c*h
		end do
		do i=1,n
			g=v(i,p)
			h=v(i,q)
			v(i,p)=c*g-s*h
			v(i,q)=s*g+c*h
		end do
	end if
end do; end do
end do
end function jacobi_evd

end module jacobi