      PROGRAM BENCH
C
C     THIS IS A BENCHMARK TEST PROGRAM TO TEST THE SPEED OF A
C     TYPICAL LINEAR EQUATION SOLUTION PROGRAM COMPILER UNDER
C     VARIOUS OPTIONS.
C
C     TEST PROGRAM BY HARRY M. MURPHY,  24 JANUARY 1980.
C     UPDATED BY H.M.M. ON 7 JANUARY 1981.
C
      INTEGER I,J,K,MXSIZE,NEQNS
C
      LOGICAL ERROR
C
      REAL A(70,70),X(70),Y(70),S(70)
      REAL BIGABS,BIGREL,SEED,SUM,SUMSQR,T1,T2,T3,TEVL,TSOL
C
      DATA MXSIZE/70/
      DATA NEQNS/70/
      DATA SEED/4H1981/
C
C                   START.  INITIALIZE RANDOM NUMBER GENERATOR.
    1 CALL SETRAN(SEED)
C                   EXERCISE RANDOM GENERATOR.
      DO 2 I=1,64
      V=RANF(0)
    2 CONTINUE
C                   OPEN THE PRINT FILE.
      OPEN(UNIT=6,NAME='LP:',TYPE='NEW',CARRIAGECONTROL='FORTRAN')
C                   REPEAT THIS TEST TEN TIMES
      DO 9 NTEST=1,10
C                   MARK START TIME.
      CALL MRKTIM
C                   GENERATE SOLUTION VECTOR, S.
      DO 3 I=1,NEQNS
      S(I)=AINT(200.0*(RANF(0)-0.5))
    3 CONTINUE
C                   GENERATE MATRIX, A.
      DO 5 J=1,NEQNS
      DO 4 I=1,NEQNS
      A(I,J)=0.0
      IF (RANF(1).GT.0.1) A(I,J)=AINT(200.0*(RANF(0)-0.5))
    4 CONTINUE
    5 CONTINUE
C                   CALCULATE GIVEN VECTOR, Y.
      DO 7 I=1,NEQNS
      SUM=0.0
      DO 6 J=1,NEQNS
      SUM=SUM+S(J)*A(I,J)
    6 CONTINUE
      Y(I)=SUM
    7 CONTINUE
C                   GET END OF SETUP TIME.  (START OF SOLUTION.)
      CALL RUNTIM(T1)
C                   CALL LINSYS TO SOLVE SYSTEM OF EQUATIONS.
      CALL LINSYS(A,X,Y,NEQNS,MXSIZE,ERROR)
      CALL RUNTIM(T2)
      IF (ERROR) STOP
C                   VALIDATE SOLUTION BY COMPARING X WITH S.
C                   GET BIGGEST ABSOLUTE AND RELATIVE ERRORS.
      BIGABS=0.0
      BIGREL=0.0
      SUMSQR=0.0
      DO 8 I=1,NEQNS
      ERR=X(I)-S(I)
      BIGABS=AMAX1(BIGABS,ABS(ERR))
      IF (S(I).NE.0.0) BIGREL=AMAX1(BIGREL,ABS(ERR/S(I)))
      SUMSQR=SUMSQR+ERR**2
    8 CONTINUE
C                   GET END OF VALIDATION TIME.
      CALL RUNTIM(T3)
C                   CALCULATE AND PRINT ELAPSED TIMES AND RESULTS.
      TSOL=T2-T1
      TEVL=T3-T2
C
      WRITE(6,10) NTEST,SEED,NEQNS,T1,TSOL,TEVL,T3,BIGABS,BIGREL,
     .            SUMSQR
C                   LOOP TEN TIMES
    9 CONTINUE
C                   CLOSE PRINT UNIT.
      CLOSE(UNIT=6)
C                   AND STOP.
      STOP
C
   10 FORMAT (1H1/'   Test number ',I2,' for seed "',A4,'".'/1X/
     . '   Number of equations =',I3,'.'/1X/
     . '   Eqn setup time  =',F6.2,' seconds.'/
     . '   Solution  time  =',F6.2,' seconds.'/
     . '   Evaluation time =',F6.2,' seconds.'/
     . '   Total run time  =',F6.2,' seconds.'/1X/
     . '   Max abs  error  =',1PE11.4/
     . '   Max rel  error  =',1PE11.4/
     . '   Sum sqr  error  =',1PE11.4/1X)
      END
      SUBROUTINE LINSYS (A,X,Y,N,M,ERROR)
C
C     SINGLE-PRECISION LINEAR EQUATION SOLVER.
C
C     SOLVES A SYSTEM OF N BY N LINEAR EQUATIONS OF THE FORM --
C
C                         A X = Y
C
C     THE ALGORITHM IS GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING.
C
C     NOTES --
C
C     (1)  THIS ROUTINE DESTROYS THE CONTENTS OF THE MATRIX, A, AND THE
C          CONTENTS OF THE GIVEN VECTOR, Y.
C     (2)  THE SOLUTION IS RETURNED IN THE UNKNOWN VECTOR, X.
C     (3)  ERROR MESSAGES ARE WRITTEN TO UNIT, ERR.
C     (4)  IF AN ERROR OCCURS, ERROR IS RETURNED AS .TRUE., OTHERWISE
C          ERROR IS RETURNED AS .FALSE.
C     (5)  ON THE MINC-11, RUN TIME FOR LARGE ARRAYS IS APPROXIMATELY
C          GIVEN BY    T = 1.52E-4 * N**3  SECONDS.
C
C     ROUTINE BY HARRY M. MURPHY,  AFWL/DYVA,  28 SEPTEMBER 1979.
C
      DIMENSION A(M,M), X(M), Y(M)
