# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "r8row.h"

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

int i4_log_10 ( int i )

/******************************************************************************/
/*
  Purpose:

    I4_LOG_10 returns the integer part of the logarithm base 10 of an I4.

  Example:

        I  I4_LOG_10
    -----  --------
        0    0
        1    0
        2    0
        9    0
       10    1
       11    1
       99    1
      100    2
      101    2
      999    2
     1000    3
     1001    3
     9999    3
    10000    4

  Discussion:

    I4_LOG_10 ( I ) + 1 is the number of decimal digits in I.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 October 2007

  Author:

    John Burkardt

  Parameters:

    Input, int I, the number whose logarithm base 10 is desired.

    Output, int I4_LOG_10, the integer part of the logarithm base 10 of
    the absolute value of X.
*/
{
  int i_abs;
  int ten_pow;
  int value;

  if ( i == 0 )
  {
    value = 0;
  }
  else
  {
    value = 0;
    ten_pow = 10;

    i_abs = abs ( i );

    while ( ten_pow <= i_abs )
    {
      value = value + 1;
      ten_pow = ten_pow * 10;
    }
  }
  return value;
}
/******************************************************************************/

int i4_max ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MAX returns the maximum of two I4's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, are two integers to be compared.

    Output, int I4_MAX, the larger of I1 and I2.
*/
{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

int i4_min ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MIN returns the smaller of two I4's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, two integers to be compared.

    Output, int I4_MIN, the smaller of I1 and I2.
*/
{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

int i4_power ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    I4_POWER returns the value of I^J.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 October 2007

  Author:

    John Burkardt

  Parameters:

    Input, int I, J, the base and the power.  J should be nonnegative.

    Output, int I4_POWER, the value of I^J.
*/
{
  int k;
  int value;

  if ( j < 0 )
  {
    if ( i == 1 )
    {
      value = 1;
    }
    else if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J negative.\n" );
      exit ( 1 );
    }
    else
    {
      value = 0;
    }
  }
  else if ( j == 0 )
  {
    if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J = 0.\n" );
      exit ( 1 );
    }
    else
    {
      value = 1;
    }
  }
  else if ( j == 1 )
  {
    value = i;
  }
  else
  {
    value = 1;
    for ( k = 1; k <= j; k++ )
    {
      value = value * i;
    }
  }
  return value;
}
/******************************************************************************/

void i4mat_print ( int m, int n, int a[], char *title )

/******************************************************************************/
/*
  Purpose:

    I4MAT_PRINT prints an I4MAT.

  Discussion:

    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, int A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void i4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    I4MAT_PRINT_SOME prints some of an I4MAT.

  Discussion:

    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    20 August 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, int A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 10

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of INCX.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:" );
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %6d", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to INCX) entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %6d", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

void i4vec_print ( int n, int a[], char *title )

/******************************************************************************/
/*
  Purpose:

    I4VEC_PRINT prints an I4VEC.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 November 2003

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, int A[N], the vector to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );

  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %6d: %8d\n", i, a[i] );
  }
  return;
}
/******************************************************************************/

int r8row_compare ( int m, int n, double a[], int i, int j )

/******************************************************************************/
/*
  Purpose:

    R8ROW_COMPARE compares rows in an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Example:

    Input:

      M = 4, N = 3, I = 2, J = 4

      A = (
        1  5  9
        2  6 10
        3  7 11
        4  8 12 )

    Output:

      R8ROW_COMPARE = -1

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A[M*N], the M by N array.

    Input, int I, J, the rows to be compared.
    I and J must be between 0 and M - 1.

    Output, int R8ROW_COMPARE, the results of the comparison:
    -1, row I < row J,
     0, row I = row J,
    +1, row J < row I.
*/
{
  int k;
  int value;
/*
  Check.
*/
  if ( i < 0 || m <= i )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8ROW_COMPARE - Fatal error!\n" );
    fprintf ( stderr, "  Row index I is out of bounds.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    exit ( 1 );
  }

  if ( j < 0 || m <= j )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8ROW_COMPARE - Fatal error!\n" );
    fprintf ( stderr, "  Row index J is out of bounds.\n" );
    fprintf ( stderr, "  J = %d\n", j );
    exit ( 1 );
  }

  value = 0;

  if ( i == j )
  {
    return value;
  }

  k = 0;

  while ( k < n )
  {
    if ( a[i+k*m] < a[j+k*m] )
    {
      value = -1;
      return value;
    }
    else if ( a[j+k*m] < a[i+k*m] )
    {
      value = +1;
      return value;
    }
    k = k + 1;
  }

  return value;
}
/******************************************************************************/

