program linpack_bench_main
!
!*******************************************************************************
!
!! LINPACK_BENCH_MAIN drives the LINPACK benchmark program.
!
  implicit none
!
  integer, parameter :: n = 1000
  integer, parameter :: lda = n + 1
!
  real a(1001,1000)
  real a_max
  real b(1000)
  real b_max
  real, parameter :: cray = 0.056E+00
  real eps
  integer i
  integer info
  integer ipvt(1000)
  real ops
  real resid(n)
  real resid_max
  real residn
  real t1
  real t2
  real time(6)
  real total
  real x(1000)
!
!     this program was updated on 10/12/92 to correct a
!     problem with the random number generator. The previous
!     random number generator had a short period and produced
!     singular matrices occasionally.
!
      ops = (2.0d0*float(n)**3)/3.0d0 + 2.0d0*float(n)**2

         call matgen(a,lda,n)

  a_max = maxval ( maxval ( a(1:n,1:n), dim = 2 ), dim = 1 )

  x(1:n) = 1.0E+00
  b(1:n) = matmul ( a(1:n,1:n), x(1:n) )

  call cpu_time ( t1 )

  call sgefa(a,lda,n,ipvt,info)

  call cpu_time ( t2 )
  time(1) = t2 - t1

  call cpu_time ( t1 )

  call sgesl(a,lda,n,ipvt,b,0)

  call cpu_time ( t2 )
  time(2) = t2 - t1

  total = time(1) + time(2)
!
!     compute a residual to verify results.
!
  call matgen ( a, lda, n )

  x(1:n) = 1.0E+00
  b(1:n) = matmul ( a(1:n,1:n), x(1:n) )

  resid(1:n) = b(1:n) - matmul ( a(1:n,1:n), x(1:n) )

  resid_max = maxval ( abs ( resid(1:n) ) )
  b_max = maxval ( abs ( b(1:n) ) ) 

  eps = epsilon ( eps )

         residn = resid_max /  ( n * a_max * b_max * eps )
         write(*,40)
   40    format('     norm. resid      resid           machep', &
               '         x(1)          x(n)')
         write(*,50) residn, resid_max, eps, x(1), x(n)
   50    format(1p5e16.8)

         write(*,60) n
   60    format(//'    times are reported for matrices of order ',i5)
         write(*,70)
   70    format(6x,'factor',5x,'solve',6x,'total',5x,'mflops',7x,'unit', &
              6x,'ratio')

         time(3) = total
         time(4) = ops/(1.0e6*total)
         time(5) = 2.0e0/time(4)
         time(6) = total/cray
         write(*,80) lda
   80    format(' times for array with leading dimension of',i4)
         write(*,110) (time(i),i=1,6)
  110    format(6(1pe11.3))
         write(*,*)' end of tests -- this version dated 10/12/92'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'LINPACK_BENCH'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
function isamax ( n, x, incx )
!
!*******************************************************************************
!
!! ISAMAX finds the index of the vector element of maximum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer ISAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer isamax
  integer ix
  integer n
  real samax
  real x(*)
!
  if ( n <= 0 ) then

    isamax = 0

  else if ( n == 1 ) then

    isamax = 1

  else if ( incx == 1 ) then

    isamax = 1
    samax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > samax ) then
        isamax = i
        samax = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    isamax = 1
    samax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        isamax = i
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
subroutine matgen ( a, lda, n )
!
!*******************************************************************************
!
!! MATGEN generates a random matrix and suitable right hand side.
!
  implicit none
!
      integer lda,n,i,j,init(4)
      real a(lda,1),ran
!
      init(1) = 1
      init(2) = 2
      init(3) = 3
      init(4) = 1325

  do j = 1,n
     do i = 1,n
       a(i,j) = ran(init) - .5
    end do
  end do

  return
end
subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real values.
!
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none
!
  real x
  real y
  real z
!
  z = x
  x = y
  y = z

  return
end
subroutine saxpy ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SAXPY adds a constant times one vector to another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), the vector to which a multiple of X is to
!    be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( sa == 0.0E+00 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = y(1:n) + sa * x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine sgefa ( a, lda, n, ipvt, info )
!
!*******************************************************************************
!
!! SGEFA factors a real matrix.
!
!
!  Modified:
!
!    07 March 2001
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Parameters:
!
!    Input/output, real A(LDA,N).
!    On intput, the matrix to be factored.
!    On output, an upper triangular matrix and the multipliers used to obtain
!    it.  The factorization can be written A=L*U, where L is a product of
!    permutation and unit lower triangular matrices, and U is upper triangular.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Output, integer IPVT(N), the pivot indices.
!
!    Output, integer INFO, singularity indicator.
!    0, normal value.
!    K, if U(K,K) == 0.  This is not an error condition for this subroutine,
!    but it does indicate that SGESL or SGEDI will divide by zero if called.
!    Use RCOND in SGECO for a reliable indication of singularity.
!
  implicit none
!
  integer lda
  integer n
!
  real a(lda,n)
  integer info
  integer ipvt(n)
  integer isamax
  integer j
  integer k
  integer l
  real t
!
!  Gaussian elimination with partial pivoting.
!
  info = 0

  do k = 1, n - 1
!
!  Find L = pivot index.
!
    l = isamax ( n-k+1, a(k,k), 1 ) + k - 1
    ipvt(k) = l
!
!  Zero pivot implies this column already triangularized.
!
    if ( a(l,k) == 0.0E+00 ) then
      info = k
      cycle
    end if
!
!  Interchange if necessary.
!
    if ( l /= k ) then
      call r_swap ( a(l,k), a(k,k) )
    end if
!
!  Compute multipliers.
!
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k)
!
!  Row elimination with column indexing.
!
    do j = k+1, n
      t = a(l,j)
      if ( l /= k ) then
        a(l,j) = a(k,j)
        a(k,j) = t
      end if
      call saxpy ( n-k, t, a(k+1,k), 1, a(k+1,j), 1 )
    end do

  end do

  ipvt(n) = n

  if ( a(n,n) == 0.0E+00 ) then
    info = n
  end if

  return
