c     ..................................................................
c
c        subroutine deigen
c
c        purpose
c           compute eigenvalues and eigenvectors of a real symmetric
c           matrix
c
c        usage
c           call deigen(a,r,n,mv)
c
c        description of parameters
c           a - original matrix (symmetric), destroyed in computation.
c               resultant eigenvalues are developed in diagonal of
c               matrix a in descending order.
c           r - resultant matrix of eigenvectors (stored columnwise,
c               in same sequence as eigenvalues)
c           n - order of matrices a and r
c           mv- input code
c                   0   compute eigenvalues and eigenvectors
c                   1   compute eigenvalues only (r need not be
c                       dimensioned but must still appear in calling
c                       sequence)
c
c        remarks
c           original matrix a must be real symmetric (storage mode=1)
c           matrix a cannot be in the same location as matrix r
c
c        subroutines and function subprograms required
c           none
c
c        method
c           diagonalization method originated by jacobi and adapted
c           by von neumann for large computers as found in 'mathematical
c           methods for digital computers', edited by a. ralston and
c           h.s. wilf, john wiley and sons, new york, 1962, chapter 7
c
c     ..................................................................
c
      subroutine deigen (a,r,n,mv)
      dimension a(1), r(1)
c
c        ...............................................................
c
c        if a double precision version of this routine is desired, the
c        c in column 1 should be removed from the double precision
c        statement which follows.
c
      double precision a, r, anorm, anrmx, thr, x, y, sinx, sinx2, cosx,
     1 cosx2, sincs, range
c
c        the c must also be removed from double precision statements
c        appearing in other routines used in conjunction with this
c        routine.
c
c        the double precision version of this subroutine must also
c        contain double precision fortran functions.  sqrt in statements
c        40, 68, 75, and 78 must be changed to dsqrt.  abs in statement
c        62 must be changed to dabs. the constant in statement 5 should
c        be changed to 1.0d-12.
c
c        ...............................................................
c
c        generate identity matrix
c
      range=1.0d-12
      if (mv-1) 10,40,10
   10 iq=-n
      do 30 j=1,n
      iq=iq+n
      do 30 i=1,n
      ij=iq+i
      r(ij)=0.0
      if (i-j) 30,20,30
   20 r(ij)=1.0d0
   30 continue
c
c        compute initial and final norms (anorm and anormx)
c
   40 anorm=0.0d0
      do 60 i=1,n
      do 60 j=i,n
      if (i-j) 50,60,50
   50 ia=i+(j*j-j)/2
      anorm=anorm+a(ia)*a(ia)
   60 continue
      if (anorm) 320,320,70
   70 anorm=1.414213562d0*dsqrt(anorm)
      anrmx=anorm*range/n
c
c        initialize indicators and compute threshold, thr
c
      ind=0
      thr=anorm
   80 thr=thr/n
   90 l=1
  100 m=l+1
c
c        compute sin and cos
c
  110 mq=(m*m-m)/2
      lq=(l*l-l)/2
      lm=l+mq
      if (dabs(a(lm))-thr) 250,120,120
  120 ind=1
      ll=l+lq
      mm=m+mq
      x=.5d0*(a(ll)-a(mm))
      y=-a(lm)/dsqrt(a(lm)*a(lm)+x*x)
      if (x) 130,140,140
  130 y=-y
  140 sinx=y/dsqrt(2.d0*(1.d0+(dsqrt(1.d0-y*y))))
      sinx2=sinx*sinx
      cosx=dsqrt(1.d0-sinx2)
      cosx2=cosx*cosx
      sincs=sinx*cosx
c
c        rotate l and m columns
c
      ilq=n*(l-1)
      imq=n*(m-1)
      do 240 i=1,n
      iq=(i*i-i)/2
      if (i-l) 150,220,150
  150 if (i-m) 160,220,170
  160 im=i+mq
      go to 180
  170 im=m+iq
  180 if (i-l) 190,200,200
  190 il=i+lq
      go to 210
  200 il=l+iq
  210 x=a(il)*cosx-a(im)*sinx
      a(im)=a(il)*sinx+a(im)*cosx
      a(il)=x
  220 if (mv-1) 230,240,230
  230 ilr=ilq+i
      imr=imq+i
      x=r(ilr)*cosx-r(imr)*sinx
      r(imr)=r(ilr)*sinx+r(imr)*cosx
      r(ilr)=x
  240 continue
      x=2.d0*a(lm)*sincs
      y=a(ll)*cosx2+a(mm)*sinx2-x
      x=a(ll)*sinx2+a(mm)*cosx2+x
      a(lm)=(a(ll)-a(mm))*sincs+a(lm)*(cosx2-sinx2)
      a(ll)=y
      a(mm)=x
c
c        tests for completion
c
c        test for m = last column
c
  250 if (m-n) 260,270,260
  260 m=m+1
      go to 110
c
c        test for l = second from last column
c
  270 if (l-(n-1)) 280,290,280
  280 l=l+1
      go to 100
  290 if (ind-1) 310,300,310
  300 ind=0
      go to 90
c
c        compare threshold with final norm
c
  310 if (thr-anrmx) 320,320,80
c
c        sort eigenvalues and eigenvectors
c
  320 iq=-n
      do 360 i=1,n
      iq=iq+n
      ll=i+(i*i-i)/2
      jq=n*(i-2)
      do 360 j=i,n
      jq=jq+n
      mm=j+(j*j-j)/2
      if (a(ll)-a(mm)) 330,360,360
  330 x=a(ll)
      a(ll)=a(mm)
      a(mm)=x
      if (mv-1) 340,360,340
  340 do 350 k=1,n
      ilr=iq+k
      imr=jq+k
      x=r(ilr)
      r(ilr)=r(imr)
  350 r(imr)=x
  360 continue
      return
      end