function linpack_bench ( n, lda )

%*****************************************************************************80
%
%% linpack_bench() drives the LINPACK benchmark program.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 February 2025
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt.
%
%  Input:
%
%    integer N, the order of the matrix.  If no value is supplied,
%    the default value of 1000 is used.
%
%    integer LDA, the leading dimension of the matrix.  
%    If no value is supplied, the default value of N+1 is used.
%
  timestamp ( );

  if ( nargin < 1 ) 
    n = 1000;
  end

  if ( nargin < 2 )
    lda = n + 1;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'linpack_bench():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  The LINPACK benchmark.\n' );
  fprintf ( 1, '  Datatype: Real Double Precision\n' );
  fprintf ( 1, '  Matrix order N               = %d\n', n );
  fprintf ( 1, '  Leading matrix dimension LDA = %d\n', lda );

  A = r8_matgen ( lda, n );
  x_exact(1:n,1) = 1.0;
  b = A * x_exact;

  ALU = A;

  t = cputime ( );
  [ ALU, ipvt, info ] = dgefa ( ALU, lda, n );
  t = cputime ( ) - t;

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'linpack_bench(): Fatal error!\n' );
    fprintf ( 1, '  The matrix A is apparently singular.\n' );
    fprintf ( 1, '  Abnormal end of execution.\n' );
    error ( 'linpack_bench(): Fatal error!' );
  end

  x = b';
  job = 0;

  t2 = cputime ( );
  x = dgesl ( ALU, lda, n, ipvt, x, job );
  t2 = cputime ( ) - t2;
  t = t + t2;
%
%  Compute the residual.
%
  r = b - A * x';
%
%  Compute infinity norms.
%
  A_norm = norm ( A, Inf );
  x_norm = norm ( x, Inf );
  r_norm = norm ( r, Inf );
%
%  Report.
%
  ratio = r_norm / A_norm / x_norm / eps;

  ops = ( 2 * n * n * n ) / 3.0 + 2.0 * ( n * n );
  mflops = ops / ( 1000000 * t );
  unit = 2.0 / mflops;

  cray = 0.056;
  cray_ratio = t / cray;

  fprintf ( 1, "\n" );
  fprintf ( 1, "  Normalized residual = %14f \n", ratio );
  fprintf ( 1, "  Residual norm       = %14f \n", r_norm );
  fprintf ( 1, "  Machine epsilon     = %14f \n", eps );
  fprintf ( 1, "  First X[]           = %14f \n", x(1) );
  fprintf ( 1, "  Last X[]            = %14f \n", x(n) );
  fprintf ( 1, "  Time in seconds     = %14f \n", t );
  fprintf ( 1, "  MegaFLOPS           = %14f \n", mflops );
  fprintf ( 1, "  Unit                = %14f \n", unit );
  fprintf ( 1, "  Cray ratio          = %14f \n", cray_ratio );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'linpack_bench():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function y = daxpy ( n, sa, x, incx, y, incy )

%*****************************************************************************80
%
%% daxpy() adds a constant times one vector to another.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt
%
%  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
%
%  Input:
%
%    integer N, the number of entries in the vector.
%
%    real SA, the multiplier.
%
%    Input, real X(*), the vector to be scaled and added to Y.
%
%    integer INCX, the increment between successive entries of X.
%
%    real Y(*), the vector to which a multiple of X is to be added.
%
%    integer INCY, the increment between successive entries of Y.
%
%  Output:
%
%    real Y(*), the updated vector.
%
  if ( n <= 0 )

  elseif ( sa == 0.0 )

  elseif ( incx == 1 && incy == 1 )

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

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      y(iy) = y(iy) + sa * x(ix);
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end
function [ a, ipvt, info ] = dgefa ( a, lda, n )

%*****************************************************************************80
%
%% dgefa() factors a real matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt
%
%  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
%
%  Input:
%
%    real A(LDA,N), the matrix to be factored.
%
%    integer LDA, the leading dimension of A.
%
%    integer N, the order of the matrix A.
%
%  Output:
%
%    real A(LDA,N), 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.
%
%    integer IPVT(N), the pivot indices.
%
%    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 DGESL or DGEDI will divide by zero if called.
%    Use RCOND in DGECO for a reliable indication of singularity.
%

%
%  Gaussian elimination with partial pivoting.
%
  info = 0;

  for k = 1 : n - 1
