C JACOBI11.F OF N.TAJIMA'S FORTRAN BENCH-MARK TESTS (VER.2)
C REPEATS A PROCEDURE ALMOST LIKE JACOBI DIAGONALIZATION FOR DIM 11.
C HISOTRY: 97/6/26
C
      PROGRAM JACOB1
      IMPLICIT REAL*8 (A-H,O-Z)
      CALL JACOB2(1.0D0, 0.7D0, 5, 200000)
      END
C
      SUBROUTINE JACOB2(PARM1,PARM2,NCYCLE,LMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(N=11,M1=1289,M2=104701,M3=1048576)
      DIMENSION H0(N,N),H(N,N),V(N,N),E(N)
C
      L=1
      DO 10 I=1,N
      DO 10 J=1,I
        IF(I.EQ.J) THEN
          H0(I,J)=I-1
        ELSE
          L=MOD(M1*L+M2,M3)
          H0(I,J)=L*PARM1/M3+(I+J-2*N)*PARM1/(2*N)
          H0(J,I)=H0(I,J)
        END IF
   10 CONTINUE
C
      DO 20 I=1,N
      DO 20 J=1,N
   20 H(I,J)=H0(I,J)
C
      FCT=(N/6.0D0)**2
      DO 50 L=1,LMAX-1
        CALL JACOBI(H,N,N,NCYCLE, E,V)
        CK=DBLE(N*L)/LMAX
        DO 40 I=1,N
        DO 40 J=1,I
          S=0
          DO 30 K=1,N
   30     S=S+V(K,I)*V(K,J)*FCT/((K-CK)**2+FCT)
          H(I,J)=H0(I,J)+PARM2*S*ABS(E(I)-E(J))
          IF(I.NE.J) H(J,I)=H(I,J)
   40   CONTINUE
   50 CONTINUE
      CALL JACOBI(H,N,N,NCYCLE, E,V)
C
      WRITE(6,950) N,NCYCLE,LMAX,E(1),E(N)
  950 FORMAT(' DIMENSION=',I3,'  #JACOBI ROT CYC=',I2,'*',I7
     & ,'  EMIN=',F11.6,'  EMAX=',F11.6)
C
      END
C
C----- JACOBI --------------------------------------------------------
      SUBROUTINE JACOBI(A,N,N1,ITERMX, E,V)
C
C THIS IS A SIMPLIFIED VERSION OF SUBR. JACOBI FOR BENCH MARK TEST
C
C JACOBI DIAGONALIZATION (CYCLIC)
C   SORTING EIGENVALUES IN THE ASCENDING ORDER
C
C *** HISTORY ***
C 3/APR/97 : SIMPLIFIED FOR A BENCH MARK TEST
C 4/APR/95 : JACOBI1.FOR WAS CREATED FROM SCRATCH
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(MXN=11,EPS1=1.0D-20)
      DIMENSION A(N1,N),E(N),V(N1,N),IDX(MXN),W(MXN)
C
      IF(N.GT.MXN) STOP 1
      DO 20 J=1,N
      DO 10 I=1,N
   10 V(I,J)=0
   20 V(J,J)=1
      ITER=0
   30 ITER=ITER+1
      IF(ITER.GT.ITERMX) GO TO 90
C TO MAKE #OPERATIONS UNCHANGED, FOLLOWING 7 LINES ARE INSERTED
      DO 40 I=2,N
      DO 40 J=1,I-1
   40 A(I,J)=A(I,J)-N*EPS1/((I-J)*ITER)
      DO 87 IP=1,N-1
      DO 86 IQ=IP+1,N
      IF(ABS(A(IQ,IP)).LT.EPS1)
     & A(IQ,IP)=A(IQ,IP)+(A(IQ,IQ)-A(IP,IP))/1000+EPS1*10
C
      THETA=(A(IQ,IQ)-A(IP,IP))/(2*A(IQ,IP))
      X=ABS(THETA)
      IF(X.LT.1.0D+8) THEN
        T=SIGN(1/(X+SQRT(X**2+1)),THETA)
      ELSE
        T=1/(2*THETA)
      END IF
      C2=1/(T**2+1)
      C=SQRT(C2)
      S=T*C
      S2=S**2
      TAU=S/(1+C)
C
      DO 50 J=1,IP-1
        TMP    =A(IP,J)-S*(A(IQ,J)+TAU*A(IP,J))
        A(IQ,J)=A(IQ,J)+S*(A(IP,J)-TAU*A(IQ,J))
        A(IP,J)=TMP
   50 CONTINUE
C
      DO 60 K=IP+1,IQ-1
        TMP    =A(K,IP)-S*(A(IQ,K)+TAU*A(K,IP))
        A(IQ,K)=A(IQ,K)+S*(A(K,IP)-TAU*A(IQ,K))
        A(K,IP)=TMP
   60 CONTINUE
C
      DO 70 I=IQ+1,N
        TMP    =A(I,IP)-S*(A(I,IQ)+TAU*A(I,IP))
        A(I,IQ)=A(I,IQ)+S*(A(I,IP)-TAU*A(I,IQ))
        A(I,IP)=TMP
   70 CONTINUE
C
      G       =2*A(IQ,IP)*S*C
      TMP     =C2*A(IP,IP)+S2*A(IQ,IQ) -G
      A(IQ,IQ)=S2*A(IP,IP)+C2*A(IQ,IQ) +G
      A(IP,IP)=TMP
      A(IQ,IP)=0
C
      DO 80 I=1,N
        TMP    =C*V(I,IP)-S*V(I,IQ)
        V(I,IQ)=S*V(I,IP)+C*V(I,IQ)
        V(I,IP)=TMP
   80 CONTINUE
C
   85 CONTINUE
   86 CONTINUE
   87 CONTINUE
C
      GO TO 30
   90 CONTINUE
C
      DO 100 I=1,N
  100 E(I)=A(I,I)
C
      DO 110 I=1,N
  110 IDX(I)=I
      DO 130 I=1,N-1
        I0=I
        E0=E(I)
        DO 120 J=I+1,N
          IF(E(J).LT.E0) THEN
            I0=J
            E0=E(J)
          END IF
  120   CONTINUE
        IF(I0.NE.I) THEN
          E(I0)=E(I)
          E(I)=E0
          IT     =IDX(I0)
          IDX(I0)=IDX(I)
          IDX(I) =IT
        END IF
  130 CONTINUE
C
      DO 160 I=1,N
        DO 140 J=1,N
  140   W(J)=V(I,IDX(J))
        DO 150 J=1,N
  150   V(I,J)=W(J)
  160 CONTINUE
C
C     WRITE(6,900) E(1),E(2),E(3),E(N-1),E(N)
C 900 FORMAT(1P5E15.7)
C
      END