c performs 
c  matrix multiplication test
c  double precision version
c
c reference: W.F. Eddy, "The DECstation 3100 - UNIX for Power Users"
c            Chance, V.3,N.2 (1990), pp.42-47
c
c auxiliary routine: TIMER - returns elapsed cpu time in tics of .01 secs
C
      PARAMETER (N=1024)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL A(n,n), B(n,n), C(n,n)
      REAL CPTIM, random
      COMMON /RAND/ IX, IY, IZ
      CALL TIMER(ITICS0)
      IX = 1333
      IY = 1777
      IZ = 3555
      DO I = 1, N
       DO J = 1, N
        A(J,I) = DBLE(RANDOM(JUNK))
       ENDDO
       DO J = 1, N
        B(I,J) = DBLE(RANDOM(JUNK))
       ENDDO
      ENDDO
      DO I = 1, N
       DO J = 1, N
        C(I,J) = 0.0
        DO K = 1, N
         C(I,J) = C(I,J) + A(I,K) * B(K,J)
        ENDDO
       ENDDO
      ENDDO
      CKSUM = 0.0D0
      DO I = 1, N
       DO J = 1, N
        CKSUM = CKSUM + C(I,J)
       ENDDO
      ENDDO
      CALL TIMER(ITICS1)
      CPTIM = FLOAT(ITICS1-ITICS0) / 100.0


c output results to file 
      OPEN (UNIT=6,FILE='matmul.out')
      WRITE (6,10000)
      WRITE (6,10005) CKSUM
      WRITE (6,10010) N
      WRITE (6,10015) CPTIM
      CLOSE (6)

      STOP
10000 FORMAT (/'  Matrix Multiplication Timing Test')
10005 FORMAT ('  check sum = ',1PD20.5)
10010 FORMAT ('  dimension of matrix = ',I4)
10015 FORMAT ('  cpu time (sec)      = ', F10.2 )
      END
c version: Feb 23 1992
      FUNCTION RANDOM(L)
c
c        ALGORITHM AS 183  APPL. STATIST. (1982) VOL.31, NO.2
c
c        RETURNS A PSEUDO-RANDOM NUMBER RECTANGULARLY DISTRIBUTED
c        BETWEEN 0 AND 1
c
c        IX, IY AND IZ SHOULD BE SET TO INTEGER VALUES BETWEEN
c        1 AND 30000 BEFORE FIRST ENTRY
c
c        INTEGER ARITHMETIC UP TO 30323 IS REQUIRED
c
      COMMON /RAND/ IX, IY, IZ
      IX = 171 * MOD(IX,177) - 2 * (IX/177)
      IY = 172 * MOD(IY,176) - 35 * (IY/176)
      IZ = 170 * MOD(IZ,178) - 63 * (IZ/178)
c
      IF (IX.LT.0) IX = IX + 30269
      IF (IY.LT.0) IY = IY + 30307
      IF (IZ.LT.0) IZ = IZ + 30323
c
c        IF INTEGER ARITHMETIC UP TO 5212632 IS AVAILABLE,
c        THE PRECEDING 6 STATEMENTS MAY BE REPLACED BY
c
c        IX = MOD(171 * IX, 30269)
c        IY = MOD(172 * IY, 30307)
c        IZ = MOD(170 * IZ, 30323)
c
c        ON SOME MACHINES, THIS MAY SLIGHTLY INCREASE
c        THE SPEED. THE RESULTS WILL BE IDENTICAL.
c
      RANDOM = AMOD(FLOAT(IX)/30269.0+FLOAT(IY)/30307.0+FLOAT(IZ)/
     +  30323.0,1.0)
c
      IF (RANDOM.GT.0.0) RETURN
      RANDOM = DMOD(DBLE(FLOAT(IX))/30269.0D0+DBLE(FLOAT(IY))/30307.0D0+
     +  DBLE(FLOAT(IZ))/30323.0D0,1.0D0)
      IF (RANDOM.GE.1.0) RANDOM = 0.999999
      RETURN
      END