      subroutine jacobi(a,u,im,n,eigen,eps1,eps2,eps3,itmax)
      dimension a(10,10),u(10,10),aik(10),eigen(10)
       nm1=n-1
      do 10 i=im,n
      do 10 j=im,n
   10 u(i,j)=0.
      do 20 i=im,n
   20 u(i,i) =1.
      if(im.ne.n) go to 30
      eigen(n)=a(n,n)
      return
   30 continue
       sigma1=0.0
      do 40 i=im,n
   40  sigma1=sigma1+a(i,i)**2
       do 160  iter=1,itmax
      do 130 i=im,nm1
       ip1=i+1
        do 130  j=ip1,n
       q= abs(a(i,i)-a(j,j))
       if( q .le. eps1) goto 50
       if(  abs(a(i,j)) .le. eps2)   goto 130
       p=2.0*a(i,j)*q/(a(i,i)-a(j,j))
       spq=  sqrt(p*p+q*q)
       csa=  sqrt((1.0+q/spq)/2.0)
       sna= p/(2.0*csa*spq)
       goto 60
   50   csa=1.0/sqrt(2.0)
       sna=csa
   60  continue
      do 70 k=im,n
      holdki=u(k,i)
      u(k,i)= holdki*csa+u(k,j)*sna
   70  u(k,j)= holdki*sna-u(k,j)*csa
      do 100 k=i,n
      if( k .gt. j) goto 90
      aik(k)=a(i,k)
      a(i,k)=csa*aik(k)+sna*a(k,j)
      if( k .ne. j)  goto 80
      a(j,k)=sna*aik(k)-csa*a(j,k)
   80  goto 100
   90   holdik=a(i,k)
      a(i,k)=csa*holdik+sna*a(j,k)
      a(j,k)=sna*holdik-csa*a(j,k)
  100  continue
      aik(j)=sna*aik(i)-csa*aik(j)
      do 120 k=im,j
      if( k .le. i) goto 110
      a(k,j)=sna*aik(k)-csa*a(k,j)
      goto 120
  110   holdki=a(k,i)
      a(k,i)=csa*holdki+sna*a(k,j)
      a(k,j)=sna*holdki-csa*a(k,j)
  120   continue
  130   a(i,j)=0.0
      sigma2=0.0
      do 140 i=im,n
      eigen(i)=a(i,i)
  140  sigma2=sigma2+eigen(i)**2
      eps = 1.0 - sigma1/sigma2
  150 format(1h ,'iter=',i5,'  jacobi eps=',e12.4)
      !write(6,150) iter,eps
      if(eps.lt.eps3) go to 170
  160  sigma1=sigma2
  170  return
      end
