lLinpackBenchmark2017.ino

//https://gist.githubusercontent.com/projectgus/8279947/raw/5746ada75a0442d14affa789fbff0f1a5ec5d52f/linpack_bench.ino

# include <stdlib.h>

# include <stdio.h>

# include <math.h>

int do_benchmark ( void );

double cpu_time ( void );

void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy );

double ddot ( int n, double dx[], int incx, double dy[], int incy );

int dgefa ( double a[], int lda, int n, int ipvt[] );

void dgesl ( double a[], int lda, int n, int ipvt[], double b[], int job );

void dscal ( int n, double sa, double x[], int incx );

int idamax ( int n, double dx[], int incx );

double r8_abs ( double x );

double r8_epsilon ( void );

double r8_max ( double x, double y );

double r8_random ( int iseed[4] );

double *r8mat_gen ( int lda, int n );

static FILE uartout = {0} ;

static int uart_putchar (char c, FILE *stream)

{

  Serial.write(c) ;

  return 0 ;

}

void setup() {

  Serial.begin(9600);

//  fdev_setup_stream (&uartout, uart_putchar, NULL, _FDEV_SETUP_WRITE);

//  stdout = &uartout ;

}

void loop() {

  printf("Starting benchmark...\n");

  do_benchmark();

}

/******************************************************************************/

int do_benchmark ( void )

/******************************************************************************/

/*

  Purpose:

    MAIN is the main program for LINPACK_BENCH.

  Discussion:

    LINPACK_BENCH drives the double precision LINPACK benchmark program.

  Modified:

    25 July 2008

  Parameters:

    N is the problem size.

*/

{

# define N 8

# define LDA ( N + 1 )

  static double a[N * LDA];

  static double a_max;

  static double b[N];

  static double b_max;

  const double cray = 0.056;

  static double eps;

  int i;

  int info;

  static int ipvt[N];

  int j;

  int job;

  double ops;

  static double resid[N];

  double resid_max;

  double residn;

  static double rhs[N];

  double t1;

  double t2;

  static double time[6];

  double total;

  double x[N];

  printf ( "\n" );

  printf ( "LINPACK_BENCH\n" );

  printf ( "  C version\n" );

  printf ( "\n" );

  printf ( "  The LINPACK benchmark.\n" );

  printf ( "  Language: C\n" );

  printf ( "  Datatype: Double precision real\n" );

  printf ( "  Matrix order N               = %d\n", N );

  printf ( "  Leading matrix dimension LDA = %d\n", LDA );

  ops = ( double ) ( 2L * N * N * N ) / 3.0 + 2.0 * ( double ) ( (long)N * N );

/*

  Allocate space for arrays.

*/

  r8mat_gen ( LDA, N, a);

  a_max = 0.0;

  for ( j = 0; j < N; j++ )

  {

    for ( i = 0; i < N; i++ )

    {

      a_max = r8_max ( a_max, a[i+j*LDA] );

    }

  }

  for ( i = 0; i < N; i++ )

  {

    x[i] = 1.0;

  }

  for ( i = 0; i < N; i++ )

  {

    b[i] = 0.0;

    for ( j = 0; j < N; j++ )

    {

      b[i] = b[i] + a[i+j*LDA] * x[j];

    }

  }

  t1 = cpu_time ( );

  info = dgefa ( a, LDA, N, ipvt );

  if ( info != 0 )

  {

    printf ( "\n" );

    printf ( "LINPACK_BENCH - Fatal error!\n" );

    printf ( "  The matrix A is apparently singular.\n" );

    printf ( "  Abnormal end of execution.\n" );

    return 1;

  }

  t2 = cpu_time ( );

  time[0] = t2 - t1;

  t1 = cpu_time ( );

  job = 0;

  dgesl ( a, LDA, N, ipvt, b, job );

  t2 = cpu_time ( );

  time[1] = t2 - t1;

  total = time[0] + time[1];

/*

  Compute a residual to verify results.

*/

  r8mat_gen ( LDA, N, a );

  for ( i = 0; i < N; i++ )

  {

    x[i] = 1.0;

  }

  for ( i = 0; i < N; i++ )

  {

    rhs[i] = 0.0;

    for ( j = 0; j < N; j++ )

    {

      rhs[i] = rhs[i] + a[i+j*LDA] * x[j];

    }

  }

  for ( i = 0; i < N; i++ )

  {

    resid[i] = -rhs[i];

    for ( j = 0; j < N; j++ )

    {

      resid[i] = resid[i] + a[i+j*LDA] * b[j];

    }

  }

  resid_max = 0.0;

  for ( i = 0; i < N; i++ )

  {

    resid_max = r8_max ( resid_max, r8_abs ( resid[i] ) );

  }

  b_max = 0.0;

  for ( i = 0; i < N; i++ )

  {

    b_max = r8_max ( b_max, r8_abs ( b[i] ) );

  }

  eps = r8_epsilon ( );

  residn = resid_max / ( double ) N / a_max / b_max / eps;

  time[2] = total;

  if ( 0.0 < total )

  {

    time[3] = ops / ( 1.0E+06 * total );

  }

  else

  {

    time[3] = -1.0;

  }

  time[4] = 2.0 / time[3];

  time[5] = total / cray;

  printf ( "\n" );

  printf ( "     Norm. Resid      Resid           MACHEP         X[1]          X[N]\n" );

  printf ( "\n" );

  Serial.print("  ");

  Serial.print(residn, 14);

  Serial.print("  ");

  Serial.print(resid_max, 14);

  Serial.print("  ");

  Serial.print(eps, 14);

  Serial.print("  ");

  Serial.print(b[0], 14);

  Serial.print("  ");

  Serial.print(b[N-1], 14);

  Serial.print("  ");

  //printf ( "  %14f  %14f  %14e  %14f  %14f\n", residn, resid_max, eps, b[0], b[N-1] );

  printf ( "\n" );

  printf ( "      Factor     Solve      Total     MFLOPS       Unit      Cray-Ratio\n" );

  printf ( "\n" );

  for(int i =0; i<6;i++) {

      Serial.print("  ");

      Serial.print(time[i], 9);

  }

  //printf ( "  %9f  %9f  %9f  %9f  %9f  %9f\n", 

  //  time[0], time[1], time[2], time[3], time[4], time[5] );

/*

  Terminate.

*/

  printf ( "\n" );

  printf ( "LINPACK_BENCH\n" );

  printf ( "  Normal end of execution.\n" );

  printf ( "\n" );

  return 0;

# undef LDA

# undef N

}

