module jacobi
implicit none
contains

function jacobi_evd(a,v) result(sweeps)
real(8) a(:,:),v(:,:),theta,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; 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),a(q,q)-a(p,p))
	c=cos(theta); s=sin(theta)
	a1pp=c*c*a(p,p)-2*s*c*a(p,q)+s*s*a(q,q)
	a1qq=s*s*a(p,p)+2*s*c*a(p,q)+c*c*a(q,q)
	if(a1pp.ne.a(p,p) .or. a1qq.ne.a(q,q))then
		changed=1
		xp=c*a(p,:)-s*a(q,:)
		xq=s*a(p,:)+c*a(q,:)
		a(p,:)=xp
		a(q,:)=xq
		xp=c*a(:,p)-s*a(:,q)
		xq=s*a(:,p)+c*a(:,q)
		a(:,p)=xp
		a(:,q)=xq
		xp=c*v(:,p)-s*v(:,q)
		xq=s*v(:,p)+c*v(:,q)
		v(:,p)=xp
		v(:,q)=xq
	end if
	end do; end do
end do
end function jacobi_evd

end module jacobi