      implicit real*8(a-h,o-z)
      integer ndim, n
      parameter (ndim=100)
      real*8 a(ndim*(ndim+1)/2), c(ndim*ndim), eigen(ndim)
c
c     read dimension of matrix to be solved
c
   10 continue
      read(5,'(i5)') n
      if( n .gt. ndim ) then
         write(6,'(a)')
     $        'dimension n is too big. Please retry to enter n'
         go to 10
      end if
c
      call makmat(n,a,c)
c
c     diagonarization matrix "a", matrix "c" is eigen vectors
c
      call erduw(a,c,n,1.0d-15)
c
      inum = 0
      do i = 1 , n
         do j = 1 , i
            inum = inum + 1
            if( i .eq. j ) eigen(i) = a(inum)
         end do
      end do
c
      call sort(n,eigen,c)
      call wrimat(n,eigen,c)
c
      end
c==================================================================
      subroutine makmat(n,a,c)
      implicit real*8(a-h,o-z)
      integer n
      real*8 a(*), c(*)
c
      zero=0.0d0
      one=1.0d0
c
      inum = 0
      do i = 1 , n
         do j = 1 , i
            inum = inum + 1
            a(inum) = zero
            if( abs(i-j) .eq. 1 ) then
               a(inum) = one
            end if
         end do
      end do
c
      inum = 0
      do i = 1 , n
         do j = 1 , n
            inum = inum + 1
            if( i .eq. j ) c(inum) = one
         end do
      end do
c
      return
      end
c==================================================================
      SUBROUTINE ERDUW(A,B,NA,EPSLON)
C TAKEN MAINLY FROM VERSION 4, AUG., 1971, OF JACSCF, U. OF WA.
C SQUARE-ROOT PROCEDURE IMPROVED.
C   MATRIX DIAGONALIZATION BY THE JACOBI METHOD.
C     A = REAL SYMMETRIC MATRIX TO BE DIAGONALIZED. IT IS STORED
C       BY COLUMNS WITH ALL SUB-DIAGONAL ELEMENTS OMITTED, SO
C       A(I,J) IS STORED AS A((J*(J-1))/2+I).
C     B = MATRIX TO BE MULTIPLIED BY THE MATRIX OF EIGENVECTORS.
C     NA = DIMENSION OF THE MATRICES.
C     EPSLON IS THE CONVERGENCE CRITERION FOR OFF-DIAGONAL ELEMENTS.
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(*), B(*)
      DATA A0, A1S2 /0.0D0, 0.5D0/
      IF(NA.EQ.1) RETURN
      LOOPC=128
      NN=(NA*(NA-1))/16
      AMAX=A0
      SUM=A0
      IJ=1
      DO 24 I=1,NA
      DO 16 J=1,I
      IF(I.EQ.J) GO TO 8
      AMAX=MAX(AMAX,ABS(A(IJ)))
    8 TERM=A(IJ)*A(IJ)
      SUM=SUM+TERM+TERM
   16 IJ=IJ+1
   24 SUM=SUM-TERM
      THRSHG=SQRT(SUM/NA)*EPSLON
      IF(THRSHG.GE.AMAX) RETURN
      THRESH=MAX(A1S2*AMAX,THRSHG)
   32 N=0
      IJ=2
      JJ=1
      DO 112 J=2,NA
      JJ=JJ+J
      JM1=J-1
      II=0
      DO 104 I=1,JM1
      II=II+I
      IF(ABS(A(IJ)).LT.THRESH) GO TO 104
      N=N+1
      SUM=A1S2*(A(JJ)+A(II))
      TERM=A1S2*(A(JJ)-A(II))
      AMAX=SIGN(SQRT(TERM*TERM+A(IJ)*A(IJ)),TERM)
      C=SQRT((AMAX+TERM)/(AMAX+AMAX))
      S=A(IJ)/(C*(AMAX+AMAX))
      A(II)=SUM-AMAX
      A(JJ)=SUM+AMAX
      A(IJ)=A0
      IM1=I-1
      IF(IM1) 40,56,40
   40 KI=II-I
      KJ=JJ-J
      DO 48 K=1,IM1
      KI=KI+1
      KJ=KJ+1
      TERM=C*A(KI)-S*A(KJ)
      A(KJ)=S*A(KI)+C*A(KJ)
   48 A(KI)=TERM
   56 IF(JM1.EQ.I) GO TO 72
      IP1=I+1
      IK=II+I
      KJ=IJ
      DO 64 K=IP1,JM1
      KJ=KJ+1
      TERM=C*A(IK)-S*A(KJ)
      A(KJ)=S*A(IK)+C*A(KJ)
      A(IK)=TERM
   64 IK=IK+K
   72 IF(J.EQ.NA) GO TO 88
      JP1=J+1
      IK=JJ+I
      JK=JJ+J
      DO 80 K=JP1,NA
      TERM=C*A(IK)-S*A(JK)
      A(JK)=S*A(IK)+C*A(JK)
      A(IK)=TERM
      IK=IK+K
   80 JK=JK+K
   88 KI=IM1*NA
      KJ=JM1*NA
      DO 96 K=1,NA
      KI=KI+1
      KJ=KJ+1
      TERM=C*B(KI)-S*B(KJ)
      B(KJ)=S*B(KI)+C*B(KJ)
   96 B(KI)=TERM
  104 IJ=IJ+1
  112 IJ=IJ+1
      LOOPC=LOOPC-1
      IF(LOOPC) 120,999,120
  120 IF(N.GT.NN) GO TO 32
      IF(THRESH.EQ.THRSHG) GO TO 128
      THRESH=MAX(A1S2*THRESH,THRSHG)
      GO TO 32
  128 IF(N) 32,136,32
  136 RETURN
C
C ********** ERROR CONDITIONS **********
C
  999 CONTINUE
      WRITE (6,2048)
      STOP 999
C
C END ERROR CONDITIONS
C
 2048 FORMAT(27H0 ERROR IN SUBROUTINE ERDUW)
      END
c==================================================================
      subroutine sort(n,eigen,c)
      implicit real*8(a-h,o-z)
      real*8 eigen(n), c(n,n)
c
      do i = 1 , n-1
         do j = i+1 , n
            IF( eigen(i) .lt. eigen(j) ) then
               dummy = eigen(i)
               eigen(i) = eigen(j)
               eigen(j) = dummy
               do k = 1 , n
                  dummy = c(k,i)
                  c(k,i) = c(k,j)
                  c(k,j) = dummy
               end do
            end if
         end do
      end do
c
      return
      end
c==================================================================
      subroutine wrimat(n,eigen,vec)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION eigen(n), vec(n,n)
c
      WRITE(6,'(a)') 'EIGEN VALUES AND EIGEN VECTORS'
      INDEX = (N-1)/8 + 1
      DO I = 1 , INDEX
         I1 = 8*(I-1) + 1
         I2 = 8*I
         IF( N .LT. I2 ) I2 = N
C
         WRITE(6,'(8f10.5)') ( eigen(K),K=I1,I2 )
         WRITE(6,*)
C
         DO J = 1 , N
            WRITE(6,'(8f10.5)') ( vec(j,k),K=I1,I2 )
         end do
         WRITE(6,*)
      end do
c
      return
      end