/******************************************************************************/

double cpu_time ( void )

/******************************************************************************/

/*

  Purpose:

    CPU_TIME returns the current reading on the CPU clock.

  Discussion:

    The CPU time measurements available through this routine are often

    not very accurate.  In some cases, the accuracy is no better than

    a hundredth of a second.  

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 June 2005

  Author:

    John Burkardt

  Parameters:

    Output, double CPU_TIME, the current reading of the CPU clock, in seconds.

*/

{

  double value;

  value = ( double ) micros ( ) 

        / ( double ) 1000000;

  return value;

}

/******************************************************************************/

void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )

/******************************************************************************/

/*

  Purpose:

    DAXPY computes constant times a vector plus a vector.

  Discussion:

    This routine uses unrolled loops for increments equal to one.

  Modified:

    30 March 2007

  Author:

    FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.

    C version by John Burkardt

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,

    LINPACK User's Guide,

    SIAM, 1979.

    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,

    Basic Linear Algebra Subprograms for Fortran Usage,

    Algorithm 539, 

    ACM Transactions on Mathematical Software, 

    Volume 5, Number 3, September 1979, pages 308-323.

  Parameters:

    Input, int N, the number of elements in DX and DY.

    Input, double DA, the multiplier of DX.

    Input, double DX[*], the first vector.

    Input, int INCX, the increment between successive entries of DX.

    Input/output, double DY[*], the second vector.

    On output, DY[*] has been replaced by DY[*] + DA * DX[*].

    Input, int INCY, the increment between successive entries of DY.

*/

{

  int i;

  int ix;

  int iy;

  int m;

  if ( n <= 0 )

  {

    return;

  }

  if ( da == 0.0 )

  {

    return;

  }

/*

  Code for unequal increments or equal increments

  not equal to 1.

*/

  if ( incx != 1 || incy != 1 )

  {

    if ( 0 <= incx )

    {

      ix = 0;

    }

    else

    {

      ix = ( - n + 1 ) * incx;

    }

    if ( 0 <= incy )

    {

      iy = 0;

    }

    else

    {

      iy = ( - n + 1 ) * incy;

    }

    for ( i = 0; i < n; i++ )

    {

      dy[iy] = dy[iy] + da * dx[ix];

      ix = ix + incx;

      iy = iy + incy;

    }

  }

/*

  Code for both increments equal to 1.

*/

  else

  {

    m = n % 4;

    for ( i = 0; i < m; i++ )

    {

      dy[i] = dy[i] + da * dx[i];

    }

    for ( i = m; i < n; i = i + 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];

    }

  }

  return;

}

/******************************************************************************/

double ddot ( int n, double dx[], int incx, double dy[], int incy )

/******************************************************************************/

