PROGRAM DIAGOT
C  TEST OF DIAGONALISATION OF EIGENVALUES AND
C  EIGENVECTORS OF SMALL SYMMETRIC REAL MATRICES
C  UP TO DIMENSION 20
      COMMON/VN01B2/A(20,20),IDIM,S(20,20)
      READ(1,1)IDIM
 1    FORMAT(I2)
      WRITE(3,2)IDIM
 2    FORMAT(1H ,' DIMENSION OF MATRIX=',I4/)
      DO 5 I=1,IDIM
      READ(1,3)(A(I,J),J=1,I)
 3    FORMAT(10F8.5/10F8.5)
      WRITE(3,4)(A(I,J),J=1,I)
 4    FORMAT(1H ,'MATRIX',10F8.5/10F8.5)
 5    CONTINUE
      DO 10 I=1,IDIM
      DO 10 J=1,I
      A(J,I)=A(I,J)
 10   CONTINUE
      CALL VN0108(IDIM)
      DO 20 I=1,IDIM
      WRITE(3,6)A(I,I)
      WRITE(3,7)(S(J,I),J=1,IDIM)
 7    FORMAT(1H ,'EIGENVECTOR'/10F9.5/10F9.5)
 6    FORMAT(1H ,'EIGENVALUE',F9.5)
 20   CONTINUE
      STOP
      END
      SUBROUTINE VN0108(N)
C  VN0108 DIAGONALISES SQUARE MATRICES UP TO THE ORDER OF 20*20
C
      COMMON/VN01B2/A(20,20),IDIM,S(20,20)
       IF(N-1)21,21,22
22    DO 1 I=1,N
      S(I,I)=1.
      I1=I+1
       IF(I1.GT.N)GOTO 60
       DO 1 J=I1,N
      S(I,J)=0.
1     S(J,I)=S(I,J)
60    NUL=0
      NBMAX=1000
      EPS=1.0E-3
      IN=3
      EF=10.
       DO 2 I=1,IN
      EPS=EPS/EF
       DO 2 NB=NUL,NBMAX
      IF(NB)4,5,4
4     IF(IDR+IDI)2,2,5
5      IDR=0
      IDI=0
      DO 33 II=2,N
       JI=II-1
      DO 33 JJ=1,JI
      C=A(II,JJ)+A(JJ,II)
      D=A(II,II)-A(JJ,JJ)
      IF(ABS(C)-EPS)6,7,7
6     CC=1.
      SS=0.
      GOTO 8
7      CC=D/C
       IF(CC)10,9,9
9      SIG=1.
       GOTO 12
10    SIG=-1.
12    COT=CC+SIG*SQRT(1.+CC*CC)
      SS=SIG/SQRT(1.+COT*COT)
      CC=SS*COT
      IDR=IDR+1
8     E=A(II,JJ)-A(JJ,II)
      CH=1.
      SH=0.
      IF(ABS(E)-EPS)14,13,13
13     CO=CC*CC-SS*SS
       SI=2.*SS*CC
      H=0.
       G=0.
       HJ=0.
      DO 15 K=1,N
      IF(K-II)16,15,16
16     IF(K-JJ)18,15,18
18     H=H+A(II,K)*A(JJ,K)-A(K,II)*A(K,JJ)
      S1=A(II,K)*A(II,K)+A(K,JJ)*A(K,JJ)
       S2=A(JJ,K)*A(JJ,K)+A(K,II)*A(K,II)
       G=G+S1+S2
      HJ=HJ+S1-S2
15    CONTINUE
       D=D*CO+C*SI
      H=2.*H*CO-HJ*SI
      F=(2.*E*D-H)/(4.*(E*E+D*D)+2.*G)
       IF(ABS(F)-EPS)14,17,17
17     CH=1./SQRT(1.-F*F)
       SH=F*CH
      IDI=IDI+1
14    C1=CH*CC-SH*SS
       C2=CH*CC+SH*SS
       S1=CH*SS+SH*CC
       S2=SH*CC-CH*SS
       IF(ABS(S1)+ABS(S2))20,33,20
20    DO 30 L=1,N
      A1=A(L,II)
       A2=A(L,JJ)
       A(L,II)=C2*A1-S2*A2
      A(L,JJ)=C1*A2-S1*A1
      A1=S(L,II)
       A2=S(L,JJ)
       S(L,II)=C2*A1-S2*A2
30    S(L,JJ)=C1*A2-S1*A1
      DO 31 L=1,N
      A1=A(II,L)
       A2=A(JJ,L)
       A(II,L)=C1*A1+S1*A2
31    A(JJ,L)=C2*A2+S2*A1
33    CONTINUE
2      CONTINUE
      TAUSQ=0.
      N1=N-1
      DO 150 I=1,N1
      II=I+1
      DO 150 J=II,N
150   TAUSQ=TAUSQ+(A(I,J)+A(J,I))**2
       RETURN
21     S(1,1)=1.
      RETURN
       END
C
C  TEST OF DIAGONALIZATION FOR A 10X10 MATRIX
C  READING OF A COMMAND FILE TO TEST
C
$ ASSIGN SYS$INPUT FOR001
$ ASSIGN SYS$OUTPUT FOR003
R DIAGOT.EXE
10
 3.00000
 0.05000 2.80000
 0.00000 0.05000 2.90000
 0.00000 0.00000 0.05000 2.85000
 0.00000 0.00000 0.00000 0.05000 2.75000
 0.00000 0.00000 0.00000 0.00000 0.05000 2.70000
 0.00000 0.00000 0.00000 0.00000 0.00000 0.05000 2.60000
 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.05000 2.55000
 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.05000 2.50000
 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.05000
 2.00000