end
subroutine sgesl ( a, lda, n, ipvt, b, job )
!
!*******************************************************************************
!
!! SGESL solves a real general linear system A * X = B.
!
!
!  Discussion:
!
!    SGESL can solve either of the systems A * X = B or A' * X = B.
!
!    The system matrix must have been factored by SGECO or SGEFA.
!
!    A division by zero will occur if the input factor contains a
!    zero on the diagonal.  Technically this indicates singularity
!    but it is often caused by improper arguments or improper
!    setting of LDA.  It will not occur if the subroutines are
!    called correctly and if SGECO has set RCOND > 0.0E+00
!    or SGEFA has set INFO == 0.
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Modified:
!
!    07 March 2001
!
!  Parameters:
!
!    Input, real A(LDA,N), the output from SGECO or SGEFA.
!
!    Input, integer LDA, the leading dimension of A.
!
!    Input, integer N, the order of the matrix A.
!
!    Input, integer IPVT(N), the pivot vector from SGECO or SGEFA.
!
!    Input/output, real B(N).
!    On input, the right hand side vector.
!    On output, the solution vector.
!
!    Input, integer JOB.
!    0, solve A * X = B;
!    nonzero, solve A' * X = B.
!
  implicit none
!
  integer lda
  integer n
!
  real a(lda,n)
  real b(n)
  integer ipvt(n)
  integer job
  integer k
  integer l
  real t
!
!  Solve A * X = B.
!
  if ( job == 0 ) then

    do k = 1, n-1

      l = ipvt(k)
      t = b(l)

      if ( l /= k ) then
        b(l) = b(k)
        b(k) = t
      end if

      call saxpy ( n-k, t, a(k+1,k), 1, b(k+1), 1 )

    end do

    do k = n, 1, -1
      b(k) = b(k) / a(k,k)
      t = -b(k)
      call saxpy ( k-1, t, a(1,k), 1, b(1), 1 )
    end do

  else
!
!  Solve A' * X = B.
!
    do k = 1, n
      t = dot_product ( a(1:k-1,k), b(1:k-1) )
      b(k) = ( b(k) - t ) / a(k,k)
    end do

    do k = n-1, 1, -1

      b(k) = b(k) + dot_product ( a(k+1:n,k), b(k+1:n) )
      l = ipvt(k)

      if ( l /= k ) then
        call r_swap ( b(l), b(k) )
      end if

    end do

  end if

  return
end
FUNCTION RAN( ISEED )
!
!*******************************************************************************
!
!! RAN returns a uniformly distributed random number between 0 and 1.
!
!
!  Arguments
!  =========
!
!  ISEED   (input/output) INTEGER array, dimension (4)
!          On entry, the seed of the random number generator; the array
!          elements must be between 0 and 4095, and ISEED(4) must be
!          odd.
!          On exit, the seed is updated.
!
!  Further Details
!  ===============
!
!  This routine uses a multiplicative congruential method with modulus
!  2**48 and multiplier 33952834046453 (see G.S.Fishman,
!  'Multiplicative congruential random number generators with modulus
!  2**b: an exhaustive analysis for b = 32 and a partial analysis for
!  b = 48', Math. Comp. 189, pp 331-344, 1990).
!
!  48-bit integers are stored in 4 integer array elements with 12 bits
!  per element. Hence the routine is portable across machines with
!  integers of 32 bits or more.
!
  implicit none
!
      INTEGER            ISEED( 4 )
!     .. Parameters ..
      INTEGER            M1, M2, M3, M4
      PARAMETER          ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 )
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
      INTEGER            IPW2
      REAL               R
      PARAMETER          ( IPW2 = 4096, R = ONE / IPW2 )

      INTEGER            IT1, IT2, IT3, IT4
  real ran
!
!  Multiply the seed by the multiplier modulo 2**48
!
      IT4 = ISEED( 4 )*M4
      IT3 = IT4 / IPW2
      IT4 = IT4 - IPW2*IT3
      IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3
      IT2 = IT3 / IPW2
      IT3 = IT3 - IPW2*IT2
      IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2
      IT1 = IT2 / IPW2
      IT2 = IT2 - IPW2*IT1
      IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 +ISEED( 4 )*M1
      IT1 = MOD( IT1, IPW2 )
!
!  Return updated seed
!
      ISEED( 1 ) = IT1
      ISEED( 2 ) = IT2
      ISEED( 3 ) = IT3
      ISEED( 4 ) = IT4
!
!  Convert 48-bit integer to a real number in the interval (0,1)
!
      RAN = R*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R* &
              ( REAL( IT4 ) ) ) ) )

  RETURN
end