/*

  Purpose:

    DDOT forms the dot product of two vectors.

  Discussion:

    This routine uses unrolled loops for increments equal to one.

  Modified:

    30 March 2007

  Author:

    FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.

    C version by John Burkardt

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,

    LINPACK User's Guide,

    SIAM, 1979.

    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,

    Basic Linear Algebra Subprograms for Fortran Usage,

    Algorithm 539, 

    ACM Transactions on Mathematical Software, 

    Volume 5, Number 3, September 1979, pages 308-323.

  Parameters:

    Input, int N, the number of entries in the vectors.

    Input, double DX[*], the first vector.

    Input, int INCX, the increment between successive entries in DX.

    Input, double DY[*], the second vector.

    Input, int INCY, the increment between successive entries in DY.

    Output, double DDOT, the sum of the product of the corresponding

    entries of DX and DY.

*/

{

  double dtemp;

  int i;

  int ix;

  int iy;

  int m;

  dtemp = 0.0;

  if ( n <= 0 )

  {

    return dtemp;

  }

/*

  Code for unequal increments or equal increments

  not equal to 1.

*/

  if ( incx != 1 || incy != 1 )

  {

    if ( 0 <= incx )

    {

      ix = 0;

    }

    else

    {

      ix = ( - n + 1 ) * incx;

    }

    if ( 0 <= incy )

    {

      iy = 0;

    }

    else

    {

      iy = ( - n + 1 ) * incy;

    }

    for ( i = 0; i < n; i++ )

    {

      dtemp = dtemp + dx[ix] * dy[iy];

      ix = ix + incx;

      iy = iy + incy;

    }

  }

/*

  Code for both increments equal to 1.

*/

  else

  {

    m = n % 5;

    for ( i = 0; i < m; i++ )

    {

      dtemp = dtemp + dx[i] * dy[i];

    }

    for ( i = m; i < n; i = i + 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];

    }

  }

  return dtemp;

}

/******************************************************************************/

int dgefa ( double a[], int lda, int n, int ipvt[] )

/******************************************************************************/

/*

  Purpose:

    DGEFA factors a real general matrix.

  Modified:

    16 May 2005

  Author:

    C version by John Burkardt.

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch and Pete 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, double 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, int LDA, the leading dimension of A.

    Input, int N, the order of the matrix A.

    Output, int IPVT[N], the pivot indices.

    Output, int DGEFA, 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.

*/

{

  int info;

  int j;

  int k;

  int l;

  double t;

/*

  Gaussian elimination with partial pivoting.

*/

  info = 0;

  for ( k = 1; k <= n-1; k++ )

  {

/*

  Find L = pivot index.

*/

    l = idamax ( n-k+1, a+(k-1)+(k-1)*lda, 1 ) + k - 1;

    ipvt[k-1] = l;

/*

  Zero pivot implies this column already triangularized.

*/

    if ( a[l-1+(k-1)*lda] == 0.0 )

    {

      info = k;

      continue;

    }

/*

  Interchange if necessary.

*/

    if ( l != k )

    {

      t = a[l-1+(k-1)*lda];

      a[l-1+(k-1)*lda] = a[k-1+(k-1)*lda];

      a[k-1+(k-1)*lda] = t;

    }

/*

  Compute multipliers.

*/

    t = -1.0 / a[k-1+(k-1)*lda];

    dscal ( n-k, t, a+k+(k-1)*lda, 1 );

/*

  Row elimination with column indexing.

*/

    for ( j = k+1; j <= n; j++ )

    {

      t = a[l-1+(j-1)*lda];

      if ( l != k )

      {

        a[l-1+(j-1)*lda] = a[k-1+(j-1)*lda];

        a[k-1+(j-1)*lda] = t;

      }

      daxpy ( n-k, t, a+k+(k-1)*lda, 1, a+k+(j-1)*lda, 1 );

    }

  }

  ipvt[n-1] = n;

  if ( a[n-1+(n-1)*lda] == 0.0 )

  {

    info = n;

  }

  return info;

}

/******************************************************************************/

void dgesl ( double a[], int lda, int n, int ipvt[], double b[], int job )

/******************************************************************************/

/*

  Purpose:

    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.

  Modified:

    16 May 2005

  Author:

    C version by John Burkardt.

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch and Pete 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, double A[LDA*N], the output from DGECO or DGEFA.

    Input, int LDA, the leading dimension of A.

    Input, int N, the order of the matrix A.

    Input, int IPVT[N], the pivot vector from DGECO or DGEFA.

    Input/output, double B[N].

    On input, the right hand side vector.

    On output, the solution vector.

    Input, int JOB.

    0, solve A * X = B;

    nonzero, solve A' * X = B.

*/