double *r8row_indicator_new ( int m, int n )

/******************************************************************************/
/*
  Purpose:

    R8ROW_INDICATOR_NEW sets up an "indicator" R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    The value of each entry suggests its location, as in:

      11  12  13  14
      21  22  23  24
      31  32  33  34

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Output, double R8ROW_INDICATOR_NEW[M*N], the table.
*/
{
  double *a;
  int fac;
  int i;
  int j;

  a = ( double * ) malloc ( m * n * sizeof ( double ) );

  fac = i4_power ( 10, i4_log_10 ( n ) + 1 );

  for ( i = 1; i <= m; i++ )
  {
    for ( j = 1; j <= n; j++ )
    {
      a[i-1+(j-1)*m] = ( double ) ( fac * i + j );
    }
  }
  return a;
}
/******************************************************************************/

double *r8row_max ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_MAX returns the row maximums of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Example:

    A =
      1  2  3
      2  6  7

    MAX =
      3
      7

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 May 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns in the array.

    Input, double A[M*N], the array to be examined.

    Output, double R8ROW_MAX[M], the maximums of the rows.
*/
{
  int i;
  int j;
  double *amax;

  amax = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    amax[i] = a[i+0*m];

    for ( j = 1; j < n; j++ )
    {
      if ( amax[i] < a[i+j*m] )
      {
        amax[i] = a[i+j*m];
      }
    }
  }

  return amax;
}
/******************************************************************************/

double *r8row_mean ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_MEAN returns the row means of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Example:

    A =
      1  2  3
      2  6  7

    MEAN =
      2
      5

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    30 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns in the array.

    Input, double A[M*N], the array to be examined.

    Output, double R8ROW_MEAN[M], the means, or averages, of the rows.
*/
{
  int i;
  int j;
  double *mean;

  mean = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    mean[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      mean[i] = mean[i] + a[i+j*m];
    }
    mean[i] = mean[i] / ( double ) ( n );
  }

  return mean;
}
/******************************************************************************/

double *r8row_min ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_MIN returns the row minimums of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Example:

    A =
      1  2  3
      2  6  7

    MIN =
      1
      2

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 May 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns in the array.

    Input, double A[M*N], the array to be examined.

    Output, double R8ROW_MIN[M], the minimums of the rows.
*/
{
  int i;
  int j;
  double *amin;

  amin = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    amin[i] = a[i+0*m];
    for ( j = 1; j < n; j++ )
    {
      if ( a[i+j*m] < amin[i] )
      {
        amin[i] = a[i+j*m];
      }
    }
  }

  return amin;
}
/******************************************************************************/

void r8row_print ( int m, int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8ROW_PRINT prints an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, double A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  r8row_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void r8row_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    R8ROW_PRINT_SOME prints some of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, double A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:  ");
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %7d     ", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to) 5 entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %14g", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

void r8row_reverse ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_REVERSE reverses the order of the rows of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    To reverse the rows is to start with something like

      11 12 13 14 15
      21 22 23 24 25
      31 32 33 34 35
      41 42 43 44 45
      51 52 53 54 55

    and return

      51 52 53 54 55
      41 42 43 44 45
      31 32 33 34 35
      21 22 23 24 25
      11 12 13 14 15

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    06 May 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input/output, double A[M*N], the matrix whose rows are to be flipped.
*/
{
  int i;
  int j;
  double t;

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < ( m / 2 ); i++ )
    {
      t            = a[    i+j*m];
      a[    i+j*m] = a[m-1-i+j*m];
      a[m-1-i+j*m] = t;
    }
  }
  return;
}
/******************************************************************************/

