        PROGRAM BENCH2
C
C    -- These Fortran kernels are taken from some well known benchmarks.
C
C    -- Note: Parts of this code that are (or may be) VMS-specific are preceded
C    -- by a comment line containing the string "-- VMS Specific --".
C
C
        COMMON /HEAD/ NAME, NKERN, TIMES, BFNAME, CPUNAM
        COMMON /LUNS/ IFILE, ITTYIN, ITTYOU, TFILE

        CHARACTER*80 BFNAME
        CHARACTER*6 NAME(50)
        CHARACTER*15 CPUNAM
        INTEGER TFILE
        REAL TIMES(50)
C
        DATA          IFILE, ITTYIN, ITTYOU, TFILE /4, 5, 6, 7/
C               
C      INITIALIZE THE KU VAX VERSION OF SECOND() BY CALLING INITTIME.
C
      CALL INITTIME
C
        NKERN = 0
        DO 5 I = 1, 50
          TIMES (I) = 0.0
    5   CONTINUE
C
        CALL INFO
   15   FORMAT (A)
C               
C--------------------------------------------------
C Fast fourier transform (real)
C
  10    NKERN = NKERN + 1
       NAME(NKERN) = 'FFTR  '           
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL FFTR
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - Hydro excerpt
C
  20    NKERN = NKERN + 1
       NAME(NKERN) = 'LAHYDR'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAHYDR
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - banded linear equations
C
  30    NKERN = NKERN + 1
       NAME(NKERN) = 'LABLE '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LABLE
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - tridiagonal elimination below diagonal
C
  40    NKERN = NKERN + 1
       NAME(NKERN) = 'LATRD1'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LATRD1
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - tridiagonal elimination above diagonal
C
  50    NKERN = NKERN + 1
       NAME(NKERN) = 'LATRD2'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LATRD2
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - pde integration
C
  60    NKERN = NKERN + 1
       NAME(NKERN) = 'LAPDE '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAPDE
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - difference predictors
C
  70    NKERN = NKERN + 1
       NAME(NKERN) = 'LADIFF'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LADIFF
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - first sums
C
  80    NKERN = NKERN + 1
       NAME(NKERN) = 'LAFSUM'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAFSUM
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - 2D particle pusher
C
  90    NKERN = NKERN + 1
       NAME(NKERN) = 'LA2DPP'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LA2DPP
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C L.L.N.L. - 1D particle pusher
C
  100   NKERN = NKERN + 1
       NAME(NKERN) = 'LA1DPP'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LA1DPP
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C New L.L.N.L. - 2D Hydro Excerpt
C
  110  NKERN = NKERN + 1
       NAME(NKERN) = 'LA2DHY'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LA2DHY
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C New L.L.N.L. - First minimum
C
  120  NKERN = NKERN + 1
       NAME(NKERN) = 'LAFIMI'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAFIMI
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C New L.L.N.L. - Implicit conditional computation
C
  130  NKERN = NKERN + 1
       NAME(NKERN) = 'LAIMCO'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAIMCO
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C New L.L.N.L. - Matrix product
C
  140  NKERN = NKERN + 1
       NAME(NKERN) = 'LAMAPR'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LAMAPR
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Kernels from linpack benchmark
C
  150  NKERN = NKERN + 1
       NAME(NKERN) = 'LINPAC'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LINPAC
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C AIRCO real nested loop
C
  160  NKERN = NKERN + 1
       NAME(NKERN) = 'AIRREL'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL AIRREL
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Double precision Whetstone
C
  180  NKERN = NKERN + 1
       NAME(NKERN) = 'DWHET '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL DWHET
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Double precision nested loops
C
  190  NKERN = NKERN + 1
       NAME(NKERN) = 'DOUBLE'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL DOUBLE
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Single precision matrix inversion
C
  200  NKERN = NKERN + 1
       NAME(NKERN) = 'EGYPT '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL EGYPT
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Dingle precision matrix inversion
C
  210  NKERN = NKERN + 1
       NAME(NKERN) = 'EGYPTD'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL EGYPTD
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Guass quadrature integration
C
  220  NKERN = NKERN + 1
       NAME(NKERN) = 'GAUSS '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL GAUSS
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Towers of hanoi integer problem
C
  230  NKERN = NKERN + 1
       NAME(NKERN) = 'HANOI '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL HANOI
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Real frenel lens calculation
C
  240  NKERN = NKERN + 1
       NAME(NKERN) = 'NEFF  '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL NEFF
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Double precision matrix inversion
C
  250  NKERN = NKERN + 1
       NAME(NKERN) = 'PHILCO'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL PHILCO
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Integer prime number calculation
C
  260  NKERN = NKERN + 1
       NAME(NKERN) = 'PRIME '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL PRIME
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Real steady state simulation
C
  270  NKERN = NKERN + 1
       NAME(NKERN) = 'RR2   '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL RR2
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Real rigs calibration
C
  280  NKERN = NKERN + 1
       NAME(NKERN) = 'RR3   '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL RR3
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Real blade cooling calculation
C
  290  NKERN = NKERN + 1
       NAME(NKERN) = 'RR4   '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL RR4
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Single precision nested loops
C
  300  NKERN = NKERN + 1
       NAME(NKERN) = 'SINGLE'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL SINGLE
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Single precision Whetstone
C
  310  NKERN = NKERN + 1
       NAME(NKERN) ='SWHET '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL SWHET
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Alamos18 -- vector calculations like r=a*b+c*d where a,b,c,d are vectors
C
  320  NKERN = NKERN + 1
       NAME(NKERN) ='ALAM18'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL ALAM18
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Gamma function; parallelizable
C
  330  NKERN = NKERN + 1
       NAME(NKERN) ='GAMSIM'
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL GAMSIM
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Double precision LU decomposition
C
  340  NKERN = NKERN + 1
       NAME(NKERN) ='LUDD  '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LUDD
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C Single precision LU decomposition
C
  350  NKERN = NKERN + 1
       NAME(NKERN) ='LUSD  '
       WRITE (*, 15) NAME(NKERN)
       CALL INIT(NAME(NKERN))
       CPU0 = SECOND()
       CALL LUSD
       TIMES (NKERN) = SECOND() - CPU0
C
C--------------------------------------------------
C
C -- this can be re-included to dump the array out before calling the analysis
C -- routine, in case the analysis routine isn't behaving.
C        TYPE *,'Contents of TIMES() array:'
C        DO 1990 I = 1, NKERN
C          TYPE *,TIMES (I)
C1990    CONTINUE

        CALL ANALYSIS

        CLOSE (UNIT=IFILE)
        END

C-----------------------------------------------------------------------
        SUBROUTINE INFO
C
C Get information from the user for the report header
C
        COMMON /HEAD/ NAME, NKERN, TIMES, BFNAME, CPUNAM
        COMMON /LUNS/ IFILE, ITTYIN, ITTYOU, TFILE

        CHARACTER*80 BFNAME
        CHARACTER*6 NAME(50)
        CHARACTER*15 CPUNAM
        INTEGER NKERN, TFILE
        REAL TIMES(50)
        CHARACTER*72 LINE, FILENA, TFILEN
C
        WRITE (ITTYOU, 10)
   10   FORMAT (1X, 'Benchmark Program', /, 
     1          1X, 'File name for timing report:')
C
        READ  (ITTYIN, 20) FILENA
   20   FORMAT (A)
C
        WRITE (ITTYOU, 25)
   25   FORMAT (1X, 'File name for file to be imported into 20/20:')
        READ  (ITTYIN, 20) TFILEN

        OPEN  (UNIT=IFILE, FILE=FILENA, STATUS='NEW',
     1         ERR=1000)
        OPEN  (UNIT=TFILE, FILE=TFILEN, STATUS='NEW',
     1         ERR=1010)
C
        WRITE (IFILE, 30) 
   30   FORMAT (1X, 'Digital Review Fortran Benchmark', /,
     1          1X, '--------------------------------', /)
C
        WRITE (ITTYOU, 40) 
   40   FORMAT (1X, 'Operating system (include version number):')
        READ  (ITTYIN, 20) LINE
        WRITE (IFILE, 50) LINE
   50   FORMAT (1X, 'System:    ', A)
C
        WRITE (ITTYOU, 60) 
   60   FORMAT (1X, 'Processor description (72 chars max):')
        READ  (ITTYIN, 20) LINE
        WRITE (IFILE, 70) LINE
   70   FORMAT (1X, 'Hardware:  ', A)
        WRITE (ITTYOU, 72)
   72   FORMAT (1X, 'Short Processor description (15 chars max):')
        READ  (ITTYIN, 20) CPUNAM
C
        WRITE (ITTYOU, 80) 
   80   FORMAT (1X, 'Your name:')
        READ  (ITTYIN, 20) LINE
        WRITE (IFILE, 90) LINE
   90   FORMAT (1X, 'Run by:    ', A)
C
        WRITE (ITTYOU, 92) 
   92   FORMAT (1X, 'Location of base times file:')
        READ  (ITTYIN, 20) BFNAME
C
        WRITE (ITTYOU, 100) 
  100   FORMAT (1X, 'Running benchmarks.  This will take a while...')
C
        RETURN
C
 1000   STOP 'Could not create report file.'
 1010   STOP 'Could not create file for 20/20 import.'
        END
        SUBROUTINE ANALYSIS
      
C Subroutine Analysis V2.0, from CPU.PAS V1.0
C 		Theoretical work by Jack F. Fegreus
C 		Pascal programming by Jeffrey A. Steinberg 4/11/88
C               Conversion to Fortran (ugh) by Denis W. Haskin 9/3/88
C 		DR Labs, Digital Review.  Copyright 1988 Ziff-Davis.
C 
C This analysis subroutine takes the CPU times of the individual benchmark
C kernels and performs various analyses on them, as well as comparing them
C to the DR Labs' base times for the MicroVAX II and VAX Fortran compiler.
C 
C The output of this program is defined by report_file and 
C contains the results for the MicroVAX II, the cpu in 
C question, the relative performance (MVII/CPU), the Pseudo
C MVUPS, and the scaled residuals for each of the modules.  At
C the bottom is overall geometric means, a Jacknife value  for 
C the test run and then a summary that contains the standard 
C deviation, the variance and the range of the jackknife based 
C on a 95% confidence interval.  To change the confidence 
C interval, the Case statement must be changed with the 
C correct values obtained from a statistical text book for the 
C T values.
C     Modification history:
C       
C       V2.1 January 23 1989 Denis Haskin Change geometric algorithm
C            to avoid large numbers which may cause overflow.  Allows routine
C            to be more portable.  Thanks to Bart Kessler at Sequent Computer
C            Systems.

      COMMON /HEAD/ NAME, NKERN, TIMES, BFNAME, CPUNAM
      COMMON /LUNS/ IFILE, ITTYIN, ITTYOU, TFILE

      CHARACTER*80 BFNAME
      CHARACTER*6 BNAME(50), NAME(50)
      CHARACTER*15 CPUNAM
      LOGICAL TOOFEW

      REAL TIMES(50)
      REAL BTIME(50), BTIMET(50), RPERF(50),
     1     PSOA, PBTIME, PRP,
     1     GSOA, GBTIME, GMVUP, PPS,
     1     GPS(50), PMVUP(50), JKN, JKNT, JKNT1, S2, SSTAR, VARIAN,
     1     DRSTAN, CONINT, LOWLIM, UPPLIM, SR(50)
      INTEGER BFILE, BKERNS, I, J, NKERN, TFILE
                              
      DATA IFILE, ITTYIN, ITTYOU, TFILE, BFILE /4, 5, 6, 7, 8/
      DATA PSOA/0.0/
      DATA PBTIME/0.0/, PRP/0.0/
      DATA PPS/0.0/, JKN/0.0/, JKNT/0.0/, JKNT1/0.0/
      DATA S2/0.0/
      DATA SSTAR/0.0/
      DATA VARIAN/0.0/, DRSTAN/0.0/, CONINT/0.0/
      DATA LOWLIM/0.0/
      DATA UPPLIM/0.0/
              
C open and read in a similar array with the names and times for our base
C CPU (MicroVAX II)

      OPEN(UNIT=BFILE, FILE=BFNAME, STATUS='OLD', FORM='FORMATTED')
      DO 100 BKERNS = 1,50
         READ (BFILE,*,END=110) BNAME(BKERNS), BTIMET(BKERNS)
 100  CONTINUE
 110  CLOSE(UNIT=BFILE)
      BKERNS = BKERNS - 1

      DO 140 I=1,NKERN
         DO 120 J=1,BKERNS
            IF (NAME(I).EQ.BNAME(J)) GO TO 130
 120     CONTINUE
C     -- If we get here we can't find our test in the base set; there's
C     -- some sort of problem.  Flag out and leave.
C     -- Eventually should have better error message.
         TYPE *,I,J,NAME(I),BNAME(J)
         STOP 'Test not found in base set.'
C     -- Match up times and accumulate into product for
C     -- later calculation of geometric mean.
 130     BTIME(I) = BTIMET(J)
         RPERF(I) = BTIME(I)/TIMES(I)
C     PSOA = PSOA * TIMES(I)
C     PBTIME = PBTIME * BTIME(I)
C     PRP = PRP * RPERF(I)
         PSOA = PSOA + LOG(TIMES(I))
         PBTIME = PBTIME + LOG(BTIME(I))
         PRP = PRP + LOG(RPERF(I))
 140  CONTINUE
      
C     -- Do geometric means for all this stuff.
C     -- VMS Specific --
C     GSOA = PSOA ** (1.0/QEXT(NKERN))
      GSOA = EXP((1.0/NKERN)*PSOA)
C     -- VMS Specific --
C     GBTIME = PBTIME ** (1.0/QEXT(NKERN))
      GBTIME = EXP((1.0/NKERN)*PBTIME)
C     -- VMS Specific --
C     GMVUP = PRP ** (1.0/QEXT(NKERN))
      GMVUP = EXP((1.0/NKERN)*PRP)
      
C     -- Figure pseudo mvups for each kernel
      DO 160 I = 1,NKERN
         PPS = 0.0
         GPS(I) = 0.0
         PMVUP(I) = 0.0
         DO 150 J = 1,NKERN
C     IF (J.NE.I) PPS = PPS * RPERF(J)
            IF (J.NE.I) PPS = PPS + LOG(RPERF(J))
 150     CONTINUE
C     GPS(I) = PPS ** (1.0/(NKERN-1.0))
         GPS(I) = EXP((1.0/(NKERN-1.0)) * PPS)
         PMVUP(I) = NKERN*GMVUP - (NKERN-1.0)*GPS(I)
C     -- VMS Specific --
C     JKN = JKN + QEXT(1.0/NKERN)*PMVUP(I)
         JKN = JKN + (1.0/NKERN)*PMVUP(I)
         JKNT = JKNT + PMVUP(I)**2.0
         JKNT1 = JKNT1 + PMVUP(I)
 160  CONTINUE
      
C     -- Do stat'l calculations for confidence intervals and variance.
C     -- VMS Specific --
C     S2 = ABS((JKNT - QEXT(1.0/NKERN)*JKNT1**2.0)/
C     1    QEXT(NKERN-1.0))
      S2 = ABS((JKNT - (1.0/NKERN)*JKNT1**2.0)/
     1     (NKERN-1.0))