%
%  Find L = pivot index.
%
    l = idamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;
%
%  Zero pivot implies this column already triangularized.
%
    if ( a(l,k) == 0.0 )
      info = k;
      continue
    end
%
%  Interchange if necessary.
%
    if ( l ~= k )
      t      = a(l,k);
      a(l,k) = a(k,k);
      a(k,k) = t;
    end
%
%  Compute multipliers.
%
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k);
%
%  Row elimination with column indexing.
%
    for j = k+1 : n
      t = a(l,j);
      if ( l ~= k )
        a(l,j) = a(k,j);
        a(k,j) = t;
      end
      a(k+1:n,j) = daxpy ( n-k, t, a(k+1:n,k), 1, a(k+1:n,j), 1 );
    end

  end

  ipvt(n) = n;

  if ( a(n,n) == 0.0 )
    info = n;
  end

  return
end
function b = dgesl ( a, lda, n, ipvt, b, job )

%*****************************************************************************80
%
%% dgesl() solves a real general linear system A * X = B.
%
%  Discussion:
%
%    DGESL can solve either of the systems A * X = B or A' * X = B.
%
%    The system matrix must have been factored by DGECO or DGEFA.
%
%    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 DGECO has set 0.0 < RCOND 
%    or DGEFA has set INFO == 0.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt.
%
%  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
%
%  Input:
%
%    real A(LDA,N), the output from DGECO or DGEFA.
%
%    integer LDA, the leading dimension of A.
%
%    integer N, the order of the matrix A.
%
%    integer IPVT(N), the pivot vector from DGECO or DGEFA.
%
%    real B(N), the right hand side vector.
%
%    integer JOB.
%    0, solve A * X = B;
%    nonzero, solve A' * X = B.
%
%  Output:
%
%    real B(N), the solution vector.
%

%
%  Solve A * X = B.
%
  if ( job == 0 )

    for k = 1 : n-1

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

      if ( l ~= k )
        b(l) = b(k);
        b(k) = t;
      end

      b(k+1:n) = daxpy ( n-k, t, a(k+1:n,k)', 1, b(k+1:n), 1 );

    end

    for k = n : -1 : 1
      b(k) = b(k) / a(k,k);
      t = -b(k);
      b(1:k-1) = daxpy ( k-1, t, a(1:k-1,k)', 1, b(1:k-1), 1 );
    end

  else
%
%  Solve A' * X = B.
%
    for k = 1 : n
      t = a(1:k-1,k)' * b(1:k-1);
      b(k) = ( b(k) - t ) / a(k,k);
    end

    for k = n-1 : -1 : 1

      b(k) = b(k) + a(k+1:n,k)' * b(k+1:n);
      l = ipvt(k);

      if ( l ~= k )
        t    = b(l);
        b(l) = b(k);
        b(k) = t;
      end

    end

  end

  return
end
function value = idamax ( n, x, incx )

%*****************************************************************************80
%
%% idamax() finds the index of the vector element of maximum absolute value.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt
%
%  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
%
%  Input:
%
%    integer N, the number of entries in the vector.
%
%    real X(*), the vector to be examined.
%
%    integer INCX, the increment between successive entries of SX.
%
%  Output:
%
%    integer VALUE, the index of the element of SX of maximum
%    absolute value.
%
  if ( n <= 0 )

    value = 0;

  elseif ( n == 1 )

    value = 1;

  elseif ( incx == 1 )

    value = 1;
    samax = abs ( x(1) );

    for i = 2 : n

      if ( samax < abs ( x(i) ) )
        value = i;
        samax = abs ( x(i) );
      end

    end

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    value = 1;
    samax = abs ( x(ix) );

    ix = ix + incx;

    for i = 2 : n
      if ( samax < abs ( x(ix) ) )
        value = i;
        samax = abs ( x(ix) );
      end
      ix = ix + incx;
    end

  end

  return
end
function a = r8_matgen ( lda, n )

%*****************************************************************************80
%
%% r8_matgen() generates a random matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.
%    This version by John Burkardt
%
%  Input:
%
%    integer LDA, the leading dimension of the matrix.
%
%    integer N, the order of the matrix.
%
%  Output:
%
%    real A(LDA,N), the N by N matrix.
%
  a = rand ( n, n ) - 0.5;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
 