double *r8row_running_average ( int m, int n, double v[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_RUNNING_AVERAGE computes the running averages of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows.

    Input, int N, the number of items in each row.

    Input, double V[M*N], the data.

    Output, double R8ROW_RUNNING_AVERAGE[M*(N+1)], the running average.  
    A(I,J) is the average value of V(I,1:J-1).
*/
{
  double *a;
  int i;
  int j;

  a = ( double * ) malloc ( m * ( n + 1 ) * sizeof ( double ) );
/*
  Sum.
*/
  for ( i = 0; i < m; i++ )
  {
    a[i+0*m] = 0.0;
    for ( j = 1; j < n + 1; j++ )
    {
      a[i+j*m] = a[i+(j-1)*m] + v[i+(j-1)*m];
    }
  }
/*
  Average.
*/
  for ( i = 0; i < m; i++ )
  {
    for ( j = 1; j < n + 1; j++ )
    {
      a[i+j*m] = a[i+j*m] / ( double ) ( j );
    }
  }

  return a;
}
/******************************************************************************/

double *r8row_running_sum ( int m, int n, double v[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_RUNNING_SUM computes the running sum of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows.

    Input, int N, the number of items in each row.

    Input, double V[M,N], the data.

    Output, double R8ROW_RUNNING_SUM[M*(N+1)], the running sums.  
    S(I,J) is the sum of V(i,1:J-1).
*/
{
  int i;
  int j;
  double *s;

  s = ( double * ) malloc ( m * ( n + 1 ) * sizeof ( double ) );
/*
  Sum.
*/
  for ( i = 0; i < m; i++ )
  {
    s[i+0*m] = 0.0;
    for ( j = 1; j < n + 1; j++ )
    {
      s[i+j*m] = s[i+(j-1)*m] + v[i+(j-1)*m];
    }
  }

  return s;
}
/******************************************************************************/

void r8row_sort_heap_a ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_SORT_HEAP_A ascending heapsorts an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    In lexicographic order, the statement "X < Y", applied to two real
    vectors X and Y of length M, means that there is some index I, with
    1 <= I <= M, with the property that

      X(J) = Y(J) for J < I,
    and
      X(I) < Y(I).

    In other words, the first time they differ, X is smaller.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    25 May 2012

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input/output, double A[M*N].
    On input, the array of N columns of M-vectors.
    On output, the rows of A have been sorted in lexicographic order.
*/
{
  int i;
  int indx;
  int isgn;
  int j;

  if ( m <= 0 )
  {
    return;
  }

  if ( n <= 1 )
  {
    return;
  }
/*
  Initialize.
*/
  i = 0;
  indx = 0;
  isgn = 0;
  j = 0;
/*
  Call the external heap sorter.
*/
  for ( ; ; )
  {
    sort_heap_external ( m, &indx, &i, &j, isgn );
/*
  Interchange the I and J objects.
*/
    if ( 0 < indx )
    {
      r8row_swap ( m, n, a, i - 1, j - 1 );
    }
/*
  Compare the I and J objects.
*/
    else if ( indx < 0 )
    {
      isgn = r8row_compare ( m, n, a, i - 1, j - 1 );
    }
    else if ( indx == 0 )
    {
      break;
    }
  }

  return;
}
/******************************************************************************/

double *r8row_sum ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_SUM returns the sums of the rows of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 May 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A[M*N], the M by N array.

    Output, double ROWSUM[M], the sum of the entries of
    each row.
*/
{
  int i;
  int j;
  double *rowsum;

  rowsum = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    rowsum[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      rowsum[i] = rowsum[i] + a[i+j*m];
    }
  }

  return rowsum;
}
/******************************************************************************/

void r8row_swap ( int m, int n, double a[], int irow1, int irow2 )

/******************************************************************************/
/*
  Purpose:

    R8ROW_SWAP swaps two rows of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    The row indices are 1 based, NOT 0 based.  However, a preprocessor
    variable, called OFFSET, can be reset from 1 to 0 if you wish to
    use 0-based indices.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 September 2003

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input/output, double A[M*N], an array of data.

    Input, int IROW1, IROW2, the two rows to swap.
    These indices should be between 0 and M - 1.
*/
{
  int j;
  double t;
/*
  Check.
*/
  if ( irow1 < 0 || m <= irow1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8ROW_SWAP - Fatal error!\n" );
    fprintf ( stderr, "  IROW1 is out of range.\n" );
    fprintf ( stderr, "  IROW1 = %d\n", irow1 );
    fprintf ( stderr, "  M = %d\n", m );
    exit ( 1 );
  }

  if ( irow2 < 0 || m <= irow2 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8ROW_SWAP - Fatal error!\n" );
    fprintf ( stderr, "  IROW2 is out of range.\n" );
    fprintf ( stderr, "  IROW2 = %d\n", irow2 );
    fprintf ( stderr, "  M = %d\n", m );
    exit ( 1 );
  }

  if ( irow1 == irow2 )
  {
    return;
  }

  for ( j = 0; j < n; j++ )
  {
    t              = a[irow1+j*m];
    a[irow1+j*m] = a[irow2+j*m];
    a[irow2+j*m] = t;
  }

  return;
# undef OFFSET
}
/******************************************************************************/

double *r8row_to_r8vec ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_TO_R8VEC converts an R8ROW into an R8VEC.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    M = 3, N = 4

    A =
      11 12 13 14
      21 22 23 24
      31 32 33 34

    R8ROW_TO_R8VEC = ( 11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34 )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 May 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A[M*N], the M by N array.

    Output, double R8ROW_TO_R8VEC[M*N], a vector containing the M rows of A.
*/
{
  int i;
  int j;
  int k;
  double *x;

  x = ( double * ) malloc ( m * n * sizeof ( double ) );

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      x[k] = a[i+j*m];
      k = k + 1;
    }
  }

  return x;
}
/******************************************************************************/

void r8row_transpose_print ( int m, int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8ROW_TRANSPOSE_PRINT prints an R8ROW, transposed.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A[M*N], an M by N matrix to be printed.

    Input, char *TITLE, a title.
*/
{
  r8row_transpose_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void r8row_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
  int ihi, int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    R8ROW_TRANSPOSE_PRINT_SOME prints some of an R8ROW, transposed.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 February 2016

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A[M*N], an M by N matrix to be printed.

    Input, int ILO, JLO, the first row and column to print.

    Input, int IHI, JHI, the last row and column to print.

    Input, char *TITLE, a title.
*/
{
# define INCX 5

  int i;
  int i2;
  int i2hi;
  int i2lo;
  int i2lo_hi;
  int i2lo_lo;
  int inc;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }

  if ( ilo < 1 )
  {
    i2lo_lo = 1;
  }
  else
  {
    i2lo_lo = ilo;
  }

  if ( ihi < m )
  {
    i2lo_hi = m;
  }
  else
  {
    i2lo_hi = ihi;
  }

  for ( i2lo = i2lo_lo; i2lo <= i2lo_hi; i2lo = i2lo + INCX )
  {
    i2hi = i2lo + INCX - 1;

    if ( m < i2hi )
    {
      i2hi = m;
    }
    if ( ihi < i2hi )
    {
      i2hi = ihi;
    }

    inc = i2hi + 1 - i2lo;

    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row:" );
    for ( i = i2lo; i <= i2hi; i++ )
    {
      fprintf ( stdout, "  %7d     ", i - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Col\n" );
    fprintf ( stdout, "\n" );

    if ( jlo < 1 )
    {
      j2lo = 1;
    }
    else
    {
      j2lo = jlo;
    }
    if ( n < jhi )
    {
      j2hi = n;
    }
    else
    {
      j2hi = jhi;
    }
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "%5d:", j - 1 );
      for ( i2 = 1; i2 <= inc; i2++ )
      {
        i = i2lo - 1 + i2;
        fprintf ( stdout, "  %14g", a[(i-1)+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

double *r8row_uniform_ab_new ( int m, int n, double a, double b, int *seed )

/******************************************************************************/
/*
  Purpose:

    R8ROW_UNIFORM_AB_NEW returns a scaled pseudorandom R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

    This routine implements the recursion

      seed = 16807 * seed mod ( 2^31 - 1 )
      unif = seed / ( 2^31 - 1 )

    The integer arithmetic never requires more than 32 bits,
    including a sign bit.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 October 2005

  Author:

    John Burkardt

  Reference:

    Paul Bratley, Bennett Fox, Linus Schrage,
    A Guide to Simulation,
    Second Edition,
    Springer, 1987,
    ISBN: 0387964673,
    LC: QA76.9.C65.B73.

    Bennett Fox,
    Algorithm 647:
    Implementation and Relative Efficiency of Quasirandom
    Sequence Generators,
    ACM Transactions on Mathematical Software,
    Volume 12, Number 4, December 1986, pages 362-376.

    Pierre L'Ecuyer,
    Random Number Generation,
    in Handbook of Simulation,
    edited by Jerry Banks,
    Wiley, 1998,
    ISBN: 0471134031,
    LC: T57.62.H37.

    Peter Lewis, Allen Goodman, James Miller,
    A Pseudo-Random Number Generator for the System/360,
    IBM Systems Journal,
    Volume 8, Number 2, 1969, pages 136-143.

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A, B, the limits of the pseudorandom values.

    Input/output, int *SEED, the "seed" value.  Normally, this
    value should not be 0.  On output, SEED has 
    been updated.

    Output, double R8ROW_UNIFORM_AB_NEW[M*N], a matrix of pseudorandom values.
*/
{
  int i;
  const int i4_huge = 2147483647;
  int j;
  int k;
  double *r;

  if ( *seed == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8ROW_UNIFORM_AB_NEW - Fatal error!\n" );
    fprintf ( stderr, "  Input value of SEED = 0.\n" );
    exit ( 1 );
  }

  r = ( double * ) malloc ( m * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      k = *seed / 127773;

      *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

      if ( *seed < 0 )
      {
        *seed = *seed + i4_huge;
      }
      r[i+j*m] = a + ( b - a ) * ( double ) ( *seed ) * 4.656612875E-10;
    }
  }

  return r;
}
/******************************************************************************/

double *r8row_variance ( int m, int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8ROW_VARIANCE returns the variances of an R8ROW.

  Discussion:

    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
    each of length N.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 August 2011

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns in the array.

    Input, double A[M*N], the array whose variances are desired.

    Output, double R8ROW_VARIANCE[M], the variances of the rows.
*/
{
  int i;
  int j;
  double mean;
  double *variance;

  variance = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    mean = 0.0;
    for ( j = 0; j < n; j++ )
    {
      mean = mean + a[i+j*m];
    }
    mean = mean / ( double ) ( n );

    variance[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      variance[i] = variance[i] + pow ( ( a[i+j*m] - mean ), 2 );
    }

    if ( 1 < n )
    {
      variance[i] = variance[i] / ( double ) ( n - 1 );
    }
    else
    {
      variance[i] = 0.0;
    }

  }

  return variance;
}
/******************************************************************************/

int r8vec_eq ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_EQ is true if two R8VEC's are equal.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 August 2003

  Author:

    John Burkardt

  Parameters:

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

    Input, double A1[N], A2[N], two vectors to compare.

    Output, int R8VEC_EQ, is TRUE if every pair of elements A1(I) and A2(I) are equal,
    and FALSE otherwise.
*/
{
  int i;
  int value;

  for ( i = 0; i < n; i++ )
  {
    if ( a1[i] != a2[i] )
    {
      value = 0;
      return value;
    }
  }
  value = 1;
  return value;
}
/******************************************************************************/

int r8vec_gt ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_GT == ( A1 > A2 ) for two R8VEC's.

  Discussion:

    An R8VEC is a vector of R8's.

    The comparison is lexicographic.

    A1 > A2  <=>                              A1(1) > A2(1) or
                 ( A1(1)     == A2(1)     and A1(2) > A2(2) ) or
                 ...
                 ( A1(1:N-1) == A2(1:N-1) and A1(N) > A2(N)

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 January 2010

  Author:

    John Burkardt

  Parameters:

    Input, int N, the dimension of the vectors.

    Input, double A1[N], A2[N], the vectors to be compared.

    Output, int R8VEC_GT, is TRUE if and only if A1 > A2.
*/
{
  int i;

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

    if ( a2[i] < a1[i] )
    {
       return 1;
    }
    else if ( a1[i] < a2[i] )
    {
      return 0;
    }

  }

  return 0;
}
/******************************************************************************/

int r8vec_lt ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_LT == ( A1 < A2 ) for two R8VEC's.

  Discussion:

    An R8VEC is a vector of R8's.

    The comparison is lexicographic.

    A1 < A2  <=>                              A1(1) < A2(1) or
                 ( A1(1)     == A2(1)     and A1(2) < A2(2) ) or
                 ...
                 ( A1(1:N-1) == A2(1:N-1) and A1(N) < A2(N)

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 January 2010

  Author:

    John Burkardt

  Parameters:

    Input, int N, the dimension of the vectors.

    Input, double A1[N], A2[N], the vectors to be compared.

    Output, int R8VEC_LT, is TRUE if and only if A1 < A2.
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    if ( a1[i] < a2[i] )
    {
      return 1;
    }
    else if ( a2[i] < a1[i] )
    {
      return 0;
    }

  }

  return 0;
}
/******************************************************************************/

void r8vec_print ( int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8VEC_PRINT prints an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 April 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, double A[N], the vector to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );
  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %8d: %14g\n", i, a[i] );
  }

  return;
}
/******************************************************************************/

void sort_heap_external ( int n, int *indx, int *i, int *j, int isgn )

/******************************************************************************/
/*
  Purpose:

    SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.

  Discussion:

    The actual list is not passed to the routine.  Hence it may
    consist of integers, reals, numbers, names, etc.  The user,
    after each return from the routine, will be asked to compare or
    interchange two items.

    The current version of this code mimics the FORTRAN version,
    so the values of I and J, in particular, are FORTRAN indices.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    05 February 2004

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int N, the length of the input list.

    Input/output, int *INDX.
    The user must set INDX to 0 before the first call.
    On return,
      if INDX is greater than 0, the user must interchange
      items I and J and recall the routine.
      If INDX is less than 0, the user is to compare items I
      and J and return in ISGN a negative value if I is to
      precede J, and a positive value otherwise.
      If INDX is 0, the sorting is done.

    Output, int *I, *J.  On return with INDX positive,
    elements I and J of the user's list should be
    interchanged.  On return with INDX negative, elements I
    and J are to be compared by the user.

    Input, int ISGN. On return with INDX negative, the
    user should compare elements I and J of the list.  If
    item I is to precede item J, set ISGN negative,
    otherwise set ISGN positive.
*/
{
  static int i_save = 0;
  static int j_save = 0;
  static int k = 0;
  static int k1 = 0;
  static int n1 = 0;
/*
  INDX = 0: This is the first call.
*/
  if ( *indx == 0 )
  {

    i_save = 0;
    j_save = 0;
    k = n / 2;
    k1 = k;
    n1 = n;
  }
/*
  INDX < 0: The user is returning the results of a comparison.
*/
  else if ( *indx < 0 )
  {
    if ( *indx == -2 )
    {
      if ( isgn < 0 )
      {
        i_save = i_save + 1;
      }
      j_save = k1;
      k1 = i_save;
      *indx = -1;
      *i = i_save;
      *j = j_save;
      return;
    }

    if ( 0 < isgn )
    {
      *indx = 2;
      *i = i_save;
      *j = j_save;
      return;
    }

    if ( k <= 1 )
    {
      if ( n1 == 1 )
      {
        i_save = 0;
        j_save = 0;
        *indx = 0;
      }
      else
      {
        i_save = n1;
        j_save = 1;
        n1 = n1 - 1;
        *indx = 1;
      }
      *i = i_save;
      *j = j_save;
      return;
    }
    k = k - 1;
    k1 = k;
  }
/*
  0 < INDX: the user was asked to make an interchange.
*/
  else if ( *indx == 1 )
  {
    k1 = k;
  }

  for ( ; ; )
  {

    i_save = 2 * k1;

    if ( i_save == n1 )
    {
      j_save = k1;
      k1 = i_save;
      *indx = -1;
      *i = i_save;
      *j = j_save;
      return;
    }
    else if ( i_save <= n1 )
    {
      j_save = i_save + 1;
      *indx = -2;
      *i = i_save;
      *j = j_save;
      return;
    }

    if ( k <= 1 )
    {
      break;
    }

    k = k - 1;
    k1 = k;
  }

  if ( n1 == 1 )
  {
    i_save = 0;
    j_save = 0;
    *indx = 0;
    *i = i_save;
    *j = j_save;
  }
  else
  {
    i_save = n1;
    j_save = 1;
    n1 = n1 - 1;
    *indx = 1;
    *i = i_save;
    *j = j_save;
  }

  return;
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}