C     -- VMS Specific --
C     SSTAR = SQRT(QEXT(S2/NKERN))
      SSTAR = SQRT((S2/NKERN))
      VARIAN = SSTAR**2.0      
      DRSTAN = GBTIME/GSOA
      
      TOOFEW = .FALSE.
      IF (NKERN.LT.20) THEN
 162     FORMAT(/'* Warning: Less than 20 kernels were run.')
 164     FORMAT( '* Variance, standard deviation, and 
     1        confidence intervals will'/'not be calculated.')
         WRITE (IFILE, 162)
         WRITE (IFILE, 164)
         TOOFEW = .TRUE.
      ELSE IF ((NKERN.GE.20).AND.(NKERN.LE.24)) THEN
         CONINT = 1.73 * SSTAR
      ELSE IF ((NKERN.GE.25).AND.(NKERN.LE.29)) THEN
         CONINT = 1.71 * SSTAR
      ELSE
         CONINT = 1.68 * SSTAR
      END IF                            
      
      IF (.NOT.TOOFEW) THEN
         LOWLIM = JKN - CONINT             
         UPPLIM = JKN + CONINT
         DO 170 I=1,NKERN                    
            SR(I) = (JKN - PMVUP(I))/SSTAR
 170     CONTINUE                 
      END IF
      
C     -- And now pump it all out the files.
C     -- First the text output...
 172  FORMAT(///T13,'Statistical Analysis of Benchmarks for ',A//)
 174  FORMAT('Kernel',T14,'MicroVAX II',T27,A,T42,A,T59,
     1     'Pseudo',T72,'Scaled')
 176  FORMAT('Name',T14,'Time (Sec.)',T27,'Time (Sec.)',T42,
     1     'MVUPs*',T59,'MVUPs*',T72,'Residuals'/)
      WRITE(IFILE, 172) CPUNAM
      WRITE(IFILE, 174) CPUNAM, CPUNAM
      WRITE(IFILE, 176)
      
C     -- And the 20/20 file.
 178  FORMAT('##'/'##'/'##'/'# # Statistical Analysis of'
     1     ,' Benchmarks for ',A,'##'/'##'/'##'/'##')
 180  FORMAT('Kernel#MicroVAX II#',A,'#',A,'#Pseudo#Scaled##')
 182  FORMAT('Name#Time (Sec.)#Time (Sec.)#MVUPs*#MVUPs*#Residuals##')
      WRITE(TFILE, 178) CPUNAM
      WRITE(TFILE, 180) CPUNAM, CPUNAM
      WRITE(TFILE, 182)
      
C     -- Do them together in this loop.
 184  FORMAT(A,T12,F10.4,T24,F10.4,T40,F10.4,T57,F10.4,T71,F10.4)
 186  FORMAT(A,'#',F10.4,'#',F10.4,'#',F10.4,'#',F10.4,'#',F10.4,'##')
      DO 200 I=1,NKERN
         WRITE (IFILE, 184) NAME(I), BTIME(I), TIMES(I), RPERF(I),
     1        PMVUP(I), SR(I)
         WRITE (TFILE, 186) NAME(I), BTIME(I), TIMES(I), RPERF(I),
     1        PMVUP(I), SR(I)
 200  CONTINUE
      
C     -- Finish off the text report file.
 210  FORMAT(/T13,'Geo. Mean',T25,'Geo. Mean',T41,'Geo. Mean',
     1     T53,'Jacknife Value')
 220  FORMAT(T12,F10.4,T24,F10.4,T40,F10.4,T57,F10.4)
 230  FORMAT(//'Summary of Statistical Analysis'/)
 240  FORMAT('DR Labs Standard MVUPs*:',F11.4/)
 250  FORMAT('Variance               :',F11.4)
 260  FORMAT('Standard Deviation     :',F11.4)
 270  FORMAT('95% Confidence Interval:',F11.4,' - ',F11.4/)
 280  FORMAT('* MVUPs = MicroVAX II Units of Processing.')
      WRITE (IFILE, 210)               
      WRITE (IFILE, 220) GBTIME, GSOA, GMVUP, JKN
      WRITE (IFILE, 230)
      WRITE (IFILE, 240) DRSTAN
      WRITE (IFILE, 250) VARIAN
      WRITE (IFILE, 260) SSTAR
      WRITE (IFILE, 270) LOWLIM, UPPLIM
      WRITE (IFILE, 280)
      
C     -- Finish off the 20/20 report file.
 300  FORMAT(/'#Geo. Mean#Geo. Mean#Geo. Mean#Jacknife Value##')
 310  FORMAT('#',F10.4,'#',F10.4,'#',F10.4,'#',F10.4,'##')
 320  FORMAT('##'/'Summary of Statistical Analysis##')
 330  FORMAT('##'/'DR Labs Standard MVUPs*:#',F11.4,'##'/'##')
 340  FORMAT('Variance               :#',F11.4,'##')
 350  FORMAT('Standard Deviation     :#',F11.4,'##')
 360  FORMAT('95% Confidence Interval:#',F11.4,'#',F11.4,'##')
 370  FORMAT('##'/'* MVUPs = MicroVAX II Units of Processing.##')
      WRITE (TFILE, 300)               
      WRITE (TFILE, 310) GBTIME, GSOA, GMVUP, JKN
      WRITE (TFILE, 320)
      WRITE (TFILE, 330) DRSTAN
      WRITE (TFILE, 340) VARIAN
      WRITE (TFILE, 350) SSTAR
      WRITE (TFILE, 360) LOWLIM, UPPLIM
      WRITE (TFILE, 370)
      
      CLOSE(UNIT=IFILE)
      CLOSE(UNIT=TFILE)
      
      RETURN
      END
        SUBROUTINE AIRREL
C
C...... AIRCO REAL BENCHMARK
C
        INTEGER I,J,L
        REAL    X,Y
        COMMON /AIRXY/ X, Y
        COMMON /LF/ LOOPFACTOR
C
        DO 30 L = 1,(10*LOOPFACTOR)
           DO 20 I = 1,43
              X = 1.0 + (I/500) - (I/750)
              DO 10 J = 1,1000
                 X = X + .001
                 Y = 3.0 * (X+1.1) / ((X+2.2)*X)
10            CONTINUE
20         CONTINUE
30      CONTINUE
C
        RETURN
        END
C     See bottom of this file for source information - DWH
C
C     *labmk18 (variation of labmk8-orig. vector calculations)
C     Renamed Alamos18 DWH
C
C      program lbmk18(out18,tape6=out18)
      SUBROUTINE ALAM18
C
C      this is los alamos bench mark program 18.the main program
C      is a controller for three subroutines performing vector
C      calculations, no gather scatter operations.  the code was last
C      modified december,1981.
C
C            maximum field length requirements for los alamos systems:
C               ctss               121540
C               ltss               114322
C               nos                107600
C
C     ****************************************************************
C
C       nloop is the number of different vector lengths to try
C       kloop is the number of times to try each case
C       ncase, n is the vector length
C                 initialize program
      COMMON /ALAMOS18/ NCASE(10), AA(5000), BB(5000), CC(5000),
     + DD(5000), RR(5000), JJ(5000)
      COMMON /LF/ LOOPFACTOR
      nloop = 7
      kloop = (420*LOOPFACTOR)
C      k = nloop + 9
      n = ncase(nloop)
C      total = 0
      s1 = .25
      s2 = 5.
C-----Following loop is done in INIT subroutine before this is called.
C      do 10 i=1,n     
C          aa(i) = i
C          bb(i) = aa(i)**.5
C          cc(i) = (i-1)**2
C          fi = i
C          dd(i) = 1.0/fi
C          jj(i) = i - k
C   10 continue
Cc
C       do vector times scaler plus vector times scaler
Cc
      do 100  i=1,kloop
              call vvspvs(rr, aa, s1, bb, s2, n)
  100 continue
Cc
C       do vector times vector plus vector times vector
Cc
      do 200 i=1,kloop
              call vvvpvv(rr, aa, bb, cc, dd, n)
  200 continue
Cc
C       do vector times vector plus vector
Cc
      do 300 i=1,kloop
              call vvvpv(rr, aa, bb, cc, n)
  300 continue
   50 continue
      return
Cc
C99999 format (////48h  following times for vector calculations of len,
C     * 3hgth, i10)
C99998 format (24h   elapsed time in one =, f8.3/20h   elapsed time in t,
C     * 4hwo =, f8.3/26h   elapsed time in three =, f8.3/8h   total,
C     * 17h time elapsed   =, f8.3)
C99997 format (28h  average time per element =, f15.8)
C99996 format (////39h          total time for all vectors = , f15.4)
      end

      subroutine vvspvs(rr, aa, s1, bb, s2, n)
      dimension aa(1), bb(1), rr(1)
Cc
C     calculate r=a*s1+b*s2 where a,b are vectors and s1,s2 are scalers
C     input
C     n = vector length
C     a, b, and r are vectors of length n
Cc
      do 10 i=1,n
          rr(i) = aa(i)*s1 + bb(i)*s2      
   10 continue
      return
      end
      subroutine vvvpvv(rr, aa, bb, cc, dd, n)
      dimension aa(1), bb(1), cc(1), dd(1), rr(1)
Cc
C     calculate r=a*b+c*d where a,b,c and d are vectors
C     input
C     n = vector length
C     a, b, c, and r are vectors of length n
Cc
      do 10 i=1,n
          rr(i) = aa(i)*bb(i) + cc(i)*dd(i)
   10 continue
      return
      end
      subroutine vvvpv(rr, aa, bb, cc, n)
      dimension aa(1), bb(1), cc(1), rr(1)
Cc
C     calculate r=a*b+c where a,b and c are vectors
C     input
C     n = vector length
C     a, b, c, and r are vectors of length n
Cc
      do 10 i=1,n
          rr(i) = aa(i)*bb(i) + cc(i)
   10 continue
      return
      end
C      subroutine tme(ncall)
C      if(ncall.eq.1) write(6,900)
C      call clock(i)
C      write(6,990) i
C      return
C  900 format(8h labmk18)
C  990 format(1a9)
C      end
Cc
C      subroutine jobtim(t)
C      call timing(t,t1,t2)
C      return
C      end
C
C Caveat receptor. <nbslib@cmr.icst.nbs.gov> or <nbslib@icst-cmr.arpa>,
C [John Antonishek] <ant@cmr.icst.nbs.gov> or <ant@icst-cmr.arpa>.
C -- with thanks to Eric Grosse at Bell Labs and Jack Dongarra at Argonne Labs
C ------------------------------------------------------------------------
C Compliments of nbslib   Mon Aug  8 16:49:51 EDT 1988
C -- nbslib software is in the public domain or else copyrighted by a
C    U. S. government contractor which licenses the government to use,
C    reproduce, and distribute it.  Neither the United States, nor the
C    Department of Commerce, nor the National Bureau of Standards, nor
C    any of its employees, makes any warranty, express or implied, or
C    assumes any legal liability or reponsibility for the accuracy,
C    completeness, or usefulness of any information, apparatus, product,
C    process, or software disclosed or represents that its use would not
C    infringe privately owned rights.
C ------------------------------------------------------------------------
        SUBROUTINE DOUBLE
C
C...... DOUBLE PRECISION FLOATING POINT BENCHMARK
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Y)
        INTEGER I,J,LOOP
        COMMON /DOUBL/ WINT(91)
        COMMON /LF/ LOOPFACTOR
C
        F(X) = ( ((X-2.)*X+3.) * X - 4. ) *X + 5.
C
        DO 40 LOOP = 1,(10*LOOPFACTOR)
          DO 10 I = 1,91
            WINT(I) = 0.
10        CONTINUE
          DO 30 J = 1,91,10
            X2 = 0.
            UP = 1.+(J-1)/10.
            STEP = UP/1000.
            DO 20 I = 1,1000
              X1 = X2
              X2 = I*STEP
              WINT(J) = WINT(J) + 
     1                  STEP*(F(X1)+F(X2)+4*F((X1+X2)/2.))/6.
20          CONTINUE
30        CONTINUE
40      CONTINUE
C
        RETURN
        END
        SUBROUTINE DWHET
C
C       WHETSTONE SYNTHETIC BENCHMARK PROGRAM, DOUBLE-PRECISION VARIANT
C
C       THIS PROGRAM IS THE BASIS OF THE WHETSTONE-INSTRUCTIONS-PER-SECOND
C       METRIC FOR FORTRAN AND CPU PERFORMANCE.
C
C       FOR MORE INFORMATION, SEE
C           "A SYNTHETIC BENCHMARK," H. J. CURNOW AND B. A. WICHMANN,
C           THE COMPUTER JOURNAL, VOLUME 19, NO. 1, PAGES 43-49.
C
C       NOTE:  EARLIER VERSIONS OF THIS CODING OF THE WHETSTONE BENCHMARK
C              COMPUTED N4, N8, AND N9 BELOW AS INTEGERS > 32767.
C              AS A CONSEQUENCE, THEIR VALUES ON 16-BIT PROCESSORS WERE
C              INCORRECTLY COMPUTED.  THE PROBLEM HAS BEEN FIXED IN
C              THIS VERSION BY PUTTING AN ADDITIONAL LOOP OF 10 AROUND
C              THE APPROPRIATE TESTS AND DECREASING N4, N8, AND N9 BY
C              10.  ALSO ALL VARIABLES HAVE BEEN EXPLICITLY TYPED.
C               
        DOUBLE PRECISION X1,X2,X3,X4,X,Y,Z,T,T1,T2,E1
        INTEGER I,II,J,K,L,N1,N2,N3,N4,N6,N7,N8,N9,N10,N11,LP
C
        COMMON /DWHETX/ T,T1,T2,E1(4),J,K,L
        COMMON /LF/ LOOPFACTOR
C
        I  = 100
        J  = I/10
        T  = 0.499975D00
        T1 = 0.50025D00
        T2 = 2.0D00
C
C
        N1  = 0
        N2  = I * 12
        N3  = I * 14
        N4  = J * 345
        N6  = I * 210
        N7  = I * 32
        N8  = J * 899
        N9  = J * 616
        N10 = 0
        N11 = I * 93
C
        DO 2000 LP=1,(1*LOOPFACTOR)
        X1 =  1.0D0
        X2 = -1.0D0
        X3 = -1.0D0
        X4 = -1.0D0
        IF ( N1 .LE. 0 ) GO TO 120
           DO 110,I = 1,N1,1
              X1 = ( X1+X2+X3-X4) * T
              X2 = ( X1+X2-X3+X4) * T
              X4 = (-X1+X2+X3+X4) * T
              X3 = ( X1-X2+X3+X4) * T
110           CONTINUE
120        CONTINUE
C        CALL POUTD (N1,N1,N1,X1,X2,X3,X4)
C
        E1(1) =  1.0D0
        E1(2) = -1.0D0
        E1(3) = -1.0D0
        E1(4) = -1.0D0
        IF ( N2 .LE. 0 ) GO TO 220
           DO 210 I = 1,N2,1
              E1(1) = (  E1(1)+E1(2)+E1(3)-E1(4) ) * T
              E1(2) = (  E1(1)+E1(2)-E1(3)+E1(4) ) * T
              E1(3) = (  E1(1)-E1(2)+E1(3)+E1(4) ) * T
              E1(4) = ( -E1(1)+E1(2)+E1(3)+E1(4) ) * T
210           CONTINUE
220        CONTINUE
C        CALL POUTD (N2,N3,N2,E1(1),E1(2),E1(3),E1(4))
C
        IF ( N3 .LE. 0 ) GO TO 320
           DO 310 I = 1,N3,1
              CALL PAD(E1)
310           CONTINUE
320        CONTINUE
C        CALL POUTD (N3,N2,N2,E1(1),E1(2),E1(3),E1(4))
C
        J = 1
        IF ( N4 .LE. 0 ) GO TO 490
           DO 480 II = 1,10
              DO 470 I = 1,N4,1
                 IF ( J .NE. 1 ) GO TO 410
                 J = 2
                 GO TO 420
410                 J = 3
420                 IF ( J .LE. 2 ) GO TO 430
                 J = 0
                 GO TO 440
430                 J = 1
440                 IF ( J .GE. 1 ) GO TO 450
                 J = 1
                 GO TO 460
450                 J = 0
460                 CONTINUE
470              CONTINUE
480           CONTINUE
490        CONTINUE
C        CALL POUTD (N4,J,J,X1,X2,X3,X4)
C
        J = 1
        K = 2
        L = 3
        IF ( N6 .LE. 0 ) GO TO 620
           DO 610 I = 1,N6,1
              J = J * (K-J) * (L-K)
              K = L * K - (L-J) * K
              L = (L-K) * (K+J)
              E1(L-1) = J + K + L
              E1(K-1) = J * K * L
610           CONTINUE
620        CONTINUE
C        CALL POUTD (N6,J,K,E1(1),E1(2),E1(3),E1(4))
C
        X = 0.5D0
        Y = 0.5D0
        IF ( N7 .LE. 0 ) GO TO 720
           DO 710 I = 1,N7,1
           X = T*DATAN(T2*DSIN(X)*DCOS(X)/(DCOS(X+Y)+DCOS(X-Y)-1.0D0))
           Y = T*DATAN(T2*DSIN(Y)*DCOS(Y)/(DCOS(X+Y)+DCOS(X-Y)-1.0D0))
710        CONTINUE
720        CONTINUE
C        CALL POUTD (N7,J,K,X,X,Y,Y)
C
        X = 1.0D0
        Y = 1.0D0
        Z = 1.0D0
        IF ( N8 .LE. 0 ) GO TO 830
           DO 820 II = 1,10
              DO 810 I = 1,N8,1
                 CALL P3D (X,Y,Z)
810              CONTINUE
820           CONTINUE
830        CONTINUE
C        CALL POUTD (N8,J,K,X,Y,Z,Z)
C
        J = 1
        K = 2
        L = 3
        E1(1) = 1.0D0
        E1(2) = 2.0D0
        E1(3) = 3.0D0
        IF ( N9 .LE. 0 ) GO TO 930
           DO 920 II = 1,10
              DO 910 I = 1,N9,1
                 CALL P0D
910              CONTINUE
920           CONTINUE
930        CONTINUE
C        CALL POUTD (N9,J,K,E1(1),E1(2),E1(3),E1(4))
C
        J = 2
        K = 3
        IF ( N10 .LE. 0 ) GO TO 1020
           DO 1010 I = 1,N10,1
              J = J + K
              K = J + K
              J = J - K
              K = K - J - J
1010           CONTINUE
1020        CONTINUE
C        CALL POUTD (N10,J,K,X1,X2,X3,X4)
C
        X = 0.75D0
        IF ( N11 .LE. 0 ) GO TO 1120
           DO 1110 I = 1,N11,1
              X = DSQRT ( DEXP( DLOG(X)/T1 ) )
1110           CONTINUE
1120        CONTINUE
C        CALL POUTD (N11,J,K,X,X,X,X)
C
C
C...... KILOWHETSTONE INSTRUCTIONS PER SECOND IS EQUAL TO 10000/TIME
2000    CONTINUE
        END
C
C***********************************************************************
C
        SUBROUTINE PAD(E)
C
        DOUBLE PRECISION T,T1,T2,E(4),E1
        COMMON /DWHETX/ T,T1,T2,E1(4),J,K,L
        INTEGER J,K,L
C
        J = 0
1        E(1) = (  E(1)+E(2)+E(3)-E(4) ) * T
        E(2) = (  E(1)+E(2)-E(3)+E(4) ) * T
        E(3) = (  E(1)-E(2)+E(3)+E(4) ) * T
        E(4) = ( -E(1)+E(2)+E(3)+E(4) ) / T2
        J = J + 1
        IF ( J .LT. 0 ) GO TO 1
C
        END
C
C***********************************************************************
C
        SUBROUTINE P0D
C
        DOUBLE PRECISION T,T1,T2,E1
        COMMON /DWHETX/ T,T1,T2,E1(4),J,K,L
        INTEGER J,K,L
C
        E1(J) = E1(K)
        E1(K) = E1(L)
        E1(L) = E1(J)
C
        END
C
C***********************************************************************
C
        SUBROUTINE P3D (X,Y,Z)
C
        DOUBLE PRECISION T,T1,T2,X1,X,Y1,Y,Z,E1
        COMMON /DWHETX/ T,T1,T2,E1(4),J,K,L
        INTEGER J,K,L
C
        X1 = X
        Y1 = Y
        X1 = T * (X1+Y1)
        Y1 = T * (X1+Y1)
        Z  = (X1+Y1) / T2
        END
C
C***********************************************************************
C
C        SUBROUTINE POUTD(N,J,K,X1,X2,X3,X4)
C
C        DOUBLE PRECISION X1,X2,X3,X4
C
C       WRITE STATEMENT COMMENTED OUT TO IMPROVE REPEATABILITY OF TIMINGS
C
C       WRITE(2,1)N,J,K,X1,X2,X3,X4
C1      FORMAT(1H,3I7,4E12.4)
C
C        END
        SUBROUTINE EGYPT
C
C        A SIMPLE MATRIX INVERSION BENCHMARK PROGRAM TO COMPARE THE COMPUTING
C        POWER OF DIFFERENT COMPUTERS.  TAKEN FROM THE GENERAL EGYPTIAN
C        ELECTRICITY CORP. RFP FOR A NATIONAL ENERGY CONTROL CENTER.
C       
        REAL S,Z
        COMMON /EGYP/ Z(50,50)
        COMMON /LF/ LOOPFACTOR
        INTEGER I,J,JA,K,L,LP
C
        L = 50
C
        DO 100 LP=1,(1*LOOPFACTOR)

        DO 20 I = 1,L
           S = FLOAT(I)
           DO 10 J = 1,L
              Z(I,J) = S
              S = S + 1
10           CONTINUE
           S      = Z(I,L)
           Z(I,L) = Z(I,I)
           Z(I,I) = S
20        CONTINUE
C
        DO 70 JA = 1,2
           DO 60 I = 1,L
              S = Z(I,I)
              Z(I,I) = 1./S
              DO 40 J = 1,L
                 IF (J .EQ. I) GO TO 40
                 DO 30 K = 1,L
                    IF (K .EQ. I) GO TO 30
                    Z(J,K) = Z(J,K) - Z(J,I)*Z(I,K)/S
30                 CONTINUE
40              CONTINUE
              DO 50 J = 1,L
                 IF (J .EQ. I) GO TO 50
                 Z(J,I) = Z(J,I)/S
                 Z(I,J) = -Z(I,J)/S
50              CONTINUE
60           CONTINUE
70        CONTINUE
C
100     CONTINUE
        RETURN
        END
        SUBROUTINE EGYPTD
C
C       A SIMPLE MATRIX INVERSION BENCHMARK PROGRAM TO COMPARE THE COMPUTING
C       POWER OF DIFFERENT COMPUTERS.  TAKEN FROM THE GENERAL EGYPTIAN
C       ELECTRICITY CORP. RFP FOR A NATIONAL ENERGY CONTROL CENTER.
C       
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        INTEGER I,J,JA,K,L,LP
        COMMON /LF/ LOOPFACTOR
        COMMON /EGYPD/ Z(50,50)
C
        L = 50
C
        DO 100 LP=1,(1*LOOPFACTOR)
        DO 20 I = 1,L
           S = FLOAT(I)
           DO 10 J = 1,L
              Z(I,J) = S
              S = S + 1
10         CONTINUE
           S      = Z(I,L)
           Z(I,L) = Z(I,I)
           Z(I,I) = S
20      CONTINUE
C
        DO 70 JA = 1,2
           DO 60 I = 1,L
              S = Z(I,I)
              Z(I,I) = 1./S
              DO 40 J = 1,L
                 IF (J .EQ. I) GO TO 40
                 DO 30 K = 1,L
                    IF (K .EQ. I) GO TO 30
                    Z(J,K) = Z(J,K) - Z(J,I)*Z(I,K)/S
30               CONTINUE
40            CONTINUE
              DO 50 J = 1,L
                 IF (J .EQ. I) GO TO 50
                 Z(J,I) = Z(J,I)/S
                 Z(I,J) = -Z(I,J)/S
50            CONTINUE
60         CONTINUE
70      CONTINUE
C
100     CONTINUE
        RETURN
        END
      SUBROUTINE FFTR                    
C
C A FAST-FOURIER TRANSFORM.
C
      REAL X(2048), U1, U2, W1, W2, T1, T3, T4
      COMMON /LF/ LOOPFACTOR
C
      DO 1000 ILOOP = 1, (20*LOOPFACTOR)
        M = 10
        N = 2**M
        TWON = 2*N
        PI = 3.14159265358979
        F=128.0
C
        DO 5 J=1,2048,2
          X(J) = SIN(3.14159*FLOAT(J)/F)
          X(J+1)=0.0
 5      CONTINUE
C
C COSINE GENERATION
C
        DO 20 L = 1,M
          LE = 2**(M+1-L)
          TWOLE = 2*LE
          LE1 = LE/2
          U1 = 1.0
          U2 = 0.0
          W1 = COS(PI/FLOAT(LE1))
          W2 = -SIN(PI/FLOAT(LE1))
C
C BUTTERFLY COMPUTATION
C
          DO 20 J = 1,LE,2
            DO 10 I = J,TWON,TWOLE
              IP = I + LE
              T3 = X(I) - X(IP)
              T4 = X(I+1) - X(IP+1)
              X(I) = X(I) + X(IP)
              X(I+1) = X(I+1) + X(IP+1)
              X(IP) = (T3*U1) - (T4*U2)
              X(IP+1) = (T3*U2) + (T4*U1)
 10         CONTINUE
            T1 = (U1*W1) - (U2*W2)
            U2 = (U1*W2) + (U2*W1)
            U1 = T1
 20     CONTINUE
C
C  BIT REVERSAL
C
        NM1 =2*(N-1)
        J = 1
        DO 30 I = 1, NM1, 2
          IF(I.GE.J) GOTO 25
          T1 = X(J)
          T2 = X(J+1)
          X(J) = X(I)
          X(J+1) = X(I+1)
          X(I) = T1
          X(I+1) = T2
 25       K = N
 26       IF(K.GE.J) GOTO 30
          J = J - K
          K = K/2
          GOTO 26
 30     J = J + K
C
1000  CONTINUE
      RETURN
      END


C       See bottom of this file for source information and subroutine
C	documentation.
C
       SUBROUTINE GAMSIM
C
C***********************************************************************
C       GLOBAL VARIABLES
       COMMON /LF/ LOOPFACTOR
       INTEGER I, MM, MS, N
       REAL    K, L, SUM
C       LOCAL VARIABLES
       REAL              ACCSUM

       DO 200 J=1,LOOPFACTOR
C       Use preset values.
         N = 4
         MM = 10
         K = 2
         MS = 3000
C                                           
         ACCSUM = 0
C+
C       Serial (uniprocessor) implementation of invoking the taret
C       routine: a DO-loop.  For parallel implementation, MM
C       concurrent processes would be started, each with its own
C       value of I and L.
C-
         DO 100 I = 1, MM
C+
C       Compute, as a function of an index I (identifying each of the
C       MM subintervals to be integrated by each invocation of the
C       target routine), the lower bound of the MM subintervals.  The
C       lower bound, L(I), is the function of I computed by K*(I-1),
C       I = 1, 2,...,MM.
C-
              L = K * REAL(I - 1)
              CALL SIMPSP (K, L, MS, N, SUM)
              ACCSUM = ACCSUM + SUM
 100     CONTINUE
 200   CONTINUE
        
       RETURN
       END
        
       SUBROUTINE SIMPSP (K, L, MS, N, SUM )
C              Last modified: 4/30/85                                                    Ver. 0.0
C****************************************************************************
C                                                                                                                                 *
C        NBS/ICST COMPUTER MEASUREMENT RESEARCH FACILITY                                      *
C                                                                                                                                 *
C        SIMPSP is a FORTRAN encoding of an algorithm that performs                     *
C        integration of a function over an interval by applying Simpson's         *
C        1/3 rule.  This subroutine is an ICST benchmark, one in a series of   *
C        computational tasks with known characteristics that are used to          *
C        measure and analyze the performance of parallel computers.                     *
C        See JORD83, GERA70.                                                                                       *
C        Ken Dymond 4/11/85                                                                                           *
C                                                                                                                                 *
C****************************************************************************
C        VAX VMS version SIMPSP.FOR
C        Language: ANSI FORTRAN X3.9-1978 (FORTRAN-77)
C
C       On entry:
C        K                 REAL
C                        The length of the interval of integration.  K > 0.
C
C        L              REAL
C                       Lower bound of the interval over which a fucntion is to
C              be integrated.
C
C       MS       INTEGER
C              The number of ordinates to be computed over the interval
C                      (L, L+K-(K/MS)).  MS is an even, positive integer.
C
C       N       INTEGER
C              A parameter required, not by Simpson's rule, but by the
C              function to be integrated, in this case the gamma
C              function for positive integer arguments.
C
C       Note: There is one other input which must be supplied at compile
C       time, and that is the function to be integrated.  It is inserted
C       as a FORTRAN statement function.
C
C       On return:
C
C       SUM       REAL
C                        The result of evaluating an integrand over an interval
C              by Simpson's 1/3 rule.
C
C       Subroutines called: None
C       Intrinsic functions invoked: EXP, REAL
C
C****************************************************************************
C       GLOBAL VARIABLES
       INTEGER        MS, N
       REAL               K, L, SUM
C
C       LOCAL VARIABLES
       INTEGER        I, MSHALF
       REAL               DELX, U, X
C+
C        STATEMENT FUNCTION
C       (GAMMA FUNCTION FOR POSITIVE INTEGER ARGUMENT N)
C-
         F(X,N) = (X**(N-1))*EXP(-X)
C
C       COMPUTE DELX, THE INCREMENT IN THE ABSCISSA AT WHICH
C       ORDINATES ARE TO BE EVALUATED.
C
       DELX = K/MS
C
C       COMPUTE UPPER BOUND, U
C
C       U = L + K - DELX
              U = L + K
C
C       COMPUTE 1ST, 2ND, AND LAST ORDINATE IN THE INTERVAL
C
       SUM = 0
       SUM = F(L,N) + 4.*F(L+DELX,N) + F(U,N)
C
C+
C       EVALUATE THE FUNCTION AT THE REMAINING MS/2 POINTS ON THE
C       ABSCISSA AND ADD THEM TO THE RUNNING SUM.
C-
C
C       COMPUTE ONE HALF OF MS: NUMBER OF LOOP ITERATIONS.
C
       MSHALF = MS/2
C
C       COMPUTE INITIAL ABSCISSA VALUE FOR LOOP, i.e., 3RD ORDINATE
C
       X = L + (2.*DELX)
C
C       MAJOR LOOP: THE NUMBER OF FLOATING POINT OPERATIONS IS
C       O(7 * Q), WHERE Q IS THE NUMBER OF FLOATING POINT
C       OPERATIONS EXECUTED WHENEVER THE TARGET FUNCTION
C       IS EVALUATED.
C
       DO 100 I = 2, MSHALF
              SUM = SUM + 2.*F(X,N) + 4.*F(X+DELX,N)
              X = X + 2.*DELX
100       CONTINUE
C
C       Apply the 1/3rd part of Simpson's rule
C
       SUM = DELX / 3. * SUM
C
C
       RETURN
       END
        
        
C          
C  Caveat receptor. <nbslib@cmr.icst.nbs.gov> or <nbslib@icst-cmr.arpa>,
C  [John Antonishek] <ant@cmr.icst.nbs.gov> or <ant@icst-cmr.arpa>.
C  -- with thanks to Eric Grosse at Bell Labs and Jack Dongarra at Argonne Labs
C  ------------------------------------------------------------------------
C  Compliments of nbslib   Mon Aug  8 17:02:58 EDT 1988
C  -- nbslib software is in the public domain or else copyrighted by a
C        U. S. government contractor which licenses the government to use,
C        reproduce, and distribute it.  Neither the United States, nor the
C        Department of Commerce, nor the National Bureau of Standards, nor
C        any of its employees, makes any warranty, express or implied, or
C        assumes any legal liability or reponsibility for the accuracy,
C        completeness, or usefulness of any information, apparatus, product,
C        process, or software disclosed or represents that its use would not
C        infringe privately owned rights.
C  ------------------------------------------------------------------------
C       From the Sequent Balance 8000 version, NBS, 2/12/86.
C       Last modified 2/12/86
C**********************************************************************
C        Sequent version GAMSIMP.FOR -- Gamma function benchmark
C        Language: ANSI FORTRAN X3.9-1978 (FORTRAN-77)
C
C       Requirements statement:
C       Given parameters of a subroutine (which is to be executed by
C       a number of processors concurrently on an MIMD machine),
C       invoke the target routine (for multiple execution).
C       Include timing calls to measure elapsed time of subroutine
C       execution.
C       Structure the driver to be easily updated for various methods
C       of invoking the target routine, including concurrent operation
C       of processes, one per invocation of the subroutine.  (Such
C       calls will include mechanisms for partitioning the problem
C       space among processes to minimize their interaction.)
C       The driver is to be run in batch mode and written in FORTRAN-77.
C
C        Algorithm:
C       1. Obtain parameters as follows:
C       N, N > 0, value for which Gamma (N) is to be computed.
C       MM, MM > 0, an integer, the number of intervals over which
C              integration by Simpson's rule is to be performed.
C       K, K > 0, the length of each of the MM intervals, K << MM.
C       MS, MS > 0, an even integer, the number of ordinates at which
C              Gamma (N) is to be evaluated in each interval of
C              length K.
C
C       2. Compute index I identifying lower bound of each of the MM
C       subintervals to be computed for each invocation of the target
C       routine, 0 =< I =< MM.
C
C       3. Initiate MM versions of the target routine either concurrently
C       or serially for a different value of I.  Note start time of
C       each routine.
C
C       4. Accumulate the sum of the MM integrations over each subinterval.
C       Note end time of each routine.
C       Note and report any error from each routine, including the
C       index value assigned the routine.
C
C       5. Produce the value of the resulting sum,
C                     value of total elapsed time of routines,
C                     recap values of input parameters.
C
C       Interfaces to other modules:
C
C       CALL GERRCK (K, MM, MS, N, OUTLUN, ERROR)
C              to verify parameters before starting timing run.
C
C       CALL SIMPSP (K, L, MS, N, SUM)
C              to evaluate Gamma (N) by integration by Simpson's 1/3
C              rule.
C
C       SECNDS (X)
C              to return the value of the time-of-day clock less X.
C
C Modifications:
C 2/12/86 KMD. Add prompt for input parameters.
        SUBROUTINE GAUSS
C
C        10 POINT GAUSSIAN QUADRATURE INTEGRATION
C
        REAL   AH,ANM,BET,E,R,Y(6),UP(5),VE,XL,XU,W,Z
        INTEGER I,JK
        COMMON / GAUSSX / R,BET,E,VE,AH,XU,XL,Y
        COMMON /LF/ LOOPFACTOR
C
C
C...... INITIALIZE VARIABLES
C
        R     = 24.84375
        BET   = 0.850
        E     = 0.0174533
        UP(1) = 41.67
        UP(2) = 138.33
        UP(3) = 221.67
        UP(4) = 318.33
        UP(5) = 360.00
        VE    = 13.5/R
        VE    = VE*VE
        AH    = 12.0/R
        AH    = AH*AH         
C
        DO 30 JK=1,(1000*LOOPFACTOR)
C
           Y(6) = 0.0
           XL   = 0.0
C
           DO 10 I=1,4,2
              XU = UP(I)
              CALL GAS(Z)
              Y(I) = Z
              XL   = XU
              XU   = UP(I+1)
              CALL GAST(W)
              Y(I+1) = W
              XL     = XU
10           CONTINUE
C
           XU = UP(5)
           CALL GAS(Z)
           Y(5) = Z
C
           DO 20 I=1,5
              Y(6) = Y(6) + Y(I)
20           CONTINUE
C
           ANM = (1.0 - BET*BET)/(6.28318536*BET)
           Y(6) = ANM*Y(6)
C
30        CONTINUE
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE GAS(W)
C
        COMMON / GAUSSX / R,BET,E,VE,AH,XU,XL,Y
        REAL R,BET,E,VE,AH,XU,XL,Y(6),A,B,C,AD,W
C
        EXTERNAL FUN
C
        A  = 0.5*(XU + XL)
        B  = XU - XL
        C  = 0.4869533*B
        AD = FUN(A+C,E,VE,BET) + FUN(A-C,E,VE,BET)
        W  = 0.03333567*AD
        C  = 0.4325317*B
        AD = FUN(A+C,E,VE,BET) + FUN(A-C,E,VE,BET)
        W  = W + 0.07472567*AD
        C  = 0.3397048 * B
        AD = FUN(A+C,E,VE,BET) + FUN(A-C,E,VE,BET)
        W  = W + 0.1095432*AD
        C  = 0.2166977*B
        AD = FUN(A+C,E,VE,BET) + FUN(A-C,E,VE,BET)
        W  = W + 0.1346334*AD
        C  = 0.07443717*B
        AD = FUN(A+C,E,VE,BET) + FUN(A-C,E,VE,BET)
        W  = B*E*(W + 0.1477621*AD)
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE GAST(W)
C
        COMMON / GAUSSX / R,BET,E,VE,AH,XU,XL,Y
        REAL R,BET,E,VE,AH,XU,XL,Y(6),A,B,C,AB,W
        EXTERNAL FTC
C
        A  = 0.5*(XU + XL)
        B  = XU - XL
        C  = 0.4869533*B
        AB = FTC(A+C,E,AH,BET) + FTC(A-C,E,AH,BET)
        W  = 0.03333567*AB
        C  = 0.4325317*B
        AB = FTC(A+C,E,AH,BET) + FTC(A-C,E,AH,BET)
        W  = W + 0.07472567*AB
        C  = 0.3397048*B
        AB = FTC(A+C,E,AH,BET) + FTC(A-C,E,AH,BET)
        W  = W + 0.1095432*AB
        C  = 0.2166977*B
        AB = FTC(A+C,E,AH,BET) + FTC(A-C,E,AH,BET)
        W  = W + 0.1346334*AB
        C  = 0.07443717*B
        AB = FTC(A+C,E,AH,BET) + FTC(A-C,E,AH,BET)
        W  = B*E*(W + 0.1477621*AB)
        RETURN
        END
C
C*****************************************************************************
C
        FUNCTION FUN(X,E,VE,BET)
C
        REAL X,E,VE,BET,F,FUN
C
        F = COS(E*X)
        FUN = 1.0/(1.0 - BET*F/SQRT(F*F + VE))
        RETURN
        END
C
C*****************************************************************************
C
        FUNCTION FTC(X,E,AH,BET)
C
        REAL X,E,AH,BET,G,FTC
C
        G = SIN(E*X)
        FTC = 1.0/(1.0 - BET*G/SQRT(G*G + AH))
        RETURN
        END
        SUBROUTINE HANOI
C
        IMPLICIT INTEGER (A-Z)
        COMMON /HAN/ NO(20),NF(20),NT(20)
        COMMON /LF/ LOOPFACTOR
C
        DO 10 LP=1,(1*LOOPFACTOR)
          I = 1
          K = 2
          N = 20
          JP = 1
3         IF (N.EQ.1) GO TO 7
          NO(JP) = N
          NF(JP) = I
          NT(JP) = K
          K = 6-I-K
6         N = N-1
          JP = JP+1
          GO TO 3
7         JP=JP-1
          IF (JP.EQ.0) GO TO 10
8         N=NO(JP)
          IF (N.LE.0) GO TO 7
          I=NF(JP)
          K=NT(JP)
          NO(JP)=-N
          I=6-I-K
          GO TO 6
10        CONTINUE
C
      RETURN
      END
      SUBROUTINE INIT(NAME)
C
      CHARACTER*6 NAME

      INTEGER E, F
C-----This block applies to all subroutines.
      COMMON /LF/ LOOPFACTOR
      
C-----This block is for subroutines from the Livermore Loops.
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
                               
C-----This block is for the Alam18 subroutine.
      COMMON /ALAMOS18/ NCASE(10), AA(5000), BB(5000), CC(5000),
     + DD(5000), RR(5000), JJ(5000)

C-----This block for Philco subroutine.
      DOUBLE PRECISION CCM2, TEMP7
      COMMON /PHILC/ CCM2(12,12),TEMP7(12,12)

C-----These blocks for the LU decomposition routines (LUDD, LUSD)
      COMMON /LUSDC/ LUSA(301,300),LUSX(301),LUSB(301)
      DOUBLE PRECISION LUDA(301,300),LUDX(301),LUDB(301)
      COMMON /LUDDC/ LUDA, LUDX, LUDB

C-----This block is for the LINPAC kernel.
      DOUBLE PRECISION LINA(201,200),LINB(200)
      COMMON /LINPACC/ LINA, LINB

      DATA NCASE(1), NCASE(2), NCASE(3), NCASE(4), NCASE(5),
     + NCASE(6), NCASE(7), NCASE(8), NCASE(9), NCASE(10)
     + /20,100,500,1000,2000,4000,5000,10000,50000,100000/

C
C	Split into two statements so as not to exceed f77 standard
C	limit on continuation lines.

        DATA ((TEMP7(i1, j1), i1=1,12), j1=1,6)
     +   /18318.1,-2374.5,14.9,10862.5,-623.6,327.6 ,24421.7,5084.9,
     +  758.2, 22163.9,-11529.6,-1468.6 ,-2374.5,392.1,-13.8,-1348.9,
     +  180.8,-7.6, -3106.7,-528.0,-164.8,-2732.6,1547.0,123.6 ,14.9,
     +  -13.8,197.5, 22.7,3.7,3.1 ,47.6,15.2,-3.1,47.1,-18.9,-7.2 ,
     +  10862.5,-1348.9, 22.7,6718.4,-315.,221.9 ,14894.3,3188.9,402.2,
     +  13588.9,-6963.1, -928.1 ,-623.6,180.8,3.7,-315.0,171.5,54.4 ,
     +  -771.5,-12.3,-138.4, -612.5,470.3,-60.3 ,327.6,-7.6,3.1,221.9,
     +  54.4,60. ,476.,153.9, -31.6,467.5,-186.9,-72.1 / 

        DATA ((TEMP7(i1, j1), i1=1,12), j1=7,12)
     +   /24421.7,-3106.7,47.6,14894.3,-771.5,476. ,33517.8,7071.6,
     +  1003.4,30401.1,-15707.3,-2059.9 ,5084.9,-528.,15.2,3188.9,-12.3,
     +  153.9 ,7071.6,1693.,62.2,6549.7,-3202.3,-527.7 ,758.2,-164.8,
     +  -3.1,402.2,-138.4,-31.6 ,1003.4,62.2,233.2,847.1,-542.3,26.4 ,
     +  22163.9,-2732.6,47.1,13588.9,-612.5,467.5 ,30401.1,6549.7,847.1,
     +  27846.6,-14237.4,-1944.2 ,-11529.6,1547.,-18.9,-6963.1,470.3,
     +  -186.9 ,-15707.3,-3202.3,-542.3,-14237.4,7498.1,864.4 ,-1468.6,
     +  123.6,-7.2,-928.1,-60.3,-72.1 ,-2059.9,-527.7,26.4,-1944.2,
     +  864.4,301.3/ 

C-----Now do whatever initializations are required for whatever benchmark
C-----subroutine we're about to do.
C
C-----In general,we make any random values very close to 1.0 to try and
C-----minimize the possibility of floating overflow.
      IF (NAME.EQ.'AIRREL') THEN
        LOOPFACTOR = 7
C                                         
      ELSE IF (NAME.EQ.'ALAM18') THEN                       
        LOOPFACTOR = 1 
        NLOOP = 7
        K = NLOOP + 9
        N = 5000        
        DO 10 I=1,N
          AA(I) = I
          BB(I) = AA(I)**.5
          CC(I) = (I-1)**2
          FI = I
          DD(I) = 1.0/FI
          JJ(I) = I - K
   10   CONTINUE
C
      ELSE IF (NAME.EQ.'DOUBLE') THEN
        LOOPFACTOR = 5
C
      ELSE IF (NAME.EQ.'DWHET ') THEN
        LOOPFACTOR = 5
C
      ELSE IF (NAME.EQ.'EGYPT ') THEN
        LOOPFACTOR = 14
C                 
      ELSE IF (NAME.EQ.'EGYPTD') THEN
        LOOPFACTOR = 11 
C
      ELSE IF (NAME.EQ.'FFTR  ') THEN
        LOOPFACTOR = 9
C
      ELSE IF (NAME.EQ.'GAMSIM') THEN
        LOOPFACTOR = 18
C
      ELSE IF (NAME.EQ.'GAUSS ') THEN
        LOOPFACTOR = 7
C
      ELSE IF (NAME.EQ.'HANOI ') THEN
        LOOPFACTOR=6
C
      ELSE IF (NAME.EQ.'LA1DPP') THEN
        LOOPFACTOR = 73
        DO 20 J=1,150
          VX(J) = 1.0 - 1.0/(J+12)
          XX(J) = 1.0 - 1.0/(J+13)
          GRD(J) = MOD(J*51027,65)+1
C          VX(J) = J/1E10
C          XX(J) = J/1E10
C          GRD(J) = J
  20    CONTINUE
        DO 30 J=1,67
C These were all FLOAT(J)
C          RH(J) = J/1E10
C          DEX(J) = J/1E10
C          EX(J) = -VX(J)
          RH(J) = FLOAT(J)
          DEX(J) = 1.0 - 1.0/(J+10)
          EX(J) = 1.0 - 1.0/(J+11)
  30   CONTINUE
C
      ELSE IF (NAME.EQ.'LA2DHY') THEN
        LOOPFACTOR=28
        DO 40 j= 1,7    
        DO 40 k= 1,101
C make these very close to 1 instead of close to zero
          za(k,j) = J/1E10 + 1.0
          zb(k,j) = J/1E10 + 1.0
          zr(k,j) = J/1E10 + 1.0
          zu(k,j) = J/1E10 + 1.0
          zv(k,j) = J/1E10 + 1.0
          zz(k,j) = J/1E10 + 1.0
   40   CONTINUE
C
      ELSE IF (NAME.EQ.'LA2DPP') THEN
        LOOPFACTOR = 94
        DO 50 K=1,512     
C------All computations are from MIPS' modifications to the data set
              P(1,K) = K/32.0 + 16.0
              P(2,K) = K/256.0 + 2.0
              P(3,K) = 0.0  
              P(4,K) = 0.0
   50   CONTINUE
        DO 60 K=1,1001
          Y(K) = 0.00000321
          Z(K) = 0.00000321
   60   CONTINUE
        DO 70 J=1,192
          E(J) = (32-(J-32))/3
          F(J) = (4-(J-32))/3
   70   CONTINUE            
        DO 80 J=1,64                                     
          DO 80 K=1,8              
              BR = (32.5-J)/80.0
              CR = (4.5-K)/10.0
              B(J,K) = BR*ABS(BR)
              C(J,K) = CR*ABS(CR)
              H(J,K) = 0.0
   80   CONTINUE
C
      ELSE IF (NAME.EQ.'LABLE') THEN
        LOOPFACTOR = 126      
        DO 90 K=1,1001
          X(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Y(K) = ((K/K+1.0)-0.5)/1E5+1.0
   90   CONTINUE                       
C
      ELSE IF (NAME.EQ.'LADIFF') THEN
        LOOPFACTOR = 5
        DO 100 K=1,15
          DO 100 L=1,100  
            PX(K,L) = ((K/K+1.0)-0.5)/1E5+1.0
            CX(K,L) = ((K/K+1.0)-0.5)/1E5+1.0
  100     CONTINUE
C
      ELSE IF (NAME.EQ.'LAFIMI') THEN
        LOOPFACTOR=25
        DO 110 K=1,1001
          X(K) = ((K/K+1.0)-0.5)/1E5+1.0
  110   CONTINUE
C
      ELSE IF (NAME.EQ.'LAFSUM') THEN
        LOOPFACTOR = 139       
        DO 120 K=1,1001
          X(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Y(K) = ((K/K+1.0)-0.5)/1E5+1.0
  120   CONTINUE
C
      ELSE IF (NAME.EQ.'LAHYDR') THEN
        LOOPFACTOR = 12
        DO 130 K=1,1001
          X(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Y(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Z(K) = ((K/K+1.0)-0.5)/1E5+1.0
  130   CONTINUE
C
      ELSE IF (NAME.EQ.'LAIMCO') THEN
        LOOPFACTOR = 44
        DO 140 j=1,101
          VE3(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VLIN(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VLR(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VSP(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VSTP(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VXND(j) = ((J/J+1.0)-0.5)/1E5+1.0
          VXNE(j) = ((J/J+1.0)-0.5)/1E5+1.0
  140   CONTINUE
C
      ELSE IF (NAME.EQ.'LAMAPR') THEN
        LOOPFACTOR=10         
        DO 150 j= 1,101
          DO 150 k= 1,25
            cx(k,j) = ((K/K+1.0)-0.5)/1E5+1.0
            px(k,j) = ((K/K+1.0)-0.5)/1E5+1.0
            vy(j,k) = ((K/K+1.0)-0.5)/1E5+1.0
  150   CONTINUE
C
      ELSE IF (NAME.EQ.'LAPDE') THEN
        LOOPFACTOR = 204
        DO 160 J=1,5
          DO 160 K=1,22        
            DO 160 L=1,2
              U1(J,K,L) = ((K/K+1.0)-0.5)/1E5+1.0
              U2(J,K,L) = ((K/K+1.0)-0.5)/1E5+1.0
              U3(J,K,L) = ((K/K+1.0)-0.5)/1E5+1.0
  160       CONTINUE  
C
      ELSE IF (NAME.EQ.'LATRD1') THEN
        LOOPFACTOR = 84
        DO 170 K=1,1001
          X(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Y(K) = ((K/K+1.0)-0.5)/1E5+1.0
          Z(K) = ((K/K+1.0)-0.5)/1E5+1.0
  170   CONTINUE
C
      ELSE IF (NAME.EQ.'LATRD2') THEN
        LOOPFACTOR=84
        DO 180 K=1,1001
C-----Make these very small rather than close to 1.0;
          X(K) = ((K/K+1.0)-0.5)/1E5
          Z(K) = ((K/K+1.0)-0.5)/1E5
  180   CONTINUE
C
C-----In the following three initializations values aren't important as
C-----they call the MATGEN subroutine themselves to set up data values.
C-----We just do a dummy initialization here to get them in memory so we
C-----don't pagefault at MATGEN (which is inside the timing routine).
C
      ELSE IF (NAME.EQ.'LINPAC') THEN
        LOOPFACTOR = 13
        DO 182 J = 1,200
          LINB(J) = 0.0
          DO 181 I = 1,201
            LINA(I,J) = 0.0
  181     CONTINUE
  182   CONTINUE
C
      ELSE IF (NAME.EQ.'LUDD  ') THEN
        LOOPFACTOR = 135
        DO 184 J = 1,LOOPFACTOR
          LUDB(J) = 0.0
          DO 183 I = 1,LOOPFACTOR
            LUDA(I,J) = 0.0
  183     CONTINUE
  184   CONTINUE
C                  
      ELSE IF (NAME.EQ.'LUSD  ') THEN
        LOOPFACTOR = 118
        DO 186 J = 1,LOOPFACTOR
          LUSB(J) = 0.0
          DO 185 I = 1,LOOPFACTOR
            LUSA(I,J) = 0.0
  185     CONTINUE
  186   CONTINUE
C
      ELSE IF (NAME.EQ.'NEFF  ') THEN
        LOOPFACTOR = 7
C
      ELSE IF (NAME.EQ.'PHILCO') THEN
        LOOPFACTOR=7
        DO 190 I=1,12
          DO 190 J=1,12
            CCM2(I,J) = TEMP7(J,I)
190       CONTINUE
C
      ELSE IF (NAME.EQ.'PRIME ') THEN
        LOOPFACTOR = 4
C
      ELSE IF (NAME.EQ.'RR2   ') THEN
        LOOPFACTOR = 16
C
      ELSE IF (NAME.EQ.'RR3   ') THEN
        LOOPFACTOR = 12
C
      ELSE IF (NAME.EQ.'RR4   ') THEN
        LOOPFACTOR = 21
C
      ELSE IF (NAME.EQ.'SINGLE') THEN
        LOOPFACTOR=7          
C
      ELSE IF (NAME.EQ.'SWHET ') THEN
        LOOPFACTOR=7
C
      END IF

      RETURN
      END
                    
       SUBROUTINE LA1DPP
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C***********************************************************************
C***********************************************************************
C***  LOOP 14  1-D PARTICLE PUSHER
      DO 520 KK=1,(100*LOOPFACTOR)
C          CALL INIT -- This is now called from the main program
C-----Since we're still having trouble getting well-behaved data for this
C-----thing, call INIT periodically to reset data to resonable values.
        IF (MOD(KK,15).EQ.0) CALL INIT('LA1DPP')
          FLX = 0.00000001
          DO 510 K=1,150
              IX = GRD(K)
              XI = IX
              VX(K) = VX(K) + EX(IX) + (XX(K)-XI)*DEX(IX)
              XX(K) = XX(K) + VX(K) + FLX
              IR = XX(K)  
              RI = IR
              RX1 = XX(K) - RI
              IR = IR - (IR/64)*64 + 64
              XX(K) = RI + RX1
              RH(IR) = RH(IR) + 1.0 - RX1
              RH(IR+1) = RH(IR+1) + RX1
  510     CONTINUE
  520 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LA2DHY
C     *********************************************************
C     *  2-D IMPLICIT HYDRODYNAMICS FRAGMENT (LLNL Kernel 23) *
C     *********************************************************

      COMMON /LF/ LOOPFACTOR                  
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
                      
      n = 101
      fw= 0.000175D0
      DO 23  L= 1,(100*LOOPFACTOR)
        DO 23  j= 2,6
          DO 23  k= 2,n
            QA= ZA(k,j+1)*ZR(k,j) +ZA(k,j-1)*ZB(k,j) +
     .          ZA(k+1,j)*ZU(k,j) +ZA(k-1,j)*ZV(k,j) +ZZ(k,j)
   23 ZA(k,j)= ZA(k,j) +fw*(QA -ZA(k,j))

      RETURN
      END
      SUBROUTINE LA2DPP                        
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C***********************************************************************
C***********************************************************************
C***  LOOP 13  2-D PARTICLE PUSHER
      DO 470 KK=1,(100*LOOPFACTOR)
C         CALL INIT--Now called from main program
C-----Since we're still having trouble getting well-behaved data for this
C-----thing, call INIT periodically to reset data to resonable values.
         IF (MOD(KK,500).EQ.0) CALL INIT('LA2DPP')
         DO 460 IP=1,128
              I1 = P(1,IP)                       
              J1 = P(2,IP)
              P(3,IP) = P(3,IP) + B(I1,J1)
              P(4,IP) = P(4,IP) + C(I1,J1)
              P(1,IP) = P(1,IP) + P(3,IP)
              P(2,IP) = P(2,IP) + P(4,IP)
              I2 = P(1,IP)
              J2 = P(2,IP)
              P(1,IP) = P(1,IP) + Y(I2+32)
              P(2,IP) = P(2,IP) + Z(J2+32)
              I2 = I2 + E(I2+32)
              J2 = J2 + F(J2+32)
              H(I2,J2) = H(I2,J2) + 1.0
  460     CONTINUE
  470 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LABLE
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C
C***********************************************************************
C***********************************************************************
C***  LOOP 4   BANDED LINEAR EQUATIONS
C      X100 = XEX
      DO 130 KK=1,(100*LOOPFACTOR)               
C         CALL INIT--Now called from main program
         DO 120 L=7,107,50
            LW = L
            DO 110 J=30,870,5
              X(L-1) = X(L-1) - X(LW)*Y(J)
              LW = LW + 1
  110       CONTINUE
            X(L-1) = Y(5)*X(L-1)
  120    CONTINUE
  130 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LADIFF
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C***********************************************************************
C
      DO 1000 ILOOP = 1,(10*LOOPFACTOR)
C***********************************************************************
C***  LOOP 10 DIFFERENCE PREDICTORS EXCERPT
      DO 10 L=1,100
          DO 10 I=1,400
              AR = CX(5,I) 
              BR = AR - PX(5,I)
              PX(5,I) = AR
              CR = BR - PX(6,I)
              PX(6,I) = BR
              AR = CR - PX(7,I)
              PX(7,I) = CR
              BR = AR - PX(8,I)
              PX(8,I) = AR
              CR = BR - PX(9,I)
              PX(9,I) = BR
              AR = CR - PX(10,I)
              PX(10,I) = CR
              BR = AR - PX(11,I)
              PX(11,I)= AR
              CR = BR - PX(12,I)
              PX(12,I) = BR
              PX(14,I) = CR - PX(13,I)
              PX(13,I) = CR
   10   CONTINUE
 1000 CONTINUE
C
      RETURN
      END
      SUBROUTINE LAFIMI
C     ************************************************************
C     * FIND LOCATION OF FIRST MINIMUM IN ARRAY (LLNL Kernel 24) *
C     ************************************************************

      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)

      n = 1001
      X( n/2)= -1.000D+10
      DO 24  L= 1,(500*LOOPFACTOR)
        m= 1
        DO 24  k= 2,n
          IF( X(k).LT.X(m))  m= k
   24 CONTINUE
      RETURN
      END
      SUBROUTINE LAFSUM
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C***********************************************************************
C***********************************************************************
C***  LOOP 11  FIRST SUMS
      X(1) = Y(1)
      DO 390 KK=1,(100*LOOPFACTOR)
C          CALL INIT--Now called from main program.
          DO 380 K=2,1000
              X(K) = X(K-1) + Y(K)                 
  380     CONTINUE
  390 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LAHYDR                 
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C
C***********************************************************************
C
      R = 1.0000001
      T = 1.0000001
      DO 1000 ILOOP = 1,(10*LOOPFACTOR)
C***********************************************************************
C***  LOOP 1   HYDRO EXCERPT
      Q = 0.0
      DO 30 KK=1,100
          DO 20 K=1,400
              X(K) = Q + Y(K)*(R*Z(K+10)+T*Z(K+11))
   20     CONTINUE
   30 CONTINUE
C***********************************************************************
 1000 CONTINUE
      RETURN
      END
      SUBROUTINE LAIMCO
C     *******************************************************
C     * IMPLICIT, CONDITIONAL COMPUTATION (LLNL Kernel 17)  *
C     *******************************************************
C         RECURSIVE-DOUBLING VECTOR TECHNIQUES CAN NOT BE USED
C         BECAUSE CONDITIONAL OPERATIONS APPLY TO EACH ELEMENT.

      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
      dw = 5.0000D0/3.0000D0
      fw = 1.0000D0/3.0000D0
      tw = 1.0300D0/3.0700D0
      n = 101
      DO 62 L = 1,(500*LOOPFACTOR)
        i = n
        j = 1
        INK = -1
        SCALE = dw
        XNM = fw
        E6 = tw
        GO TO 61
C                                            STEP MODEL
  60    E6 = XNM*VSP(i)+VSTP(i)
        VXNE(i) = E6
        XNM = E6
        VE3(i) = E6
        i = i+INK
        IF( i.EQ.j) GO TO  62
  61      E3 = XNM*VLR(i) +VLIN(i)
          XNEI = VXNE(i)
          VXND(i) = E6 + L/1E10
          XNC = SCALE*E3
C                                            SELECT MODEL
          IF( XNM .GT.XNC) GO TO  60
          IF( XNEI.GT.XNC) GO TO  60
C                                            LINEAR MODEL
          VE3(i) = E3
          E6 = E3+E3-XNM
          VXNE(i) = E3+E3-XNEI
          XNM = E6
          i = i+INK
          IF( i.NE.j) GO TO 61
   62 CONTINUE
      RETURN
      END

      SUBROUTINE LAMAPR
C     ******************************************
C     * MATRIX*MATRIX PRODUCT (LLNL Kernel 21) *
C     ******************************************

      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
      n = 101
      DO 21 L= 1,(10*LOOPFACTOR)
        DO 21 k= 1,25
          DO 21 i= 1,25
            DO 21 j= 1,n
              PX(i,j)= PX(i,j) +VY(i,k) * CX(k,j)
   21 CONTINUE
      RETURN
      END

      SUBROUTINE LAPDE
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C
C***********************************************************************
C***********************************************************************
C***  LOOP 8   P.D.E. INTEGRATION
      NL1 = 1          
      NL2 = 2
C These inits added in DWH 7 October 1988
      SIG = 1.00000001
      A11 = 1.00000001
      A21 = 1.00000001
      A31 = 1.00000001
      A12 = 1.00000001
      A22 = 1.00000001
      A32 = 1.00000001
      A13 = 1.00000001
      A23 = 1.00000001
      A33 = 1.00000001
      DO 270 KK=1,(100*LOOPFACTOR)
C          CALL INIT--Now called from main program
          DO 260 KX=2,3
              DO 250 KY=2,21
                 DU1 = U1(KX,KY+1,NL1) - U1(KX,KY-1,NL1)
                 DU2 = U2(KX,KY+1,NL1) - U2(KX,KY-1,NL1)
                 DU3 = U3(KX,KY+1,NL1) - U3(KX,KY-1,NL1)
                 U1(KX,KY,NL2) = U1(KX,KY,NL1) + A11*DU1 + A12*DU2 +
     *            A13*DU3 + SIG*(U1(KX+1,KY,NL1)-2.0*U1(KX,KY,NL1)
     *            +U1(KX-1,KY,NL1))
                 U2(KX,KY,NL2) = U2(KX,KY,NL1) + A21*DU1 + A22*DU2 +
     *            A23*DU3 + SIG*(U2(KX+1,KY,NL1)-2.0*U2(KX,KY,NL1)
     *            +U2(KX-1,KY,NL1))
                 U3(KX,KY,NL2) = U3(KX,KY,NL1) + A31*DU1 + A32*DU2 +
     *            A33*DU3 + SIG*(U3(KX+1,KY,NL1)-2.0*U3(KX,KY,NL1)
     *            +U3(KX-1,KY,NL1))
  250         CONTINUE
  260     CONTINUE
  270 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LATRD1
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C
C***********************************************************************
C***********************************************************************
C***  LOOP 5   TRI-DIAGONAL ELIMINATION, BELOW DIAGONAL
      DO 160 KK=1,(100*LOOPFACTOR)
C         CALL INIT--Now called from main program
         DO 150 I=2,998,3
            X(I) = Z(I)*(Y(I)-X(I-1))
            X(I+1) = Z(I+1)*(Y(I+1)-X(I))
            X(I+2) = Z(I+2)*(Y(I+2)-X(I+1))
  150    CONTINUE
  160 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LATRD2            
C
      INTEGER E, F
      COMMON /LF/ LOOPFACTOR
      COMMON /LAWLIV/ B(64,8), BNK1(5), BNK2(5), BNK3(5), BNK4(5),
     +  BNK5(5), C(64,8), CX(25,101), D(300), DEX(67), E(192), EX(67),
     +  F(192), GRD(150), H(64,8), IZONE(300), P(4,512), PLAN(300),
     +  PX(25,101), RH(67), U(1001), U1(5,22,2), U2(5,22,2), U3(5,22,2),
     +  VE3(101), VLIN(101), VLR(101), VSP(101), VSTP(101), VX(150),
     +  VXND(101), VXNE(101), VY(101,25), X(1001), XX(150), Y(1001),
     +  Z(1001), ZA(101,7), ZB(101,7), ZR(101,7), ZU(101,7), ZV(101,7),
     +  ZZ(101,7)
C
C
C***********************************************************************
C***********************************************************************
C***  LOOP 6   TRI-DIAGONAL ELIMINATION, ABOVE DIAGONAL
      DO 190 KK=1,(100*LOOPFACTOR)
C         CALL INIT--Now called from main program.
         DO 180 J=3,999,3
            I = 1000 - J + 3
            X(I) = X(I) - Z(I)*X(I+1)
            X(I-1) = X(I-1) - Z(I-1)*X(I)
            X(I-2) = X(I-2) - Z(I-2)*X(I-1)
  180    CONTINUE
  190 CONTINUE
C***********************************************************************
      RETURN
      END
      SUBROUTINE LINPAC
      DOUBLE PRECISION LINA(201,200),LINB(200)
      COMMON /LINPACC/ LINA, LINB
      DOUBLE PRECISION NORMA
      INTEGER IPVT(200)
      COMMON /LF/ LOOPFACTOR
      LDA = 201
C
      N = 100
C
      DO 10 I=1,(1*LOOPFACTOR)
         CALL MATGEN(LINA,LDA,N,LINB,NORMA)
         CALL DGEFA(LINA,LDA,N,IPVT,INFO)
         CALL DGESL(LINA,LDA,N,IPVT,LINB,0)
   10 CONTINUE
C
      RETURN
      END

      SUBROUTINE MATGEN(A,LDA,N,B,NORMA)
      DOUBLE PRECISION A(LDA,1),B(1),NORMA
C
      INIT = 1325
      NORMA = 0.0
      DO 30 J = 1,N
         DO 20 I = 1,N
            INIT = MOD(3125*INIT,65536)
            A(I,J) = (INIT - 32768.0)/16384.0
            NORMA = DMAX1(A(I,J), NORMA)
   20    CONTINUE
   30 CONTINUE
      DO 35 I = 1,N
          B(I) = 0.0
   35 CONTINUE
      DO 50 J = 1,N
         DO 40 I = 1,N
            B(I) = B(I) + A(I,J)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
      INTEGER LDA,N,IPVT(1),INFO
      DOUBLE PRECISION A(LDA,1)
C
C     DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
C
C     DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C     (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
C                     IF CALLED.  USE  RCOND  IN DGECO FOR A RELIABLE
C                     INDICATION OF SINGULARITY.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DSCAL,IDAMAX
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION T
      INTEGER IDAMAX,J,K,KP1,L,NM1
C
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      INFO = 0
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 70
      DO 60 K = 1, NM1
         KP1 = K + 1
C
C        FIND L = PIVOT INDEX
C
         L = IDAMAX(N-K+1,A(K,K),1) + K - 1
         IPVT(K) = L
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C           INTERCHANGE IF NECESSARY
C
            IF (L .EQ. K) GO TO 10
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   10       CONTINUE
C
C           COMPUTE MULTIPLIERS
C
            T = -1.0D0/A(K,K)
            CALL DSCAL(N-K,T,A(K+1,K),1)
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C
            DO 30 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 20
                  A(L,J) = A(K,J)
                  A(K,J) = T
   20          CONTINUE
               CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
   30       CONTINUE
         GO TO 50
   40    CONTINUE
            INFO = K   
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D0) INFO = N
      RETURN
      END
      SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)
      INTEGER LDA,N,IPVT(1),JOB
      DOUBLE PRECISION A(LDA,1),B(1)
C
C     DGESL SOLVES THE DOUBLE PRECISION SYSTEM
C     A * X = B  OR  TRANS(A) * X = B
C     USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE OUTPUT FROM DGECO OR DGEFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM DGECO OR DGEFA.
C
C        B       DOUBLE PRECISION(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C                = NONZERO   TO SOLVE  TRANS(A)*X = B  WHERE
C                            TRANS(A)  IS THE TRANSPOSE.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
C        CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
C        OR DGEFA HAS SET INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
C           IF (RCOND IS TOO SMALL) GO TO ...
C           DO 10 J = 1, P
C              CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
C        10 CONTINUE
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DDOT
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DDOT,T
      INTEGER K,KB,L,NM1
C
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE  L*Y = B
C
         IF (NM1 .LT. 1) GO TO 30
         DO 20 K = 1, NM1
            L = IPVT(K)
            T = B(L)
            IF (L .EQ. K) GO TO 10
               B(L) = B(K)
               B(K) = T
   10       CONTINUE
            CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
   20    CONTINUE
   30    CONTINUE
C
C        NOW SOLVE  U*X = Y
C
         DO 40 KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/A(K,K)
            T = -B(K)
            CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
   40    CONTINUE
      GO TO 100
   50 CONTINUE
C
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
C        FIRST SOLVE  TRANS(U)*Y = B
C
         DO 60 K = 1, N
            T = DDOT(K-1,A(1,K),1,B(1),1)
            B(K) = (B(K) - T)/A(K,K)
   60    CONTINUE
C
C        NOW SOLVE TRANS(L)*X = Y
C
         IF (NM1 .LT. 1) GO TO 90
         DO 80 KB = 1, NM1
            K = N - KB
            B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
            L = IPVT(K)
            IF (L .EQ. K) GO TO 70
               T = B(L)
               B(L) = B(K)
               B(K) = T
   70       CONTINUE
   80    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DA
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF (DA .EQ. 0.0D0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DY(IY) = DY(IY) + DA*DX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        DY(I) = DY(I) + DA*DX(I)
        DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
        DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
        DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
   50 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      DDOT = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DTEMP + DX(IX)*DY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      DDOT = DTEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DX(I)*DY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
     *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
   50 CONTINUE
   60 DDOT = DTEMP
      RETURN
      END
      SUBROUTINE  DSCAL(N,DA,DX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DA,DX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        DX(I) = DA*DX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DX(I) = DA*DX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        DX(I) = DA*DX(I)
        DX(I + 1) = DA*DX(I + 1)
        DX(I + 2) = DA*DX(I + 2)
        DX(I + 3) = DA*DX(I + 3)
        DX(I + 4) = DA*DX(I + 4)
   50 CONTINUE
      RETURN
      END
      INTEGER FUNCTION IDAMAX(N,DX,INCX)
C
C     FINDS THE INDEX OF ELEMENT HAVING MAX. DABSOLUTE VALUE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DMAX
      INTEGER I,INCX,IX,N
C
      IDAMAX = 0
      IF( N .LT. 1 ) RETURN
      IDAMAX = 1
      IF(N.EQ.1)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      IX = 1
      DMAX = DABS(DX(1))
      IX = IX + INCX
      DO 10 I = 2,N
         IF(DABS(DX(IX)).LE.DMAX) GO TO 5
         IDAMAX = I
         DMAX = DABS(DX(IX))
    5    IX = IX + INCX
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
   20 DMAX = DABS(DX(1))
      DO 30 I = 2,N
         IF(DABS(DX(I)).LE.DMAX) GO TO 30
         IDAMAX = I
         DMAX = DABS(DX(I))
   30 CONTINUE
      RETURN
      END
      SUBROUTINE MM (A, LDA, N1, N3, B, LDB, N2, C, LDC)
      DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*)
C
C   PURPOSE:
C     MULTIPLY MATRIX B TIMES MATRIX C AND STORE THE RESULT IN MATRIX A.
C
C   PARAMETERS:
C
C     A DOUBLE PRECISION(LDA,N3), MATRIX OF N1 ROWS AND N3 COLUMNS
C
C     LDA INTEGER, LEADING DIMENSION OF ARRAY A
C
C     N1 INTEGER, NUMBER OF ROWS IN MATRICES A AND B
C
C     N3 INTEGER, NUMBER OF COLUMNS IN MATRICES A AND C
C
C     B DOUBLE PRECISION(LDB,N2), MATRIX OF N1 ROWS AND N2 COLUMNS
C
C     LDB INTEGER, LEADING DIMENSION OF ARRAY B
C
C     N2 INTEGER, NUMBER OF COLUMNS IN MATRIX B, AND NUMBER OF ROWS IN
C         MATRIX C
C
C     C DOUBLE PRECISION(LDC,N3), MATRIX OF N2 ROWS AND N3 COLUMNS
C
C     LDC INTEGER, LEADING DIMENSION OF ARRAY C
C
C ----------------------------------------------------------------------
C
      DO 20 J = 1, N3
         DO 10 I = 1, N1
            A(I,J) = 0.0
   10    CONTINUE
         CALL DMXPY (N2,A(1,J),N1,LDB,C(1,J),B)
   20 CONTINUE
C
      RETURN
      END
      SUBROUTINE DMXPY (N1, Y, N2, LDM, X, M)
      DOUBLE PRECISION Y(*), X(*), M(LDM,*)
C
C   PURPOSE:
C     MULTIPLY MATRIX M TIMES VECTOR X AND ADD THE RESULT TO VECTOR Y.
C
C   PARAMETERS:
C
C     N1 INTEGER, NUMBER OF ELEMENTS IN VECTOR Y, AND NUMBER OF ROWS IN
C         MATRIX M
C
C     Y DOUBLE PRECISION(N1), VECTOR OF LENGTH N1 TO WHICH IS ADDED 
C         THE PRODUCT M*X
C
C     N2 INTEGER, NUMBER OF ELEMENTS IN VECTOR X, AND NUMBER OF COLUMNS
C         IN MATRIX M
C
C     LDM INTEGER, LEADING DIMENSION OF ARRAY M
C
C     X DOUBLE PRECISION(N2), VECTOR OF LENGTH N2
C
C     M DOUBLE PRECISION(LDM,N2), MATRIX OF N1 ROWS AND N2 COLUMNS
C
C ----------------------------------------------------------------------
C
C   CLEANUP ODD VECTOR
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   CLEANUP ODD GROUP OF TWO VECTORS
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   CLEANUP ODD GROUP OF FOUR VECTORS
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   30    CONTINUE
      ENDIF
C
C   CLEANUP ODD GROUP OF EIGHT VECTORS
C
      J = MOD(N2,16)
      IF (J .GE. 8) THEN
         DO 40 I = 1, N1
            Y(I) = ((((((( (Y(I))
     $             + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6))
     $             + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   40    CONTINUE
      ENDIF
C
C   MAIN LOOP - GROUPS OF SIXTEEN VECTORS
C
      JMIN = J+16
      DO 60 J = JMIN, N2, 16
         DO 50 I = 1, N1
            Y(I) = ((((((((((((((( (Y(I))
     $             + X(J-15)*M(I,J-15)) + X(J-14)*M(I,J-14))
     $             + X(J-13)*M(I,J-13)) + X(J-12)*M(I,J-12))
     $             + X(J-11)*M(I,J-11)) + X(J-10)*M(I,J-10))
     $             + X(J- 9)*M(I,J- 9)) + X(J- 8)*M(I,J- 8))
     $             + X(J- 7)*M(I,J- 7)) + X(J- 6)*M(I,J- 6))
     $             + X(J- 5)*M(I,J- 5)) + X(J- 4)*M(I,J- 4))
     $             + X(J- 3)*M(I,J- 3)) + X(J- 2)*M(I,J- 2))
     $             + X(J- 1)*M(I,J- 1)) + X(J)   *M(I,J)
   50    CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE LUDD

      DOUBLE PRECISION LUDA(301,300),X(301),LUDB(301)
      COMMON /LUDDC/ LUDA, X, LUDB
      COMMON /LF/ LOOPFACTOR
      INTEGER IPVT(300)
C
      LDA = 301
C
C      DO 200 N = 50,300,50
      DO 200 N = 70,110,20
C      N = LOOPFACTOR
         CALL DMATGEN(LUDA,LDA,N,LUDB)
         CALL DLU1(LUDA,LDA,N,IPVT,IERR)
         CALL DLUS(LUDA,LDA,N,IPVT,X,LUDB)
C
         CALL DMATGEN(LUDA,LDA,N,LUDB)
         CALL DLU2(LUDA,LDA,N,IPVT,IERR)
         CALL DLUS(LUDA,LDA,N,IPVT,X,LUDB)
C
         CALL DMATGEN(LUDA,LDA,N,LUDB)
         CALL DLU4(LUDA,LDA,N,IPVT,IERR)
         CALL DLUS(LUDA,LDA,N,IPVT,X,LUDB)
C
         CALL DMATGEN(LUDA,LDA,N,LUDB)
         CALL DLU8(LUDA,LDA,N,IPVT,IERR)
         CALL DLUS(LUDA,LDA,N,IPVT,X,LUDB)
C
         CALL DMATGEN(LUDA,LDA,N,LUDB)
         CALL DLUX(LUDA,LDA,N,IPVT,IERR)
         CALL DLUS(LUDA,LDA,N,IPVT,X,LUDB)
C         
  200 CONTINUE
      RETURN
      END

      SUBROUTINE DMATGEN(LUDA,LDA,N,LUDB)
      DOUBLE PRECISION LUDA(LDA,*),LUDB(*)
C
      INIT = 1325
      DO 30 J = 1,N
         DO 20 I = 1,N
            INIT = MOD(3125*INIT,65536)
            LUDA(I,J) = (INIT - 32768.0)/16384.0
   20    CONTINUE
   30 CONTINUE
      DO 35 I = 1,N
          LUDB(I) = 0.0
   35 CONTINUE
      DO 50 J = 1,N
         DO 40 I = 1,N
            LUDB(I) = LUDB(I) + LUDA(I,J)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END

      SUBROUTINE DLUS (LUDA,LDA,N,IPVT,X,LUDB)
      INTEGER IPVT(*)
      DOUBLE PRECISION LUDA(LDA,*), X(*), LUDB(*), XK
C
C   PURPOSE:
C     Solve the linear system ax = b given the LU factorization of A (as
C     computed in LU).
C
C   ADDITIONAL PARAMETERS (NOT PARAMETERS OF LU):
C
C     X double precision(N), solution of linear system
C
C     B double precision(N), right-hand-side of linear system
C
C ----------------------------------------------------------------------
C
      DO 10 K = 1, N
       X(K) = LUDB(K)
   10 CONTINUE
C
      DO 20 K = 1, N
       L = IPVT(K)
       XK = X(L)
       X(L) = X(K)
       X(K) = XK
   20 CONTINUE
C
      DO 40 K = 1, N
       XK = X(K)*LUDA(K,K)
       DO 30 I = K+1, N
          X(I) = X(I) - LUDA(I,K)*XK
   30    CONTINUE
       X(K) = XK
   40 CONTINUE
C
      DO 60 K = N, 1, -1
       XK = X(K)
       DO 50 I = 1, K-1
          X(I) = X(I) + LUDA(I,K)*XK          
   50    CONTINUE
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE DLU1(LUDA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      double precision LUDA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL DMXPY1(N-J+1,LUDA(J,J),J-1,LDA,LUDA(1,J),LUDA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUDA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUDA(I,J)) .GT. T) THEN
               T = ABS(LUDA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUDA(J,I)
            LUDA(J,I) = LUDA(K,I)
            LUDA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUDA(J,J) = 1/LUDA(J,J)
         CALL DXMPY1(N-J,LDA,LUDA(J,J+1),J-1,LDA,LUDA(J,1)
     1     ,LDA,LUDA(1,J+1))
         T = -LUDA(J,J)
         DO 30 I = J+1, N
            LUDA(J,I) = T*LUDA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE DLU2(LUDA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      double precision LUDA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL DMXPY2(N-J+1,LUDA(J,J),J-1,LDA,LUDA(1,J),LUDA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUDA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUDA(I,J)) .GT. T) THEN
               T = ABS(LUDA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUDA(J,I)
            LUDA(J,I) = LUDA(K,I)
            LUDA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUDA(J,J) = 1/LUDA(J,J)
         CALL DXMPY2(N-J,LDA,LUDA(J,J+1),J-1,LDA,LUDA(J,1)
     1     ,LDA,LUDA(1,J+1))
         T = -LUDA(J,J)
         DO 30 I = J+1, N
            LUDA(J,I) = T*LUDA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE DLU4(LUDA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      double precision LUDA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL DMXPY4(N-J+1,LUDA(J,J),J-1,LDA,LUDA(1,J),LUDA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUDA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUDA(I,J)) .GT. T) THEN
               T = ABS(LUDA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUDA(J,I)
            LUDA(J,I) = LUDA(K,I)
            LUDA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUDA(J,J) = 1/LUDA(J,J)
         CALL DXMPY4(N-J,LDA,LUDA(J,J+1),J-1,LDA,LUDA(J,1)
     1     ,LDA,LUDA(1,J+1))
         T = -LUDA(J,J)
         DO 30 I = J+1, N
            LUDA(J,I) = T*LUDA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE DLU8(LUDA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      double precision LUDA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL DMXPY8(N-J+1,LUDA(J,J),J-1,LDA,LUDA(1,J),LUDA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUDA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUDA(I,J)) .GT. T) THEN
               T = ABS(LUDA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUDA(J,I)
            LUDA(J,I) = LUDA(K,I)
            LUDA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUDA(J,J) = 1/LUDA(J,J)
         CALL DXMPY8(N-J,LDA,LUDA(J,J+1),J-1,LDA,LUDA(J,1)
     1     ,LDA,LUDA(1,J+1))
         T = -LUDA(J,J)
         DO 30 I = J+1, N
            LUDA(J,I) = T*LUDA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE DLUX(LUDA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      double precision LUDA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL DMXPYX(N-J+1,LUDA(J,J),J-1,LDA,LUDA(1,J),LUDA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUDA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUDA(I,J)) .GT. T) THEN
               T = ABS(LUDA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUDA(J,I)
            LUDA(J,I) = LUDA(K,I)
            LUDA(K,I) = T             
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUDA(J,J) = 1/LUDA(J,J)
         CALL DXMPYX(N-J,LDA,LUDA(J,J+1),J-1,LDA,LUDA(J,1)
     1     ,LDA,LUDA(1,J+1))
         T = -LUDA(J,J)
         DO 30 I = J+1, N
            LUDA(J,I) = T*LUDA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE DMXPY1 (N1,Y,N2,LDM,X,M)
C ** DMXPY01.F -- DMXPY UNROLLED TO A DEPTH OF 1
C
      INTEGER LDM,N1,N2
      double precision Y(*),X(*),M(LDM,*)
C     double precision Y(N1),X(N2),M(LDM,N2)
C
      DO 20 J = 1, N2
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE DMXPY2(N1,Y,N2,LDM,X,M)
C ** DMXPY02.F -- DMXPY UNROLLED TO A DEPTH OF 2
C
      INTEGER LDM,N1,N2
      double precision Y(*),X(*),M(LDM,*)
C     double precision Y(N1),X(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Main loop - groups of two vectors
C
      JMIN = J+2
      DO 30 J = JMIN, N2, 2
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE DMXPY4(N1,Y,N2,LDM,X,M)
C ** DMXPY04.F -- DMXPY UNROLLED TO A DEPTH OF 4
C
      INTEGER LDM,N1,N2
      double precision Y(*),X(*),M(LDM,*)
C     double precision Y(N1),X(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Main loop - groups of four vectors
C
      JMIN = J+4
      DO 40 J = JMIN, N2, 4
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE DMXPY8(N1,Y,N2,LDM,X,M)
C ** DMXPY08.F -- DMXPY UNROLLED TO A DEPTH OF 8
C
      INTEGER LDM,N1,N2
      double precision Y(*),X(*),M(LDM,*)
C     double precision Y(N1),X(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   30    CONTINUE
      ENDIF
C
C   Main loop - groups of eight vectors
C
      JMIN = J+8
      DO 50 J = JMIN, N2, 8
         DO 40 I = 1, N1
            Y(I) = ((((((( (Y(I))
     $             + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6))
     $             + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE DMXPYX(N1,Y,N2,LDM,X,M)
C ** DMXPY16.F -- DMXPY UNROLLED TO A DEPTH OF 16
C
      INTEGER LDM,N1,N2
      double precision Y(*),X(*),M(LDM,*)
C     double precision Y(N1),X(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + X(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + X(J-1)*M(I,J-1)) + X(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   30    CONTINUE
      ENDIF
C
C   Cleanup odd group of eight vectors
C
      J = MOD(N2,16)
      IF (J .GE. 8) THEN
         DO 40 I = 1, N1
            Y(I) = ((((((( (Y(I))
     $             + X(J-7)*M(I,J-7)) + X(J-6)*M(I,J-6))
     $             + X(J-5)*M(I,J-5)) + X(J-4)*M(I,J-4))
     $             + X(J-3)*M(I,J-3)) + X(J-2)*M(I,J-2))
     $             + X(J-1)*M(I,J-1)) + X(J)  *M(I,J)
   40    CONTINUE
      ENDIF
C
C   Main loop - groups of sixteen vectors
C
      JMIN = J+16
      DO 60 J = JMIN, N2, 16
         DO 50 I = 1, N1
            Y(I) = ((((((((((((((( (Y(I))
     $             + X(J-15)*M(I,J-15)) + X(J-14)*M(I,J-14))
     $             + X(J-13)*M(I,J-13)) + X(J-12)*M(I,J-12))
     $             + X(J-11)*M(I,J-11)) + X(J-10)*M(I,J-10))
     $             + X(J- 9)*M(I,J- 9)) + X(J- 8)*M(I,J- 8))
     $             + X(J- 7)*M(I,J- 7)) + X(J- 6)*M(I,J- 6))
     $             + X(J- 5)*M(I,J- 5)) + X(J- 4)*M(I,J- 4))
     $             + X(J- 3)*M(I,J- 3)) + X(J- 2)*M(I,J- 2))
     $             + X(J- 1)*M(I,J- 1)) + X(J)   *M(I,J)
   50    CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE DXMPY1(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** DXMPY01.F -- DXMPY UNROLLED TO A DEPTH OF 1
C
      INTEGER N1,LDY,N2,LDX,LDM
      double precision Y(LDY,*),X(LDX,*),M(LDM,*)
C     double precision Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
      DO 20 J = 1, N2
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE DXMPY2(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** DXMPY02.F -- DXMPY UNROLLED TO A DEPTH OF 2
C
      INTEGER N1,LDY,N2,LDX,LDM
      double precision Y(LDY,*),X(LDX,*),M(LDM,*)
C     double precision Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Main loop - groups of two vectors
C
      JMIN = J+2
      DO 30 J = JMIN, N2, 2
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $             + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE DXMPY4(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** DXMPY04.F -- DXMPY UNROLLED TO A DEPTH OF 4
C
      INTEGER N1,LDY,N2,LDX,LDM
      double precision Y(LDY,*),X(LDX,*),M(LDM,*)
C     double precision Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Main loop - groups of four vectors
C
      JMIN = J+4
      DO 40 J = JMIN, N2, 4
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)  *M(J,I)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE DXMPY8(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** DXMPY08.F -- DXMPY UNROLLED TO A DEPTH OF 8
C
      INTEGER N1,LDY,N2,LDX,LDM
      double precision Y(LDY,*),X(LDX,*),M(LDM,*)
C     double precision Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)  *M(J,I)
   30    CONTINUE
      ENDIF
C
C   Main loop - groups of eight vectors
C
      JMIN = J+8
      DO 50 J = JMIN, N2, 8
         DO 40 I = 1, N1
            Y(1,I) = ((((((( (Y(1,I))
     $               + X(1,J-7)*M(J-7,I)) + X(1,J-6)*M(J-6,I))
     $               + X(1,J-5)*M(J-5,I)) + X(1,J-4)*M(J-4,I))
     $               + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)  *M(J,I)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE DXMPYX(N1,LDY,Y,N2,LDX,X,LDM,M)
C ** DXMPY16.F -- DXMPY UNROLLED TO A DEPTH OF 16
C
      INTEGER N1,LDY,N2,LDX,LDM
      double precision Y(LDY,*),X(LDX,*),M(LDM,*)
C     double precision Y(LDY,N1),X(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + X(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)  *M(J,I)
   30    CONTINUE
      ENDIF
C
C   Cleanup odd group of eight vectors
C
      J = MOD(N2,16)
      IF (J .GE. 8) THEN
         DO 40 I = 1, N1
            Y(1,I) = ((((((( (Y(1,I))
     $               + X(1,J-7)*M(J-7,I)) + X(1,J-6)*M(J-6,I))
     $               + X(1,J-5)*M(J-5,I)) + X(1,J-4)*M(J-4,I))
     $               + X(1,J-3)*M(J-3,I)) + X(1,J-2)*M(J-2,I))
     $               + X(1,J-1)*M(J-1,I)) + X(1,J)  *M(J,I)
   40    CONTINUE
      ENDIF
C
C   Main loop - groups of sixteen vectors
C
      JMIN = J+16
      DO 60 J = JMIN, N2, 16
         DO 50 I = 1, N1
            Y(1,I) = ((((((((((((((( (Y(1,I))
     $               + X(1,J-15)*M(J-15,I)) + X(1,J-14)*M(J-14,I))
     $               + X(1,J-13)*M(J-13,I)) + X(1,J-12)*M(J-12,I))
     $               + X(1,J-11)*M(J-11,I)) + X(1,J-10)*M(J-10,I))
     $               + X(1,J- 9)*M(J- 9,I)) + X(1,J- 8)*M(J- 8,I))
     $               + X(1,J- 7)*M(J- 7,I)) + X(1,J- 6)*M(J- 6,I))
     $               + X(1,J- 5)*M(J- 5,I)) + X(1,J- 4)*M(J- 4,I))
     $               + X(1,J- 3)*M(J- 3,I)) + X(1,J- 2)*M(J- 2,I))
     $               + X(1,J- 1)*M(J- 1,I)) + X(1,J)   *M(J,I)
   50    CONTINUE
   60 CONTINUE
      RETURN
      END
C Caveat receptor.  (Jack) dongarra@anl-mcs.arpa, (Eric Grosse) research!ehg
C Careful! Anything free comes with no guarantee.
C *** from netlib, Mon Aug 22 11:25:55 CDT 1988 ***
      SUBROUTINE LUSD

      COMMON /LUSDC/ LUSA(301,300),LUSX(301),LUSB(301)
      COMMON /LF/ LOOPFACTOR
      INTEGER IPVT(300)
      REAL X
C
      LDA = 301
C
      DO 200 N = 80,120,20
C       N = LOOPFACTOR
         CALL SMATGEN(LUSA,LDA,N,LUSB)
         CALL SLU1(LUSA,LDA,N,IPVT,IERR)
         CALL SLUS(LUSA,LDA,N,IPVT,LUSX,LUSB)
C
         CALL SMATGEN(LUSA,LDA,N,LUSB)
         CALL SLU2(LUSA,LDA,N,IPVT,IERR)
         CALL SLUS(LUSA,LDA,N,IPVT,LUSX,LUSB)
C
         CALL SMATGEN(LUSA,LDA,N,LUSB)
         CALL SLU4(LUSA,LDA,N,IPVT,IERR)
         CALL SLUS(LUSA,LDA,N,IPVT,LUSX,LUSB)
C
         CALL SMATGEN(LUSA,LDA,N,LUSB)
         CALL SLU8(LUSA,LDA,N,IPVT,IERR)
         CALL SLUS(LUSA,LDA,N,IPVT,LUSX,LUSB)
C
         CALL SMATGEN(LUSA,LDA,N,LUSB)
         CALL SLUX(LUSA,LDA,N,IPVT,IERR)
         CALL SLUS(LUSA,LDA,N,IPVT,LUSX,LUSB)
C                            
  200 CONTINUE
      RETURN
      END

      SUBROUTINE SMATGEN(LUSA,LDA,N,LUSB)
      REAL LUSA(LDA,*),LUSB(*)
C
      INIT = 1325
      DO 30 J = 1,N
         DO 20 I = 1,N
            INIT = MOD(3125*INIT,65536)
            LUSA(I,J) = (INIT - 32768.0)/16384.0
   20    CONTINUE
   30 CONTINUE
      DO 35 I = 1,N
          LUSB(I) = 0.0
   35 CONTINUE
      DO 50 J = 1,N
         DO 40 I = 1,N
            LUSB(I) = LUSB(I) + LUSA(I,J)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END

      SUBROUTINE SLUS (LUSA,LDA,N,IPVT,LUSX,LUSB)
      INTEGER IPVT(*)
      REAL LUSA(LDA,*), LUSX(*), LUSB(*), XK
C
C   PURPOSE:
C     Solve the linear system Ax = b given the LU factorization of A (as
C     computed in LU).
C
C   ADDITIONAL PARAMETERS (NOT PARAMETERS OF LU):
C   
C     X REAL(N), solution of linear system
C
C     B REAL(N), right-hand-side of linear system
C
C ----------------------------------------------------------------------
C
      DO 10 K = 1, N
       LUSX(K) = LUSB(K)
   10 CONTINUE
C
      DO 20 K = 1, N
       L = IPVT(K)
       XK = LUSX(L)
       LUSX(L) = LUSX(K)
       LUSX(K) = XK
   20 CONTINUE
C
      DO 40 K = 1, N
       XK = LUSX(K)*LUSA(K,K)
       DO 30 I = K+1, N
          LUSX(I) = LUSX(I) - LUSA(I,K)*XK
   30    CONTINUE
       LUSX(K) = XK
   40 CONTINUE
C
      DO 60 K = N, 1, -1
       XK = LUSX(K)
       DO 50 I = 1, K-1
          LUSX(I) = LUSX(I) + LUSA(I,K)*XK
   50    CONTINUE
   60 CONTINUE
C
      RETURN
      END
      SUBROUTINE SLU1(LUSA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      REAL LUSA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL SMXPY1(N-J+1,LUSA(J,J),J-1,LDA,LUSA(1,J),LUSA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUSA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUSA(I,J)) .GT. T) THEN
               T = ABS(LUSA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUSA(J,I)
            LUSA(J,I) = LUSA(K,I)
            LUSA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUSA(J,J) = 1/LUSA(J,J)
         CALL SXMPY1(N-J,LDA,LUSA(J,J+1),J-1,LDA,LUSA(J,1)
     1      ,LDA,LUSA(1,J+1))
         T = -LUSA(J,J)
         DO 30 I = J+1, N
            LUSA(J,I) = T*LUSA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE SLU2(LUSA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      REAL LUSA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL SMXPY2(N-J+1,LUSA(J,J),J-1,LDA,LUSA(1,J),LUSA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUSA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUSA(I,J)) .GT. T) THEN
               T = ABS(LUSA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUSA(J,I)
            LUSA(J,I) = LUSA(K,I)
            LUSA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUSA(J,J) = 1/LUSA(J,J)
         CALL SXMPY2(N-J,LDA,LUSA(J,J+1),J-1,LDA,LUSA(J,1)
     1      ,LDA,LUSA(1,J+1))
         T = -LUSA(J,J)
         DO 30 I = J+1, N
            LUSA(J,I) = T*LUSA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE SLU4(LUSA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      REAL LUSA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL SMXPY4(N-J+1,LUSA(J,J),J-1,LDA,LUSA(1,J),LUSA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUSA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUSA(I,J)) .GT. T) THEN
               T = ABS(LUSA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUSA(J,I)
            LUSA(J,I) = LUSA(K,I)
            LUSA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUSA(J,J) = 1/LUSA(J,J)
         CALL SXMPY4(N-J,LDA,LUSA(J,J+1),J-1,LDA,LUSA(J,1)
     1     ,LDA,LUSA(1,J+1))
         T = -LUSA(J,J)
         DO 30 I = J+1, N
            LUSA(J,I) = T*LUSA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE SLU8(LUSA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      REAL LUSA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL SMXPY8(N-J+1,LUSA(J,J),J-1,LDA,LUSA(1,J),LUSA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUSA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUSA(I,J)) .GT. T) THEN
               T = ABS(LUSA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUSA(J,I)
            LUSA(J,I) = LUSA(K,I)
            LUSA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUSA(J,J) = 1/LUSA(J,J)
         CALL SXMPY8(N-J,LDA,LUSA(J,J+1),J-1,LDA,LUSA(J,1)
     1      ,LDA,LUSA(1,J+1))
         T = -LUSA(J,J)
         DO 30 I = J+1, N
            LUSA(J,I) = T*LUSA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE SLUX(LUSA,LDA,N,IPVT,INFO)
C ** LU.F -- LU DECOMPOSITION
C
C
      INTEGER LDA,N,IPVT(N),INFO
      REAL LUSA(LDA,*),t
C
      INFO = 0
      DO 40 J = 1, N
C
C       FORM J-TH COLUMN OF L
C
         CALL SMXPYX(N-J+1,LUSA(J,J),J-1,LDA,LUSA(1,J),LUSA(J,1))
C
C       SEARCH FOR PIVOT
C
         T = ABS(LUSA(J,J))
         K = J
         DO 10 I = J+1, N
            IF (ABS(LUSA(I,J)) .GT. T) THEN
               T = ABS(LUSA(I,J))
               K = I
            END IF
   10    CONTINUE
         IPVT(J) = K
C
C       TEST FOR ZERO PIVOT
C
         IF (T .EQ. 0) THEN
            INFO = J
            GO TO 50
         ENDIF
C
C       INTERCHANGE ROWS
C
         DO 20 I = 1, N
            T = LUSA(J,I)
            LUSA(J,I) = LUSA(K,I)
            LUSA(K,I) = T
   20    CONTINUE
C
C       FORM J-TH ROW OF U
C
         LUSA(J,J) = 1/LUSA(J,J)
         CALL SXMPYX(N-J,LDA,LUSA(J,J+1),J-1,LDA,LUSA(J,1)
     1      ,LDA,LUSA(1,J+1))
         T = -LUSA(J,J)
         DO 30 I = J+1, N
            LUSA(J,I) = T*LUSA(J,I)
   30    CONTINUE
   40 CONTINUE
   50 RETURN
      END
      SUBROUTINE SMXPY1 (N1,Y,N2,LDM,LUSX,M)
C ** SMXPY01.F -- SMXPY UNROLLED TO A DEPTH OF 1
C
      INTEGER LDM,N1,N2
      REAL Y(*),LUSX(*),M(LDM,*)
C     REAL Y(N1),LUSX(N2),M(LDM,N2)
C
      DO 20 J = 1, N2
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + LUSX(J)*M(I,J)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE SMXPY2(N1,Y,N2,LDM,LUSX,M)
C ** SMXPY02.F -- SMXPY UNROLLED TO A DEPTH OF 2
C
      INTEGER LDM,N1,N2
      REAL Y(*),LUSX(*),M(LDM,*)
C     REAL Y(N1),LUSX(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + LUSX(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Main loop - groups of two vectors
C
      JMIN = J+2
      DO 30 J = JMIN, N2, 2
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)*M(I,J)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE SMXPY4(N1,Y,N2,LDM,LUSX,M)
C ** SMXPY04.F -- SMXPY UNROLLED TO A DEPTH OF 4
C
      INTEGER LDM,N1,N2
      REAL Y(*),LUSX(*),M(LDM,*)
C     REAL Y(N1),LUSX(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + LUSX(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Main loop - groups of four vectors
C
      JMIN = J+4
      DO 40 J = JMIN, N2, 4
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + LUSX(J-3)*M(I,J-3)) + LUSX(J-2)*M(I,J-2))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)  *M(I,J)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE SMXPY8(N1,Y,N2,LDM,LUSX,M)
C ** SMXPY08.F -- SMXPY UNROLLED TO A DEPTH OF 8
C
      INTEGER LDM,N1,N2
      REAL Y(*),LUSX(*),M(LDM,*)
C     REAL Y(N1),LUSX(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + LUSX(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + LUSX(J-3)*M(I,J-3)) + LUSX(J-2)*M(I,J-2))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)  *M(I,J)
   30    CONTINUE
      ENDIF
C
C   Main loop - groups of eight vectors
C
      JMIN = J+8
      DO 50 J = JMIN, N2, 8
         DO 40 I = 1, N1
            Y(I) = ((((((( (Y(I))
     $             + LUSX(J-7)*M(I,J-7)) + LUSX(J-6)*M(I,J-6))
     $             + LUSX(J-5)*M(I,J-5)) + LUSX(J-4)*M(I,J-4))
     $             + LUSX(J-3)*M(I,J-3)) + LUSX(J-2)*M(I,J-2))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)  *M(I,J)
   40    CONTINUE                             
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SMXPYX(N1,Y,N2,LDM,LUSX,M)
C ** SMXPY16.F -- SMXPY UNROLLED TO A DEPTH OF 16
C                  
      INTEGER LDM,N1,N2
      REAL Y(*),LUSX(*),M(LDM,*)
C     REAL Y(N1),LUSX(N2),M(LDM,N2)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(I) = (Y(I)) + LUSX(J)*M(I,J)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(I) = ( (Y(I))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)*M(I,J)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(I) = ((( (Y(I))
     $             + LUSX(J-3)*M(I,J-3)) + LUSX(J-2)*M(I,J-2))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)  *M(I,J)
   30    CONTINUE
      ENDIF
C
C   Cleanup odd group of eight vectors
C
      J = MOD(N2,16)
      IF (J .GE. 8) THEN
         DO 40 I = 1, N1
            Y(I) = ((((((( (Y(I))
     $             + LUSX(J-7)*M(I,J-7)) + LUSX(J-6)*M(I,J-6))
     $             + LUSX(J-5)*M(I,J-5)) + LUSX(J-4)*M(I,J-4))
     $             + LUSX(J-3)*M(I,J-3)) + LUSX(J-2)*M(I,J-2))
     $             + LUSX(J-1)*M(I,J-1)) + LUSX(J)  *M(I,J)
   40    CONTINUE
      ENDIF
C
C   Main loop - groups of sixteen vectors
C
      JMIN = J+16
      DO 60 J = JMIN, N2, 16
         DO 50 I = 1, N1
            Y(I) = ((((((((((((((( (Y(I))
     $             + LUSX(J-15)*M(I,J-15)) + LUSX(J-14)*M(I,J-14))
     $             + LUSX(J-13)*M(I,J-13)) + LUSX(J-12)*M(I,J-12))
     $             + LUSX(J-11)*M(I,J-11)) + LUSX(J-10)*M(I,J-10))
     $             + LUSX(J- 9)*M(I,J- 9)) + LUSX(J- 8)*M(I,J- 8))
     $             + LUSX(J- 7)*M(I,J- 7)) + LUSX(J- 6)*M(I,J- 6))
     $             + LUSX(J- 5)*M(I,J- 5)) + LUSX(J- 4)*M(I,J- 4))
     $             + LUSX(J- 3)*M(I,J- 3)) + LUSX(J- 2)*M(I,J- 2))
     $             + LUSX(J- 1)*M(I,J- 1)) + LUSX(J)   *M(I,J)
   50    CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE SXMPY1(N1,LDY,Y,N2,LDX,LUSX,LDM,M)
C ** SXMPY01.F -- SXMPY UNROLLED TO A DEPTH OF 1
C
      INTEGER N1,LDY,N2,LDX,LDM
      REAL Y(LDY,*),LUSX(LDX,*),M(LDM,*)
C     REAL Y(LDY,N1),LUSX(LDX,N2),M(LDM,N1)
C
      DO 20 J = 1, N2
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + LUSX(1,J)*M(J,I)
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE SXMPY2(N1,LDY,Y,N2,LDX,LUSX,LDM,M)
C ** SXMPY02.F -- SXMPY UNROLLED TO A DEPTH OF 2
C
      INTEGER N1,LDY,N2,LDX,LDM
      REAL Y(LDY,*),LUSX(LDX,*),M(LDM,*)
C     REAL Y(LDY,N1),LUSX(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + LUSX(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Main loop - groups of two vectors
C
      JMIN = J+2
      DO 30 J = JMIN, N2, 2
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $             + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)*M(J,I)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
      SUBROUTINE SXMPY4(N1,LDY,Y,N2,LDX,LUSX,LDM,M)
C ** SXMPY04.F -- SXMPY UNROLLED TO A DEPTH OF 4
C
      INTEGER N1,LDY,N2,LDX,LDM
      REAL Y(LDY,*),LUSX(LDX,*),M(LDM,*)
C     REAL Y(LDY,N1),LUSX(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + LUSX(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Main loop - groups of four vectors
C
      JMIN = J+4
      DO 40 J = JMIN, N2, 4
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + LUSX(1,J-3)*M(J-3,I)) + LUSX(1,J-2)*M(J-2,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)  *M(J,I)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE SXMPY8(N1,LDY,Y,N2,LDX,LUSX,LDM,M)
C ** SXMPY08.F -- SXMPY UNROLLED TO A DEPTH OF 8
C
      INTEGER N1,LDY,N2,LDX,LDM
      REAL Y(LDY,*),LUSX(LDX,*),M(LDM,*)
C     REAL Y(LDY,N1),LUSX(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + LUSX(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + LUSX(1,J-3)*M(J-3,I)) + LUSX(1,J-2)*M(J-2,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)  *M(J,I)
   30    CONTINUE
      ENDIF
C
C   Main loop - groups of eight vectors
C
      JMIN = J+8
      DO 50 J = JMIN, N2, 8
         DO 40 I = 1, N1
            Y(1,I) = ((((((( (Y(1,I))
     $               + LUSX(1,J-7)*M(J-7,I)) + LUSX(1,J-6)*M(J-6,I))
     $               + LUSX(1,J-5)*M(J-5,I)) + LUSX(1,J-4)*M(J-4,I))
     $               + LUSX(1,J-3)*M(J-3,I)) + LUSX(1,J-2)*M(J-2,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)  *M(J,I)
   40    CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SXMPYX(N1,LDY,Y,N2,LDX,LUSX,LDM,M)
C ** SXMPY16.F -- SXMPY UNROLLED TO A DEPTH OF 16
C
      INTEGER N1,LDY,N2,LDX,LDM
      REAL Y(LDY,*),LUSX(LDX,*),M(LDM,*)
C     REAL Y(LDY,N1),LUSX(LDX,N2),M(LDM,N1)
C
C   Cleanup odd vector
C
      J = MOD(N2,2)
      IF (J .GE. 1) THEN
         DO 10 I = 1, N1
            Y(1,I) = (Y(1,I)) + LUSX(1,J)*M(J,I)
   10    CONTINUE
      ENDIF
C
C   Cleanup odd group of two vectors
C
      J = MOD(N2,4)
      IF (J .GE. 2) THEN
         DO 20 I = 1, N1
            Y(1,I) = ( (Y(1,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)*M(J,I)
   20    CONTINUE
      ENDIF
C
C   Cleanup odd group of four vectors
C
      J = MOD(N2,8)
      IF (J .GE. 4) THEN
         DO 30 I = 1, N1
            Y(1,I) = ((( (Y(1,I))
     $               + LUSX(1,J-3)*M(J-3,I)) + LUSX(1,J-2)*M(J-2,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)  *M(J,I)
   30    CONTINUE
      ENDIF
C
C   Cleanup odd group of eight vectors
C
      J = MOD(N2,16)
      IF (J .GE. 8) THEN
         DO 40 I = 1, N1
            Y(1,I) = ((((((( (Y(1,I))
     $               + LUSX(1,J-7)*M(J-7,I)) + LUSX(1,J-6)*M(J-6,I))
     $               + LUSX(1,J-5)*M(J-5,I)) + LUSX(1,J-4)*M(J-4,I))
     $               + LUSX(1,J-3)*M(J-3,I)) + LUSX(1,J-2)*M(J-2,I))
     $               + LUSX(1,J-1)*M(J-1,I)) + LUSX(1,J)  *M(J,I)
   40    CONTINUE
      ENDIF
C
C   Main loop - groups of sixteen vectors
C
      JMIN = J+16
      DO 60 J = JMIN, N2, 16
         DO 50 I = 1, N1
            Y(1,I) = ((((((((((((((( (Y(1,I))
     $               + LUSX(1,J-15)*M(J-15,I)) + LUSX(1,J-14)*M(J-14,I))
     $               + LUSX(1,J-13)*M(J-13,I)) + LUSX(1,J-12)*M(J-12,I))
     $               + LUSX(1,J-11)*M(J-11,I)) + LUSX(1,J-10)*M(J-10,I))
     $               + LUSX(1,J- 9)*M(J- 9,I)) + LUSX(1,J- 8)*M(J- 8,I))
     $               + LUSX(1,J- 7)*M(J- 7,I)) + LUSX(1,J- 6)*M(J- 6,I))
     $               + LUSX(1,J- 5)*M(J- 5,I)) + LUSX(1,J- 4)*M(J- 4,I))
     $               + LUSX(1,J- 3)*M(J- 3,I)) + LUSX(1,J- 2)*M(J- 2,I))
     $               + LUSX(1,J- 1)*M(J- 1,I)) + LUSX(1,J)   *M(J,I)
   50    CONTINUE
   60 CONTINUE
      RETURN
      END

C Caveat receptor.  (Jack) dongarra@anl-mcs.arpa, (Eric Grosse) research!ehg
C Careful! Anything free comes with no guarantee.
C *** from netlib, Mon Aug 22 11:25:36 CDT 1988 ***

        SUBROUTINE NEFF
C
        REAL    ANU,CY,CZ,D,R,RI,SY,X,Y
        INTEGER J,K,NMAX1,NMAX2
C
        COMMON /COMP/ CZ,ANU,R,RI
        COMMON /LF/ LOOPFACTOR
C
        NMAX1  = 100
        NMAX2  = 400
        CZ  = .005
        ANU = .287
C        PI  = 3.1416
C        C   = 6.383 / SQRT(ANU)
        D   = 1. / ( (4.*CZ)**3*SQRT(ANU) )
        RI = 0.          
C       N = NMAX1 + 1
C
        DO 100 J=1,(NMAX1*LOOPFACTOR)
C           X = H1*(J-1)  ... delete because h1 not init'd
           X = J-1              
           Y = D*X  
           CALL FRENEL (Y,NMAX2,SY,CY)
           R = EXP(-.25*X*X)*(10.-2.35*X*X+.1249*X**4)*(CY-SY)
           K = 1-MOD(J,2)
           R = (2+2*K)*R
           RI = RI+R
100        CONTINUE
        RETURN
        END
C
C*********************************************************************
C
        SUBROUTINE FRENEL (Y,NMAX2,SY,CY)
C
        REAL ANU,CZ,PI,R,RI
        REAL C1,S1,CF,SF,CI,SI,CY,SY
        REAL ELIM,ENMAX2,H2,YY,Z,ZZ
        INTEGER I,J,M
C
        COMMON /COMP/ CZ,ANU,R,RI
        DATA PI /3.1415926/
C
        ELIM = SQRT(Y)
        ENMAX2 = NMAX2
        H2 = ELIM/ENMAX2
        S1 = 0.0
        C1 = 1.0
        YY = AMOD(Y,2.*PI)
        SF = SIN(YY)
        CF = COS(YY)
        M  = NMAX2 + 1
        SY = 0.
        CY = 0.
        DO 100 I=1,M
           Z = H2*(I-1)*4.*CZ
           ZZ = Z*Z
           ZZ = AMOD(ZZ,2*PI)
           SI = SIN(ZZ)
           CI = COS(ZZ)
           J = 1-MOD(I,2)
           SI = (2+2*J)*SI
           CI = (2+2*J)*CI
           SY = SY+SI
           CY = CY+CI
100        CONTINUE
        SY = .3333*H2*(SY-S1-SF)
        CY = .3333*H2*(CY-C1-CF)
C                                     
        RETURN
        END  
        SUBROUTINE PHILCO
C
C        DOUBLE PRECISION BENCHMARK THAT ASSUMES 48 BITS OF
C        PRECISION IN THE MANTISSA.
C
        DOUBLE PRECISION CCM2,INV(12,12),IP(12,12),TEMP7
        COMMON /PHILC/ CCM2(12,12),TEMP7(12,12)
        COMMON /LF/ LOOPFACTOR
        REAL LABEL(12)
        INTEGER I,IDIM,IND,J,K,L,M,N,NCARG,NRARG
C------This data moved to the INIT subroutine
C        DATA TEMP7
C     1         /18318.1D0,-2374.5D0,14.9D0,10862.5D0,-623.6D0,327.6
C     1 D0,24421.7D0,5084.9D0,758.2D0,22163.9D0,-11529.6D0,-1468.6
C     1 D0,-2374.5D0,392.1D0,-13.8D0,-1348.9D0,180.8D0,-7.6
C     1 D0,-3106.7D0,-528.0D0,-164.8D0,-2732.6D0,1547.0D0,123.6
C     1 D0,14.9D0,-13.8D0,197.5D0,22.7D0,3.7D0,3.1
C     1 D0,47.6D0,15.2D0,-3.1D0,47.1D0,-18.9D0,-7.2
C     1 D0,10862.5D0,-1348.9D0,22.7D0,6718.4D0,-315.D0,221.9
C     1 D0,14894.3D0,3188.9D0,402.2D0,13588.9D0,-6963.1D0,-928.1
C     1 D0,-623.6D0,180.8D0,3.7D0,-315.0D0,171.5D0,54.4
C     1 D0,-771.5D0,-12.3D0,-138.4D0,-612.5D0,470.3D0,-60.3
C     1 D0,327.6D0,-7.6D0,3.1D0,221.9D0,54.4D0,60.
C     1 D0,476.D0,153.9D0,-31.6D0,467.5D0,-186.9D0,-72.1D0,
C     1    24421.7D0,-3106.7D0,47.6D0,14894.3D0,-771.5D0,476.
C     1 D0,33517.8D0,7071.6D0,1003.4D0,30401.1D0,-15707.3D0,-2059.9
C     1 D0,5084.9D0,-528.D0,15.2D0,3188.9D0,-12.3D0,153.9
C     1 D0,7071.6D0,1693.D0,62.2D0,6549.7D0,-3202.3D0,-527.7
C     1 D0,758.2D0,-164.8D0,-3.1D0,402.2D0,-138.4D0,-31.6
C     1 D0,1003.4D0,62.2D0,233.2D0,847.1D0,-542.3D0,26.4
C     1 D0,22163.9D0,-2732.6D0,47.1D0,13588.9D0,-612.5D0,467.5
C     1 D0,30401.1D0,6549.7D0,847.1D0,27846.6D0,-14237.4D0,-1944.2
C     1 D0,-11529.6D0,1547.D0,-18.9D0,-6963.1D0,470.3D0,-186.9
C     1 D0,-15707.3D0,-3202.3D0,-542.3D0,-14237.4D0,7498.1D0,864.4
C     1 D0,-1468.6D0,123.6D0,-7.2D0,-928.1D0,-60.3D0,-72.1
C     1 D0,-2059.9D0,-527.7D0,26.4D0,-1944.2D0,864.4D0,301.3D0/
C
C------As well as this initialiaztion loop
C        DO 20 I=1,12
C           DO 10 J=1,12
C              CCM2(I,J) = TEMP7(J,I)
C10           CONTINUE
C20        CONTINUE
C
        DO 50 K=1,(100*LOOPFACTOR)
           DO 40 I=1,12
              DO 30 J=1,12
                 INV(I,J)=CCM2(I,J)
30              CONTINUE
40           CONTINUE
           NRARG = 12
           NCARG = 12
           IDIM  = 12  
           CALL INVERT (INV,NRARG,NCARG,IDIM,LABEL)
           IND = 0
           L   = 12
           M   = 12
           N   = 12
           CALL MTMPY (IND,CCM2,INV,IP,L,M,N)
50        CONTINUE
        RETURN
C
        END
C
C*******************************************************************
C
        SUBROUTINE INVERT (A,NRARG,NCARG,IDIM,LABEL)
C
C------------------ MATRIX INVERSION --------------------------
C
        INTEGER I,IBIG,IDIM,J1,J2,J3
        INTEGER N1,NC,NCARG,NR,NRARG
        REAL    LABEL(NRARG)
        DOUBLE PRECISION A(IDIM,NCARG),TEMP
C
1        NR=NRARG
        NC=NCARG
        DO 21 J1=1,NR
           LABEL(J1)=J1
21        CONTINUE
C
        DO 291 J1=1,NR
C
C...... FIND REMAINING ROW CONTAINING LARGEST ABSOLUTE
C...... VALUE IN PIVOTAL COLUMN
C
101        TEMP=0.0D+00
        DO 121 J2=J1,NR
           IF(DABS(A(J2,J1)).LT.TEMP) GO TO 121
           TEMP=DABS(A(J2,J1))
           IBIG=J2
121        CONTINUE
        IF(IBIG.EQ.J1)GO TO 201
C
C...... REARRANGE ROWS TO PLACE LARGEST ABSOLUTE VALUE IN PIVOT POSITION
C
        DO 141 J2=1,NC
           TEMP=A(J1,J2)
           A(J1,J2)=A(IBIG,J2)
           A(IBIG,J2)=TEMP
141        CONTINUE
        I=LABEL(J1)
        LABEL(J1)=LABEL(IBIG)
        LABEL(IBIG)=I
C
C...... COMPUTE COEFFICIENTS IN PIVOTAL ROW
C
201        TEMP=A(J1,J1)
        A(J1,J1)=1.0D+00
        DO 221 J2=1,NC
           A(J1,J2)=A(J1,J2)/TEMP
221        CONTINUE
C
C...... COMPUTE COEFFICIENTS IN OTHER ROWS
C
        DO 281 J2=1,NR
           IF(J2.EQ.J1) GO TO 281
           TEMP=A(J2,J1)
           A(J2,J1)=0.0D+00
           DO 241 J3=1,NC
              A(J2,J3)=A(J2,J3)-TEMP*A(J1,J3)
241           CONTINUE
281        CONTINUE
291        CONTINUE
C
C...... INTERCHANGE COLUMNS ACCORDING TO INTERCHANGES OF
C...... ROWS OF ORIGINAL MATRIX
C
301        N1=NR-1
        DO 391 J1=1,N1
           DO 321 J2=J1,NR
              IF(LABEL(J2).NE.J1) GO TO 321
              IF(J2.EQ.J1) GO TO 391
              GO TO 341
321           CONTINUE
341           CONTINUE
           DO 361 J3=1,NR
              TEMP=A(J3,J1)
              A(J3,J1)=A(J3,J2)
              A(J3,J2)=TEMP
361           CONTINUE
           LABEL(J2)=LABEL(J1)
391        CONTINUE
5001        RETURN
        END
C
C**************************************************************
C
        SUBROUTINE MTMPY(IND,A,B,C,LARG,MARG,N)
C
C----------------- MATRIX MULTIPLICATION ----------------------
C
        DOUBLE PRECISION A(12,12),B(12,12),C(12,12),TEMP
        INTEGER I,IND,J1,J2,J3,L,LARG,M,MARG,N
C
1        L=IABS(LARG)
        M=IABS(MARG)
        I=IND+1
        GO TO (101,201,301,401),I
C
101        DO 121 J1=1,L
        DO 121 J2=1,N
        C(J1,J2)=0.0D0
        DO 121 J3=1,M
        IF(LARG)102,5001,110
102        IF(MARG)103,5001,105
103        TEMP=A(J3,J1)*B(J2,J3)
        GO TO 121
105         TEMP=A(J3,J1)*B(J3,J2)
        GO TO 121
110        IF(MARG)111,5001,115
111        TEMP=A(J1,J3)*B(J2,J3)
        GO TO 121
115        TEMP=A(J1,J3)*B(J3,J2)
121        C(J1,J2)=C(J1,J2)+TEMP
        GO TO 5001
C
201        DO 221 J1=1,L
221        C(J1,1)=A(J1,1)*B(J1,1)
        GO TO 5001
C
301        DO 321 J1=1,L
        DO 321 J2=1,M
        IF(MARG)310,5001,315
310        TEMP=A(J1,1)*B(J2,J1)
        GO TO 321
315        TEMP=A(J1,1)*B(J1,J2)
321        C(J1,J2)=TEMP
        GO TO 5001
C
401        DO 421 J1=1,L
        DO 421 J2=1,M
        IF(LARG)410,5001,415
410        TEMP=A(J2,J1)*B(J2,1)
        GO TO 421
415        TEMP=A(J1,J2)*B(J2,1)
421        C(J1,J2)=TEMP
5001        RETURN
        END
        SUBROUTINE PRIME
C
        INTEGER I,IM,IPRIME,J,K
        COMMON /PRIM/ IPRIME(1500)
        COMMON /LF/ LOOPFACTOR
C
        K = 2
        IPRIME(1) = 2
        DO 4 LP=1,(1 * LOOPFACTOR)
          DO 3 I=2,1500
2           K  = K + 1
            IM = I - 1
            DO 1 J = 1,IM
              IF ( (K/IPRIME(J))*IPRIME(J) .EQ. K ) GO TO 2
1           CONTINUE
            IPRIME(I) = K
3         CONTINUE
4       CONTINUE
        RETURN
        END
        SUBROUTINE RR2
C
C...... STEADY STATE TBA SIMULATION
C
        REAL VALOUT(5)
        COMMON /LF/ LOOPFACTOR
        DATA VALOUT / 5*1. /
C
        NTIMES = 100
C
        DO 200 LOOP = 1,(NTIMES * LOOPFACTOR)
          DO 120 I = 1,22    
            CALL ONE (VALOUT(1),VALOUT(2),VALOUT(3))
120       CONTINUE    
          DO 140 I = 1,10    
            CALL FOUR (VALOUT(4),VALOUT(5))  
140       CONTINUE    
200     CONTINUE 
C
        RETURN
        END
C
C************************************************************
C
        SUBROUTINE ONE (A,B,C) 
C
        REAL A,B,C,D
C
        CALL TWO (C,A,D,B)     
        A = B / SQRT(C*D) + 0.488596  
        CALL TWO (A,D,B,C)     
        B = SQRT((A+C)/D)      
        C = A*B - D     
        CALL TWO (B,C,A,D)     
C
        RETURN
        END      
C
C************************************************************
C
        SUBROUTINE TWO (C,A,D,B)      
C
        REAL A,B,C,D
C
        A = 0.19 
        B = 1.1  
        CALL THREE (A,B)
        IF (A) 1,1,2    
1       C = A    
        CALL THREE (D,C)
        IF (D) 3,6,4    
2       D = (0.5*(A+1.0))**((A-1.0)/A)
        CALL THREE (D,B)
        GO TO 1  
3       D = C    
        CALL THREE (C,B)
        IF (C) 2,5,5    
4       B = A    
        CALL THREE (D,B)
        GO TO 3  
5       C = B    
        IF (A*D) 1,3,6  
6       RETURN   
        END      
C
C************************************************************
C
        SUBROUTINE THREE (A,B) 
C
        REAL A,B
C
        A = A + B*B     
C
        RETURN
        END 
C
C************************************************************
C
        SUBROUTINE FOUR (D,E)  
C
        REAL A(4),D,E
C
        DATA A(1),A(2),A(3),A(4)/40.0,0.4,1.4,20.0/ 
C
        DO 200 I = 1,60 
          D = I    
          IF (D-80.0) 100,100,200
100       E = (((D+A(1))*A(2)+D)*A(3)+D) / A(4)
200     CONTINUE 
C
        RETURN
        END      
        SUBROUTINE RR3
C
C...... RIGS CALIBRATION
C
        INTEGER IDRO(562), IPSN(562)
        REAL    CVAL(562), CALDAT
C
        COMMON /LF/ LOOPFACTOR
        COMMON /CALIB/INDEX(300),CALDAT(1000)
C
        DATA INDEX / 210,  11,  32,   5,   3,   8,  70,  10,   0,   8,
     1                70,  10,   0,   8,  20,   0,   0,   0,   0,   0,
     2                 0,   0,   0,   0,   0,   0,   0,   0, -70, -15,
     3                 0,  15,  70, -70, -15,   0,  15,  70,-700,-150,
     4                 0, 150, 700,-700,-150,   0, 150, 700,-700,-150,
     5                 0, 150, 700,-700,-150,   0, 150, 700,-700,-150,
     6                 0, 150, 700,-700,-150,   0, 150, 700,-700,-150,
     7                 0, 700,2000,-700,-150,   0, 700,2000,-700,-150,
     8                 0, 700,3000,   0,   0,   0,   0,   0,   0,   0,
     9                 0,   0,   0,   0,   0,   0,   0,   0,   3,   0,
     1                  24*11, 76*0, 24*1, 76*0 /
C
        DATA CALDAT /-100., -80., -60., -40., -20.,   0.,  20.,  40.,
     1                 60.,  80., 100.,-100., -80., -60., -40., -20.,
     2                  0.,  20.,  40.,  60.,  80., 100., 978*0./
C
        DATA IDRO /
     1       0,  50,   2,   3,   4,   5,   6,   7,   8,   9, 80*10,
     2       0,  50,   2,   3,   4,   5,   6,   7,   8,   9, 80*10,
     3       0,  50,   2,   3,   4,   5,   6,   7,   8,   9, 20*10,
     4     -70, -15,   0,  15,  70,   27*100,
     5     -70, -15,   0,  15,  70,   27*200,
     6    -700,-150,   0, 150, 700,   27*300,
     7    -700,-150,   0, 150, 700,   27*400,
     8    -700,-150,   0, 150, 700,   27*500,
     9    -700,-150,   0, 150, 700,   27*600,
     1    -700,-150,   0, 150, 700,   27*700,
     2    -700,-150,   0, 150, 700,   27*800,
     3    -700,-150,   0, 700,2000,   27*900,
     4    -700,-150,   0, 700,2000,   27*1000,
     5    -700,-150,   0, 700,3000,   27*1100/
C
        NDROS  = 562
        NTIMES = 100
C
        DO 50 I = 1,NDROS
           IPSN(I) = I
50      CONTINUE
C
        DO 200 LOOP = 1,(NTIMES*LOOPFACTOR)
          DO 120 I = 1,NDROS
            CVAL(I) = 0.0
120       CONTINUE
          CALL RIGCAL (IDRO(1),IPSN(1),CVAL(1))
200     CONTINUE
C
        RETURN
        END
C
C*****************************************************************************
C
C                                                           
        SUBROUTINE RIGCAL (IDROS,IPSNS,CVAL)
C
C...... RIGS CALIBRATION ROUTINE
C
        INTEGER IDROS(1), IPSNS(1), NNPS, NSVS, NPORTS, NREFPS
        INTEGER NSCS, NCHNLS, ISPARE, ISVREF, IABS, NXTAB, IXTAB
        INTEGER MAXD, MIND, K, LMISC, ISC, ITYP, ITCTYP, NCHS
        INTEGER I, J, ID, M, NX, IX, IY, MEASP, ISV, IREF
        REAL    CVAL(1) , XTAB(5) , YTAB(5), CALDAT, REFV
        REAL    D0, D1, S, C, X, Y, ADD
C
        COMMON /CALIB/NNPS,NSVS,NPORTS,NREFPS,NSCS,NCHNLS(4,4),
     1                ISPARE(7),ISVREF(5,12),IABS(12),
     2                NXTAB(100),IXTAB(100),
     3                CALDAT(1000)
C
        DATA  MAXD,MIND/11999,-11999/
        DATA  REFV/50.0/
C
C...... CONVERT NON-PRESSURE DATA FOR EACH SCANNER
C
        K = 0
        LMISC = 0
        DO 290 ISC = 1,NSCS
C........ GET DRO TO MV CONVERSION
          K = K+2
          D0 = IDROS(K-1)
          D1 = IDROS(K)
          S = REFV/(D1-D0)
          C = -S*D0
C........ CONVERT EACH PARAMETER TYPE
          DO 280 ITYP = 1,4
            ITCTYP = ITYP-1
            NCHS = NCHNLS(ITYP,ISC)
            IF (NCHS .LE. 0) GO TO 280
C.......... CONVERT EACH VALUE OF THIS TYPE
            DO 270 I = 1,NCHS
              K = K+1
              J = IPSNS(K)
              IF (J.LE.0) GO TO 270
              ID = IDROS(K)
              IF (ID.GE.MAXD) GO TO 270
              IF (ID.LE.MIND) GO TO 270
              X = ID
              X = X*S+C
              IF(ITCTYP.GT.0) GO TO 260
              M = LMISC+I
              NX = NXTAB(M)
              IX = IXTAB(M)
              IY = IX+NX
              CALL Q966 (CALDAT(IX),CALDAT(IY),NX,X,CVAL(J))
              GO TO 270
260           CALL TCLU (ITCTYP,X,CVAL(J))
270         CONTINUE
280       CONTINUE
          LMISC = LMISC + NCHNLS(1,ISC)
290     CONTINUE
C
C...... CONVERT PRESSURE DATA FOR EACH SCANIVALVE
C
        MEASP = NPORTS-NREFPS
        DO 490 ISV = 1,NSVS
          K = NPORTS*(ISV-1) + NNPS
C........ GET PRESSURE TO BE ADDED IF ABSOLUTE OUTPUT
          ADD = 0.0
          J = IABS(ISV)
          IF (J.LE.0) GO TO 430
          J = IPSNS(J)
          IF (J.LE.0) GO TO 490
          ADD = CVAL(J)
          IF (ADD.LE.0.) GO TO 490
430     CONTINUE
C...... SET UP CAL TABLE
          DO 440 IREF = 1,NREFPS
            K = K+1
            ID = IDROS(K)
            IF (ID.GE.MAXD) GO TO 490
            IF (ID.LE.MIND) GO TO 490
            XTAB(IREF) = ID
            Y = ISVREF(IREF,ISV)
            YTAB(IREF) = Y+ADD
440       CONTINUE                                          
C ** CALCULATE PRESSURES
          DO 480 I = 1,MEASP
            K = K+1
            J = IPSNS(K)
            IF (J.LE.0) GO TO 480
            ID = IDROS(K)
            IF (ID.GE.MAXD) GO TO 480
            IF (ID.LE.MIND) GO TO 480
            X = ID
            CALL Q966 (XTAB(1),YTAB(1),NREFPS,X,CVAL(J))
480       CONTINUE
490     CONTINUE
C
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE TCLU (ITYP,XIN,Y)
C
        INTEGER IXTAB(3), LXTAB(3), ITYP, I, J, IX, LX
        REAL    XTAB(75), YTAB(75), XIN, Y, X
C
        DATA  IXTAB/ 2,27,52/
        DATA  LXTAB/25,50,75/
        DATA  XTAB/12*0.0,13*50.0,12*0.0,13*50.0,12*0.0,13*50.0/
        DATA  YTAB/12*0.0,13*50.0,12*0.0,13*100.0,12*0.0,13*150.0/
C
        IF (ITYP.LE.0) GO TO 90
        IX = IXTAB(ITYP)
        LX = LXTAB(ITYP)
        X = XIN
        DO 20 I = IX,LX
          IF (X.LE.XTAB(I)) GO TO 30
20      CONTINUE
        I = LX
30      J = I-1
        Y = YTAB(J) + (X-XTAB(J))*(YTAB(I)-YTAB(J)) / (XTAB(I)-XTAB(J))
        RETURN
90      Y = 0.0
C
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE Q966 (XTAB,YTAB,NX,X,Y)
C
        INTEGER NX, I, J
        REAL    XTAB(1), YTAB(1), X, Y
C
        DO 10 I = 2,NX
          IF (XTAB(I).GE.X) GO TO 20
10      CONTINUE
        I = NX
20      J = I-1
        Y = YTAB(J) + (X-XTAB(J))*(YTAB(I)-YTAB(J)) / (XTAB(I)-XTAB(J))
C
        RETURN
        END
        SUBROUTINE RR4
C
C       RIGS CALCULATIONS (BLADE COOLING)
C
        REAL CURVE1(20), CURVE2(20), VALOUT(34)
C
        COMMON /LF/ LOOPFACTOR
        COMMON /K564CI/ VALIN(72)
        COMMON /K564CO/ XMF,AFR,TID,RATIOS(6),TG,TC(3),CE(7,3)
        EQUIVALENCE (VALOUT(1),XMF)
C
        DATA  NOUT/34/ , NX1/10/ , NX2/10/
        DATA  CURVE1/ 20., 40., 60., 80.,100.,120.,140.,160.,180.,200.,
     1                      20., 40., 60., 80.,100.,120.,140.,160.,
     1                     180., 200./
        DATA  CURVE2/ 20., 40., 60., 80., 100., 120., 140., 160., 180.,
     1               200., 20., 40., 60., 80.,100.,120.,140.,160.,
     1               180., 200./
        DATA  VALIN /
     1     222.019 ,    0.0  ,    0.0,    0.0  ,  199.584,  206.934,
     2     160.000 ,    0.0  ,    0.0,    0.0  ,    0.0  ,    0.0  ,
     3       0.0   ,    0.0  ,    0.0,    4.392,   12.5  ,   12.0  ,
     4      12.9   ,   12.2  ,    0.0,    0.0  ,    0.0  ,    0.0  ,  
     5       0.2674,  685.039,    0.0,  409.929,  667.086,  679.000,
     6       0.0   ,    0.0  ,    0.0,  678.000,    0.0  ,    0.0  ,
     7       0.0   ,    0.0  ,    0.0,    0.0  ,    0.0  ,    0.0  ,
     8       0.0   ,    0.0  ,    0.0,    0.0  ,  970.0  ,  927.0  ,
     9     941.0   ,  930.0  ,  936.0, 1016.0  , 1109.0  ,    0.0  ,
     1       0.0   ,  941.0  ,  919.0,  946.0  , 1005.0  ,  776.0  ,
     2     978.0   ,  931.0  ,  946.0,    0.0  ,  944.0  , 1015.0  ,
     3    1166.0   , 1176.0  , 1175.0, 1161.0  , 1157.0  , 1166.0  /
C
        NTIMES = 1000

        DO 200 LOOP = 1,(NTIMES*LOOPFACTOR)
          DO 120 I = 1,NOUT
            VALOUT(I) = 0.0
120       CONTINUE
          CALL CALCS (CURVE1,NX1,CURVE2,NX2)
200     CONTINUE
C
        RETURN
        END
C
C***************************************************************************
C
        SUBROUTINE CALCS (CURVE1,NX1,CURVE2,NX2)
C
C       BLADE COOLING RIGS CALCULATION ROUTINE
C
        INTEGER NX1, NX2
        INTEGER I,J,K,N,NTT,NTU,NT14
        REAL CURVE1(1), CURVE2(1), T14M(5)
C
        COMMON /K564CI/ P1,P8,P9,P10,P12,P14,P15,P16,P20,P13,SPARE1(5),
     1                DP21,DP24(5),SPARE2(3),
     2                FF,T1,T8,T9,T9A,T14(5),T15,T15A,
     3                TR(5),T20,T20A,SPARE3(2),
     4                TCL(21),TUC(6)
        COMMON /K564CO/ XMF,AFR,TID,
     1                P14P12,P14P15,P14STG,P12P16,P14P16,PSTFDS,
     2                TG,TC(3),CE(21)
        DATA  AR/0.4225/
        DATA  CV/10370.0/
        DATA  ZERO/0.0/
        DATA  NT14/5/
        DATA  NTU,NTT/6,27/
        IF (P1.EQ.ZERO.OR.DP21.EQ.ZERO.OR.T1.EQ.ZERO.OR.FF.EQ.ZERO)
     1                                                        GO TO 900
C
C...... APPROXIMATE MASS FLOW AND AIR/FUEL RATIO
C
        DN = 3.25 + 0.00005863*(T1-293.1)
        RHO = 1.5*P1/T1
        TT1 = (T1-1125.0)/875.0
        CP = (((0.012773*TT1-0.007016)*TT1-0.021413)*TT1+0.037079)*TT1
     1                +0.277988
        CPX = 14.58596*CP
        GAMMA = CPX/(CPX-1.0)
        PR = 1.0-(DP21/P1)
        PR2G = PR**(2.0/GAMMA)
        PRGG1 = PR**((GAMMA-1.0)/GAMMA)
        XMF = DN*DN*0.49949 * SQRT(GAMMA*PR2G*(1.0-PRGG1)*RHO*DP21 /
     1                ((GAMMA-1.0)*(1.0-AR*AR*PR2G)*(1.0-PR)))
        AFR = XMF/FF
C
C...... FLAME TEMPERATURE
C
        K = 1
        IF (T9.GT.0 .AND. T9.LT.T1) K=2
        IF (T9A.GT.0 .AND. T9A.LT.T1) K = K+2
        GO TO (291,292,293,270),K
C
270     IF (T9-T9A-50.0) 280,280,292
280     IF (T9A-T9-50.0) 294,294,293
C
291     AIT = T1
        GO TO 300
C
292     AIT = T9
        GO TO 300
C
293     AIT = T9A
        GO TO 300
C
294     AIT = (T9+T9A)*0.5
C
300     FAR = FF/XMF
        CALL Q806 (FAR,AIT,CV,TID)
C
C...... PRESSURE RATIOS
C
        IF (P12.EQ.ZERO)  GO TO 400
        CALL MEANNZ (DP24(1),5,AVE)
        IF (AVE.EQ.0.0)  GO TO 400
        IF (P12 .GT. P14) AVE = -AVE
320     P14C = P12+AVE
        P14P12 = P14C / P12
        IF (P15.NE.ZERO)  P14P15 = P14C / P15
        IF (NX1.LE.0) GO TO 340
        CALL QQ966 (CURVE1(1),CURVE1(NX1+1),NX1,P12,PSTG)
        P14STG = P14C / PSTG
340     IF (P16.EQ.ZERO)  GO TO 400
        P12P16 = P12 / P16
        P14P16 = P14C/ P16
        IF (NX1.LE.0 .OR. NX2.LE.0) GO TO 400
        CALL QQ966 (CURVE2(1),CURVE2(NX2+1),NX2,P16,PFDS)
        PSTFDS = PSTG / PFDS
C
C...... COOLING EFFECTIVENESS
C
400     CALL MEANNZ (T14(1),NT14,TAVE)
        IF (TAVE.LE.0.0)  GO TO 600
        CALL MEANNZ (TCL(1),NTT,TMEAN)
        IF (TMEAN.LT.TAVE)  GO TO 600
C
        TR4 = TR(4)
        TMIN = TR4-70.0
        IF (TMIN .LE. 0.) GO TO 420
        TMAX = TR4
        TG = TR4-50.0
        GO TO 430
420     TMIN = TID-150.0
        IF (TMIN.LE.0.0)  GO TO 600
        TMAX = TID
        TG = TID
430     N = 0
        AV = 0.0
        DO 460 I = 1,NTU
          T = TUC(I)
          IF (T.LE.TMIN .OR. T.GE.TMAX) GO TO 460
          N = N+1
          AV = AV+T
460     CONTINUE
C
        IF (N.LE.0) GO TO 480
        AN = N
        TG = AV/AN
480     CONTINUE
C
        DO 490 I = 1,NT14
          T = T14(I)
          IF (T.LE.0.0)  T=TAVE
          T14M(I) = T
490     CONTINUE
C
        N = 0
        DO 540 I = 1,3
          TC(I) = (T14M(I)+T14M(I+1))*0.5
          D = TG - TC(I)
          DO 530 J = 1,7
            N = N+1
            IF (D.EQ.0) GO TO 530
            T = TCL(N)
            IF (T.NE.0) CE(N) = (TG-T)/D
530       CONTINUE
540     CONTINUE
C
600     CONTINUE
900     RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE MEANNZ (DATA,NDATA,AVGE)
C
        REAL DATA(1), AVGE, D, A, AN
C
        A = 0.0
        N = 0
        DO 120 I = 1,NDATA
          D = DATA(I)
          IF (D.EQ.0.) GO TO 120
          A = A+D
          N = N+1
120     CONTINUE
        IF (N.LE.0) GO TO 140
        AN = N
        A = A/AN
140     AVGE = A
C
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE  Q801 (F ,T ,H)
C
C...... ROUTINE TO COMPUTE TOTAL HEAT OF COMBUSTION PRODUCTS (N.R.C.)
C......        ENTER WITH F = 0.0 FOR AIR
C...... ENTER WITH F = 0.06823 FOR STOICHIOMETRIC MIXTURE
C...... COEFFICIENTS FOR T.LE.800 DEG.K ARE AX(1) - CORRECTION FOR F.A.R.
C...... COEFFICIENTS FOR T.GT.800 DEG.K ARE AX(2) - CORRECTION FOR F.A.R.
C...... COEFFICIENTS FOR T.LE.800 DEG.K ARE AX(3) - TOTAL HEAT OF AIR
C...... COEFFICIENTS FOR T.GT.800 DEG.K ARE AX(4) - TOTAL HEAT OF AIR
C
        INTEGER I
        REAL F, T, H, A0(4), A1(4), A2(4), A3(4), A4(4), A5(4)
C
        DATA A0/ .1496069E+02,-.4223715E+02,-.4056710    , .1131763E+02/
        DATA A1/-.1263568    , .1998916    , .2433633    , .1907555    /
        DATA A2/ .7697103E-03, .8712192E-04,-.1646082E-04, .6376249E-04/
        DATA A3/-.5175851E-06, .6008276E-07, .1579838E-07,-.1821733E-07/
        DATA A4/ .5959612E-10,-.2995447E-10, .2531721E-10, .2234454E-11/
        DATA A5/ .8333768E-13, .3918570E-14,-.1797673E-13,     0.0     /
C
        DELH = 0.0
        I = 1
        IF (T .GT. 800.0) I = 2
        IF (F.LE.0.0) GO TO 5
3       CONTINUE
        H = ((((A5(I)*T +A4(I))*T +A3(I))*T + A2(I))*T +A1(I))*T +A0(I)
        IF (I.GE.2) GO TO 6
        DELH = H * F / (1.0 + F)
5       I = I + 2
        GO TO 3
6       H = H + DELH
C
        RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE Q803 (XI,YI,Y,TOL,ISW,N)
C
C...... LINEAR ITERATION SUBROUTINE
C
        INTEGER N, ISW
        REAL    X1, XI, YI, Y, TOL, E1, E2
C
        E2 = Y-YI
        N = 1
        IF (ABS(E2) .LE. TOL) GO TO 5
        N = 3
        IF (ISW - 1) 5 , 6 , 2
2       IF (ISW.GE.15) GO TO 5
        E1 = E1 - E2
        IF (E1.EQ.0.0) GO TO 5
        E1 = E2 * (X1-XI) / E1
4       X1 = XI
        XI = XI - E1
        E1 = E2
        ISW = ISW + 1
        N = 2
5       RETURN
6       E1 = -.1 * XI
        GO TO 4
        END
C
C*****************************************************************************
C
        SUBROUTINE Q805 (T1,T2,CV,FAR,N)
C
CNOTE:  N=1 OK  N=2 T2 LESS THAN T1
C       FUEL TO AIR RATIO FROM T1,T2 AND CV
C       CURRENT VALUE OF STOICHIOMETRIC FUEL/AIR RATIO
C
        INTEGER N
        REAL    H1, H2, T1, T2, CV, FAR, FS, STH2
C
        DATA FS /0.06823/
C
        FAR = 0.0
        N = 2
        IF (T2 - T1) 3 , 2 , 1
1       CONTINUE
        CALL Q801 (0.0,T1,H1)
        CALL Q801 (0.0,T2,H2)
        CALL Q801 (FS ,T2,STH2)
        FAR = (H2 - H1) / (CV + 99.5 + (H2 - STH2) / FS - STH2)
2       N = 1
3       RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE Q806 (FAR,T1,CV,T2)
C
C...... TEMPERATURE AFTER BURNING S/R  (USING FUEL CAL.VAL.)
C
        INTEGER IC, N
        REAL    FAR, CV, T1, T2, DT, TOL, F
C
        T2 = T1
        IF (FAR) 3 , 4 , 1
1       DT = FAR * (44250.0 + (FAR * (125.0*T1 - 225000.0) - 11.25*T1))
        TOL = 0.00001 * FAR + 5.0E-10
        IC = 1
2       T2 = T1 + DT
        CALL Q805 (T1 ,T2 ,CV ,F ,N)
        CALL Q803 (DT ,F ,FAR ,TOL ,IC ,N)
        IF (N-2) 4 , 2 , 3
3       T2 = 0.0
4       RETURN
        END
C
C*****************************************************************************
C
        SUBROUTINE QQ966 (XTAB,YTAB,NX,X,Y)
C
        INTEGER NX, I, J
        REAL XTAB(1), YTAB(1), X, Y
C
        DO 10 I = 2,NX
          IF (XTAB(I).GE.X) GO TO 20
10      CONTINUE
        I = NX
20      J = I-1
        Y = YTAB(J) + (X-XTAB(J))*(YTAB(I)-YTAB(J)) / (XTAB(I)-XTAB(J))
C
        RETURN
        END
        SUBROUTINE SINGLE
C
C       SINGLE PRECISION FLOATING POINT BENCHMARK
C
        REAL WINT(91),X1,X2,UP,STEP
        COMMON /SINGL/ WINT
        COMMON /LF/ LOOPFACTOR
C
        F(X) = (((X-2.)*X+3.)*X-4.)*X+5.
C
        DO 40 LOOP=1,(10*LOOPFACTOR)
          DO 10 I=1,91
            WINT(I) = 0.
10        CONTINUE
          DO 30 J = 1,91,10
            X2 = 0.
            UP = 1. + (J-1) / 10.
            STEP = UP / 1000.
            DO 20 I=1,1000
              X1 = X2
              X2 = I * STEP
              WINT(J) = WINT(J) + STEP
     1                  * (F(X1)+F(X2)+4*F((X1+X2)/2.))/6.
20          CONTINUE
30        CONTINUE
40      CONTINUE
C
        RETURN
        END
        SUBROUTINE SWHET
C
C       WHETSTONE SYNTHETIC BENCHMARK PROGRAM, SINGLE-PRECISION VARIANT
C
C       THIS PROGRAM IS THE BASIS OF THE WHETSTONE-INSTRUCTIONS-PER-SECOND
C       METRIC FOR FORTRAN AND CPU PERFORMANCE.
C
C       FOR MORE INFORMATION, SEE
C           "A SYNTHETIC BENCHMARK," H. J. CURNOW AND B. A. WICHMANN,
C           THE COMPUTER JOURNAL, VOLUME 19, NO. 1, PAGES 43-49.
C
C       NOTE:  EARLIER VERSIONS OF THIS CODING OF THE WHETSTONE BENCHMARK
C               COMPUTED N4, N8, AND N9 BELOW AS INTEGERS > 32767.
C               AS A CONSEQUENCE, THEIR VALUES ON 16-BIT PROCESSORS WERE
C               INCORRECTLY COMPUTED.  THE PROBLEM HAS BEEN FIXED IN
C               THIS VERSION BY PUTTING AN ADDITIONAL LOOP OF 10 AROUND
C               THE APPROPRIATE TESTS AND DECREASING N4, N8, AND N9 BY
C               10.  ALSO ALL VARIABLES HAVE BEEN EXPLICITLY TYPED.
C
        REAL X1,X2,X3,X4,X,Y,Z,T,T1,T2,E1
        INTEGER I,II,J,K,L,N1,N2,N3,N4,N6,N7,N8,N9,N10,N11,LP
C
        COMMON /SWHETX/ T,T1,T2,E1(4),J,K,L
        COMMON /LF/ LOOPFACTOR
C
         I  = 100
        J  = I/10
        T  = 0.499975E00
        T1 = 0.50025E00
        T2 = 2.0E00
C
C
        N1  = 0
        N2  = I * 12
        N3  = I * 14
        N4  = J * 345
        N6  = I * 210
        N7  = I * 32
        N8  = J * 899
        N9  = J * 616
        N10 = 0
        N11 = I * 93
C
        DO 2000 LP=1,(1*LOOPFACTOR)
        X1 =  1.0E0
        X2 = -1.0E0
        X3 = -1.0E0
        X4 = -1.0E0
        IF ( N1 .LE. 0 ) GO TO 120
           DO 110,I = 1,N1,1
              X1 = ( X1+X2+X3-X4) * T
              X2 = ( X1+X2-X3+X4) * T
              X4 = (-X1+X2+X3+X4) * T
              X3 = ( X1-X2+X3+X4) * T
110           CONTINUE
120        CONTINUE
C        CALL POUT (N1,N1,N1,X1,X2,X3,X4)
C
        E1(1) =  1.0E0
        E1(2) = -1.0E0
        E1(3) = -1.0E0
        E1(4) = -1.0E0
        IF ( N2 .LE. 0 ) GO TO 220
           DO 210 I = 1,N2,1
              E1(1) = (  E1(1)+E1(2)+E1(3)-E1(4) ) * T
              E1(2) = (  E1(1)+E1(2)-E1(3)+E1(4) ) * T
              E1(3) = (  E1(1)-E1(2)+E1(3)+E1(4) ) * T
              E1(4) = ( -E1(1)+E1(2)+E1(3)+E1(4) ) * T
210           CONTINUE
220        CONTINUE
C        CALL POUT (N2,N3,N2,E1(1),E1(2),E1(3),E1(4))
C
        IF ( N3 .LE. 0 ) GO TO 320
           DO 310 I = 1,N3,1
              CALL PA (E1)
310           CONTINUE
320        CONTINUE
C        CALL POUT (N3,N2,N2,E1(1),E1(2),E1(3),E1(4))
C
        J = 1
        IF ( N4 .LE. 0 ) GO TO 490
           DO 480 II = 1,10
              DO 470 I = 1,N4,1
                 IF ( J .NE. 1 ) GO TO 410
                 J = 2
                 GO TO 420
410                 J = 3
420                 IF ( J .LE. 2 ) GO TO 430
                 J = 0
                 GO TO 440
430                 J = 1
440                 IF ( J .GE. 1 ) GO TO 450
                 J = 1
                 GO TO 460
450                 J = 0
460                 CONTINUE
470              CONTINUE
480           CONTINUE
490        CONTINUE
C        CALL POUT (N4,J,J,X1,X2,X3,X4)
C
        J = 1
        K = 2
        L = 3
        IF ( N6 .LE. 0 ) GO TO 620
           DO 610 I = 1,N6,1
              J = J * (K-J) * (L-K)
              K = L * K - (L-J) * K
              L = (L-K) * (K+J)
              E1(L-1) = J + K + L
              E1(K-1) = J * K * L
610           CONTINUE
620        CONTINUE
C        CALL POUT (N6,J,K,E1(1),E1(2),E1(3),E1(4))
C
        X = 0.5E0
        Y = 0.5E0
        IF ( N7 .LE. 0 ) GO TO 720
           DO 710 I = 1,N7,1
           X = T*ATAN(T2*SIN(X)*COS(X)/(COS(X+Y)+COS(X-Y)-1.0E0))
           Y = T*ATAN(T2*SIN(Y)*COS(Y)/(COS(X+Y)+COS(X-Y)-1.0E0))
710        CONTINUE
720        CONTINUE
C        CALL POUT (N7,J,K,X,X,Y,Y)
C
        X = 1.0E0
        Y = 1.0E0
        Z = 1.0E0
        IF ( N8 .LE. 0 ) GO TO 830
           DO 820 II = 1,10
              DO 810 I = 1,N8,1
                 CALL P3 (X,Y,Z)
810              CONTINUE
820           CONTINUE
830        CONTINUE
C        CALL POUT (N8,J,K,X,Y,Z,Z)
C
        J = 1
        K = 2
        L = 3
        E1(1) = 1.0E0
        E1(2) = 2.0E0
        E1(3) = 3.0E0
        IF ( N9 .LE. 0 ) GO TO 930
           DO 920 II = 1,10
              DO 910 I = 1,N9,1
                 CALL P0
910              CONTINUE
920           CONTINUE
930        CONTINUE
C        CALL POUT (N9,J,K,E1(1),E1(2),E1(3),E1(4))
C
        J = 2
        K = 3
        IF ( N10 .LE. 0 ) GO TO 1020
           DO 1010 I = 1,N10,1
              J = J + K
              K = J + K
              J = J - K
              K = K - J - J
1010           CONTINUE
1020        CONTINUE
C        CALL POUT (N10,J,K,X1,X2,X3,X4)
C
        X = 0.75E0
        IF ( N11 .LE. 0 ) GO TO 1120
           DO 1110 I = 1,N11,1
              X = SQRT ( EXP( ALOG(X)/T1 ) )
1110           CONTINUE
1120        CONTINUE
C        CALL POUT (N11,J,K,X,X,X,X)
C
C
C...... KILOWHETSTONE INSTRUCTIONS PER SECOND IS EQUAL TO 10000/TIME
2000    CONTINUE
        RETURN
        END
C
C***********************************************************************
C
        SUBROUTINE PA (E)
C
        COMMON /SWHETX/ T,T1,T2,E1(4),J,K,L
        REAL E(4),T,T1,T2
        INTEGER J
C
        J = 0
1        E(1) = (  E(1)+E(2)+E(3)-E(4) ) * T
        E(2) = (  E(1)+E(2)-E(3)+E(4) ) * T
        E(3) = (  E(1)-E(2)+E(3)+E(4) ) * T
        E(4) = ( -E(1)+E(2)+E(3)+E(4) ) / T2
        J = J + 1
        IF ( J .LT. 0 ) GO TO 1
        RETURN  
        END
C
C***********************************************************************
C
        SUBROUTINE P0
C
        COMMON /SWHETX/ T,T1,T2,E1(4),J,K,L
        REAL T,T1,T2,E1
        INTEGER J,K,L
C
        E1(J) = E1(K)
        E1(K) = E1(L)
        E1(L) = E1(J)
C
        RETURN
        END
C
C***********************************************************************
C
        SUBROUTINE P3 (X,Y,Z)
C
        COMMON /SWHETX/ T,T1,T2,E1(4),J,K,L
        REAL T,T1,T2,X,Y,Z
C
        X1 = X
        Y1 = Y
        X1 = T * (X1+Y1)
        Y1 = T * (X1+Y1)
        Z  = (X1+Y1) / T2
        RETURN
        END
C
C***********************************************************************
C
C        SUBROUTINE POUT (N,J,K,X1,X2,X3,X4)
C
C       WRITE STATEMENT COMMENTED OUT TO IMPROVE REPEATABILITY OF TIMINGS
C
C       WRITE(2,1)N,J,K,X1,X2,X3,X4
C1      FORMAT(1H,3I7,4E12.4)
C
C        RETURN
C        END
