C LEQS1K.F OF N.TAJIMA'S FORTRAN BENCH-MARK TESTS (VER.2)
C SOLVES LINEAR EQUATIONS FOR DIMENSION 1001.
C HISTORY: 97/6/26
C
      PROGRAM LEQS1
      IMPLICIT REAL*8 (A-H,O-Z)
      CALL LEQS2(1.0D0, 1.0D0, 1.0D0, 5)
      END
C
      SUBROUTINE LEQS2(PARM0,PARM1,PARM2,ITERMX)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(N=1001,M1=1289,M2=104701,M3=1048576)
      COMMON /CMNLEQ/ H(N,N),X(N),P(0:3)
C
      P(0)=PARM0
      P(1)=PARM1
      P(2)=PARM2
      L=1
      DO 10 I=1,N
   10 X(I)=(1.5D0*I-0.75D0)/N+0.5D0
      DO 60 ITER=1,ITERMX
        FCT= P(1)/(N*M3)
        SFT=-P(2)/(N*2)
        DO 20 I=1,N
        DO 20 J=1,N
          IF(I.EQ.J) THEN
            H(I,J)=1+P(0)/I
          ELSE
            L=MOD(M1*L+M2,M3)
            H(I,J)=FCT*L+SFT
          END IF
   20   CONTINUE
        CALL SLEQS(H,N,N,X)
        DO 30 K=0,3
   30   P(K)=0
        DO 40 I=1,N-3,4
        DO 40 K=0,3
   40   P(K)=P(K)+X(I+K)**2
        S=4/(P(0)+P(1)+P(2)+P(3))
        DO 45 K=0,3
   45   P(K)=P(K)*S
        DO 50 I=1,N
   50   X(I)=SIGN(2-1.5D0/(ABS(X(I))+1),X(I))
C       WRITE(6,950) N,ITER,P(0),P(1),P(2)
   60 CONTINUE
      WRITE(6,950) N,ITERMX,P(0),P(1),P(2)
  950 FORMAT(' DIM=',I5,'  #INV=',I3,'  RESULTS:',3F14.10)
C
      END
C===== SLEQS.FOR ====================================================*
C 97/4/7 SIMPLIFIED VERSION FOR A BENCH MARK TEST
C A SUBROUTINE TO SOLVE LINEAR EQUATIONS: 25/MAY/94
C THE SOLUTION OF A*X=B IS RETURNED IN B. A IS DESTROYED.
      SUBROUTINE SLEQS(A,N,NA,B)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(MXN=1001)
      DIMENSION A(NA,N), B(N), LI(MXN)
      IF((N.LE.0).OR.(N.GT.NA).OR.(N.GT.MXN)) STOP 901
      DO 10 I=1,N
   10 LI(I)=I
      DO 60 K=1,N
        I=K
        D=ABS(A(LI(I),K))
        DO 20 I1=K+1,N
          D1=ABS(A(LI(I1),K))
          IF(D1.GT.D) THEN
            I=I1
            D=D1
          END IF
   20   CONTINUE
        IF(D.LT.1.0D-8) THEN
          WRITE(6,*) 'NOT GOOD FOR BENCH MARK TEST ', K,D
          STOP 902
        END IF
        IF(I.NE.K) THEN
          L    =LI(K)
          LI(K)=LI(I)
          LI(I)=L
        END IF
        L=LI(K)
        F=1.0D0/A(L,K)
        DO 30 J=K+1,N
   30   A(L,J)=A(L,J)*F
        B(L  )=B(L  )*F
        DO 50 I1=K+1,N
          L1=LI(I1)
          F=-A(L1,K)
          DO 40 J=K+1,N
   40     A(L1,J)=A(L1,J)+A(L,J)*F
          B(L1  )=B(L1  )+B(L  )*F
   50   CONTINUE
   60 CONTINUE
      DO 80 I=N,1,-1
        L=LI(I)
        S=B(L)
        DO 70 J=I+1,N
   70   S=S-A(L,J)*A(J,1)
        A(I,1)=S
   80 CONTINUE
      DO 90 I=1,N
   90 B(I)=A(I,1)
      END