{

  int k;

  int l;

  double t;

/*

  Solve A * X = B.

*/

  if ( job == 0 )

  {

    for ( k = 1; k <= n-1; k++ )

    {

      l = ipvt[k-1];

      t = b[l-1];

      if ( l != k )

      {

        b[l-1] = b[k-1];

        b[k-1] = t;

      }

      daxpy ( n-k, t, a+k+(k-1)*lda, 1, b+k, 1 );

    }

    for ( k = n; 1 <= k; k-- )

    {

      b[k-1] = b[k-1] / a[k-1+(k-1)*lda];

      t = -b[k-1];

      daxpy ( k-1, t, a+0+(k-1)*lda, 1, b, 1 );

    }

  }

/*

  Solve A' * X = B.

*/

  else

  {

    for ( k = 1; k <= n; k++ )

    {

      t = ddot ( k-1, a+0+(k-1)*lda, 1, b, 1 );

      b[k-1] = ( b[k-1] - t ) / a[k-1+(k-1)*lda];

    }

    for ( k = n-1; 1 <= k; k-- )

    {

      b[k-1] = b[k-1] + ddot ( n-k, a+k+(k-1)*lda, 1, b+k, 1 );

      l = ipvt[k-1];

      if ( l != k )

      {

        t = b[l-1];

        b[l-1] = b[k-1];

        b[k-1] = t;

      }

    }

  }

  return;

}

/******************************************************************************/

void dscal ( int n, double sa, double x[], int incx )

/******************************************************************************/

/*

  Purpose:

    DSCAL scales a vector by a constant.

  Modified:

    30 March 2007

  Author:

    FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.

    C version by John Burkardt

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,

    LINPACK User's Guide,

    SIAM, 1979.

    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,

    Basic Linear Algebra Subprograms for Fortran Usage,

    Algorithm 539,

    ACM Transactions on Mathematical Software,

    Volume 5, Number 3, September 1979, pages 308-323.

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double SA, the multiplier.

    Input/output, double X[*], the vector to be scaled.

    Input, int INCX, the increment between successive entries of X.

*/

{

  int i;

  int ix;

  int m;

  if ( n <= 0 )

  {

  }

  else if ( incx == 1 )

  {

    m = n % 5;

    for ( i = 0; i < m; i++ )

    {

      x[i] = sa * x[i];

    }

    for ( i = m; i < n; i = i + 5 )

    {

      x[i]   = sa * x[i];

      x[i+1] = sa * x[i+1];

      x[i+2] = sa * x[i+2];

      x[i+3] = sa * x[i+3];

      x[i+4] = sa * x[i+4];

    }

  }

  else

  {

    if ( 0 <= incx )

    {

      ix = 0;

    }

    else

    {

      ix = ( - n + 1 ) * incx;

    }

    for ( i = 0; i < n; i++ )

    {

      x[ix] = sa * x[ix];

      ix = ix + incx;

    }

  }

  return;

}

/******************************************************************************/

int idamax ( int n, double dx[], int incx )

/******************************************************************************/

/*

  Purpose:

    IDAMAX finds the index of the vector element of maximum absolute value.

  Discussion:

    WARNING: This index is a 1-based index, not a 0-based index!

  Modified:

    30 March 2007

  Author:

    FORTRAN77 original by Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart.

    C version by John Burkardt

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,

    LINPACK User's Guide,

    SIAM, 1979.

    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,

    Basic Linear Algebra Subprograms for Fortran Usage,

    Algorithm 539,

    ACM Transactions on Mathematical Software,

    Volume 5, Number 3, September 1979, pages 308-323.

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double X[*], the vector to be examined.

    Input, int INCX, the increment between successive entries of SX.

    Output, int IDAMAX, the index of the element of maximum

    absolute value.

*/

{

  double dmax;

  int i;

  int ix;

  int value;

  value = 0;

  if ( n < 1 || incx <= 0 )

  {

    return value;

  }

  value = 1;

  if ( n == 1 )

  {

    return value;

  }

  if ( incx == 1 )

  {

    dmax = r8_abs ( dx[0] );

    for ( i = 1; i < n; i++ )

    {

      if ( dmax < r8_abs ( dx[i] ) )

      {

        value = i + 1;

        dmax = r8_abs ( dx[i] );

      }

    }

  }

  else

  {

    ix = 0;

    dmax = r8_abs ( dx[0] );

    ix = ix + incx;

    for ( i = 1; i < n; i++ )

    {

      if ( dmax < r8_abs ( dx[ix] ) )

      {

        value = i + 1;

        dmax = r8_abs ( dx[ix] );

      }

      ix = ix + incx;

    }

  }

  return value;

}