C
      INTEGER ERR, I, J, P, Q
C
      LOGICAL ERROR
C
      REAL BIG, T
C                   DEFINE ERROR RETURN UNIT.
      DATA ERR/7/
C
C                   START.  VALIDATE N AND M.
    1 ERROR=N.LT.1
      IF (ERROR) WRITE (ERR,12) N
      IF (ERROR) GO TO 11
C
      ERROR=N.GT.M
      IF (ERROR) WRITE (ERR,13) N,M
      IF (ERROR) GO TO 11
C                   REDUCE MATRIX, A, TO UPPER DIAGONAL.
      DO 7 P=1,N
C                   SEARCH FOR BIGGEST ABSOLUTE TERM IN PIVOT COLUMN.
      BIG=0.0
      DO 2 I=P,N
      T=ABS(A(I,P))
      IF (T.LE.BIG) GO TO 2
      BIG=T
      Q=I
    2 CONTINUE
C                   IF ONLY ZEROS IN PIVOT COLUMN ITS A SINGULAR SYSTEM.
      ERROR=BIG.EQ.0.0
      IF (ERROR) WRITE (ERR,14)
      IF (ERROR) GO TO 11
C                   IF BIG ALREADY IN PIVOT ROW, SKIP EXCHANGE.
      IF (Q.EQ.P) GO TO 4
C                   OTHERWISE, EXCHANGE ROWS.
      DO 3 J=P,N
      T=A(P,J)
      A(P,J)=A(Q,J)
      A(Q,J)=T
    3 CONTINUE
C                   EXCHANGE CORRESPONDING ITEMS IN GIVEN VECTOR, Y.
      T=Y(P)
      Y(P)=Y(Q)
      Y(Q)=T
C                   IF THIS IS THE LAST ROW, SKIP TO BACK-SOLUTION.
    4 IF (P.EQ.N) GO TO 8
C                   OTHERWISE, FORCE ZEROS BELOW THE PIVOT COLUMN.
      Q=P+1
      DO 6 I=Q,N
      T=A(I,P)/A(P,P)
      DO 5 J=Q,N
      A(I,J)=A(I,J)-T*A(P,J)
    5 CONTINUE
      Y(I)=Y(I)-T*Y(P)
    6 CONTINUE
C
    7 CONTINUE
C                   END OF UPPER DIAGONAL REDUCTION.
C
C                   BACK-SOLVE RESULTING SYSTEM FOR X.
    8 I=N
      X(I)=Y(I)/A(I,I)
      I=I-1
C                   LOOP WHILE I GREATER THAN ZERO.
    9 IF (I.EQ.0) GO TO 11
C
      T=Y(I)
      Q=I+1
      DO 10 J=Q,N
      T=T-X(J)*A(I,J)
   10 CONTINUE
      X(I)=T/A(I,I)
C                   DECREMENT I AND LOOP.
      I=I-1
      GO TO 9
C                   RETURN.
   11 RETURN
C
   12 FORMAT (41H LINSYS ERROR. LESS THAN 1 EQUATION.  N =,I3,1H.)
   13 FORMAT (53H LINSYS ERROR. NUMBER OF EQUATIONS EXCEED ARRAY SIZE./
     + 22H NUMBER OF EQUATIONS =,I3,15H.  ARRAY SIZE =,I3,3H BY,I3,1H.)
   14 FORMAT (38H LINSYS ERROR. THE SYSTEM IS SINGULAR.)
      END
      SUBROUTINE MRKTIM
C
C     THIS ROUTINE CALLS GTIM TO MARK THE START OF A TIMING RUN.
C     ELAPSED TIME OF THE RUN IS GIVEN BY A CALL TO RUNTIM.
C
C     ROUTINE BY HARRY M. MURPHY,  7 JANUARY 1981.
C
      COMMON /MKTCOM/ITZ(2)
C
    1 CALL GTIM(ITZ)
      RETURN
      END
      SUBROUTINE RUNTIM(SEC)
C
C     THIS ROUTINE, USED WITH SUBROUTINE MRKTIM, RETURNS THE ELAPSED
C     TIME, IN SECONDS, SINCE THE LAST CALL TO MRKTIM.
C
C     ROUTINE BY HARRY M. MURPHY,  7 JANUARY 1981.
C
      COMMON /MKTCOM/ITZ(2)
C
      INTEGER ITX(2),HX,HZ,MX,MZ,SX,SZ,TX,TZ
C
      REAL DH,DM,DS,DT
C
    1 CALL GTIM(ITX)
      CALL CVTTIM(ITZ,HZ,MZ,SZ,TZ)
      CALL CVTTIM(ITX,HX,MX,SX,TX)
C
      DH=HX-HZ
      DM=MX-MZ
      DS=SX-SZ
      DT=TX-TZ
      SEC=(DH*60.0+DM)*60.0+DS+DT/60.0
      RETURN
      END
