C
C      SUBROUTINE FOR DIAGONALISATION of A REAL SYMMETRIC MATRIX 
C      (single precision arithmatic only!)
C
       SUBROUTINE DIAG(N,A,U,P)
C
C Adapted from H.H GREENWOOD,"COMPUTING METHODS IN QUANTUM ORG.CHEMISTRY"
C A is the input matrix, real symmetric; P contains the eigen values.
C Columns of U are the eigen vectors. For eg., P(2) is the eigen value for
C the eigen vector {U(2,j)j=1,n}.
C mvr-2011. 
C
      DIMENSION UTEST(N),A(N,N),U(N,N),P(N)
      PARAMETER(EP=1.0D-08)
      DO 1 J=1,N
      DO 2 I=1,N
  2   A(I,J)=A(J,I)
  1   P(J)=A(J,J)
      DO 8 J=1,N
      DO 9 I=1,N
  9   U(I,J)=0.0
  8   U(J,J)=1.0
 10   AMAX=0.0
      DO 11 I=2,N
      JUP=I-1
      DO 11 J=1,JUP
      AII=P(I)
      AJJ=P(J)
      EPS=EP*EP
      AOD=A(I,J)
      ASQ=AOD*AOD
 28   IF(ASQ-AMAX)23,23,27
 27   AMAX=ASQ
 23   IF(ASQ-EPS)11,11,12
 12   DIF=AII-AJJ
      IF(DIF)13,15,15
 13   SIGN=-2.0
      DIF=-DIF
      GO TO 16
 15   SIGN=2.0
 16   TDEN=DIF+SQRT(DIF*DIF+4.0*ASQ)
      TANK=SIGN*AOD/TDEN
      C=1.0/(SQRT(1.0+TANK*TANK))
      S=C*TANK
      DO 24 K=1,N
      XJ=C*U(K,J)-S*U(K,I)
      U(K,I)=S*U(K,J)+C*U(K,I)
      U(K,J)=XJ
      IF(K-J) 17,24,18
 17   XJ=C*A(J,K)-S*A(I,K)
      A(I,K)=S*A(J,K)+C*A(I,K)
      A(J,K)=XJ
      GO TO 24
 18   IF(K-I)19,24,21
 19   XJ=C*A(K,J)-S*A(I,K)
      A(I,K)=S*A(K,J)+C*A(I,K)
      A(K,J)=XJ
      GO TO 24
 21   XJ=C*A(K,J)-S*A(K,I)
      A(K,I)=S*A(K,J)+C*A(K,I)
      A(K,J)=XJ
 24   CONTINUE
      P(I)=C*C*AII+S*S*AJJ+2.0*S*C*AOD
      P(J)=C*C*AJJ+S*S*AII-2.0*S*C*AOD
      A(I,J)=0.0
 11   CONTINUE
      IF(AMAX-EPS)20,20,10
 20   DO 40 K=1,N
      ATEST=P(K)
      JTEST=K
      DO 41 J=K,N
      IF(ABS(P(J))-ABS(ATEST))42,41,41
 42   ATEST=P(J)
      JTEST=J
 41   CONTINUE
      P(JTEST)=P(K)
      P(K)=ATEST
      DO 40 I=1,N
      UTEST(I)=U(I,JTEST)
      U(I,JTEST)=U(I,K)
 40   U(I,K)=UTEST(I)
      RETURN
      END