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

# include "r8utp.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:

    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

  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

  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;
}
/******************************************************************************/

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

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

    r8ge_print() prints an R8GE matrix.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 February 2012

  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 R8GE matrix.

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

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

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

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

    r8ge_print_some() prints some of an R8GE matrix.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 February 2012

  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 R8GE 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;

  printf ( "\n" );
  printf ( "%s\n", title );
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

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

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

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

double *r8ge_to_r8utp ( int m, int n, double a_ge[] )

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

    r8ge_to_r8utp() copies an R8GE matrix to an R8UTP matrix.

  Discussion:

    The R8GE storage format is used for a general M by N matrix.  A storage 
    space is made for each entry.  The two dimensional logical
    array can be thought of as a vector of M*N entries, starting with
    the M entries in the column 1, then the M entries in column 2
    and so on.  Considered as a vector, the entry A(I,J) is then stored
    in vector location I+(J-1)*M.

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    15 August 2022

  Author:

    John Burkardt

  Input:

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

    double A_GE(N,N), the R8GE matrix.

  Output:

    double r8ge_to_r8utp[*], the R8UTP matrix.
*/
{
  double *a_utp;
  int i;
  int j;
  int k;

  a_utp = r8utp_zeros ( m, n );

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < i4_min ( j + 1, m ); i++ )
    {
      a_utp[k] = a_ge[i+j*m];
      k = k + 1;
    }
  }

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

double r8utp_det ( int m, int n, double a[] )

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

    r8utp_det() computes the determinant of an R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    14 August 2022

  Author:

    John Burkardt

  Input:

    int M, N, the order of the matrix.

    double A[*], the R8UTP matrix.

  Output:

    double r8utp_det: the determinant of the matrix.
*/
{
  double det;
  int j;
  int k;

  if ( m != n )
  {
    printf ( "\n" );
    printf ( "r8utp_det(): Fatal error!\n" );
    printf ( "  m and n are not equal.\n" );
    exit ( 1 );
  }

  det = 1.0;
  k = 0;
  for ( j = 0; j < n; j++ )
  {
    det = det * a[k];
    k = k + j + 2;
  }

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

double *r8utp_indicator ( int m, int n )

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

    r8utp_indicator() sets up a R8UTP indicator matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 August 2022

  Author:

    John Burkardt

  Input:

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

  Output:

    double r8utp_indicator[]: the R8UTP matrix.
*/
{
  double *a;
  int fac;
  int i;
  int j;
  int k;
  int mn;

  mn = r8utp_size ( m, n );
  a = ( double * ) malloc ( mn * sizeof ( double ) );
  fac = pow ( 10, i4_log_10 ( n ) + 1 );

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

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

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

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

    r8utp_print() prints a R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN). 

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 April 2014

  Author:

    John Burkardt

  Parameters:

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

    Input, double A[*], the matrix.

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

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

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

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

    r8utp_print_some() prints some of a R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 April 2014

  Author:

    John Burkardt

  Parameters:

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

    Input, double A[*], 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

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

  printf ( "\n" );
  printf ( "%s\n", title );
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

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

    for ( i = i2lo; i <= i2hi; i++ )
    {
      printf ( "%6d  ", i );
/*
  Print out (up to) 5 entries in row I, that lie in the current strip.
*/
      for ( j = j2lo; j <= j2hi; j++ )
      {
        if ( i <= j )
        {
          aij = a[i-1+(j*(j-1))/2];
        }
        else
        {
          aij = 0.0;
        }

        printf ( "%12g  ", aij );
      }
      printf ( "\n" );
    }
  }

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

double *r8utp_random ( int m, int n )

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

    r8utp_random() randomizes an R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 August 2022

  Author:

    John Burkardt

  Input:

    int M, N, the number of rows and columns of the matrix.
    M and N must be positive.

  Output:

    double r8utp_random[*], the R8UTP matrix.
*/
{
  double *a;
  int i;
  int mn;

  mn = r8utp_size ( m, n );
  a = ( double * ) malloc ( mn * sizeof ( double ) );

  for ( i = 0; i < mn; i++ )
  {
    a[i] = drand48 ( );
  }

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

int r8utp_size ( int m, int n )

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

    r8utp_size() returns the size of an M x N R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 August 2022

  Author:

    John Burkardt

  Input:

    integer M, N: the number of rows and columns.

  Output:

    integer r8utp_size: the length of the array needed to store the matrix.
*/
{
  int value;

  if ( n < m )
  {
    value = ( n * ( n + 1 ) ) / 2;
  }
  else if ( m == n )
  {
    value = ( m * ( m + 1 ) ) / 2;
  }
  else
  {
    value = ( m * ( m + 1 ) / 2 ) + ( n - m ) * m;
  }

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

double *r8utp_to_r8ge ( int m, int n, double a_utp[] )

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

    r8utp_to_r8ge() copies an R8UTP matrix to an R8GE matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

    The R8GE storage format is used for a general M by N matrix.  A storage 
    space is made for each entry.  The two dimensional logical
    array can be thought of as a vector of M*N entries, starting with
    the M entries in the column 1, then the M entries in column 2
    and so on.  Considered as a vector, the entry A(I,J) is then stored
    in vector location I+(J-1)*M.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 August 2022

  Author:

    John Burkardt

  Input:

    int M, N, the order of the matrix.

    double A_UTP(*), the R8UTP matrix.

  Output:

    double r8utp_to_r8ge[m,n]: the R8GE matrix.
*/
{
  double *a_ge;
  int i;
  int j;
  int k;

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

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      if ( i < i4_min ( j + 1, m ) )
      {
        a_ge[i+j*m] = a_utp[k];
        k = k + 1;
      }
      else
      {
        a_ge[i+j*m] = 0.0;
      }
    }
  }

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

double *r8utp_zeros ( int m, int n )

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

    r8utp_zeros() zeros an R8UTP matrix.

  Discussion:

    The R8UTP storage format is appropriate for an upper triangular
    matrix.  Only the upper triangle of the matrix is stored,
    by successive partial columns, in an array which contains
    (A11,A12,A22,A13,A23,A33,A14,...,AMN).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 August 2022

  Author:

    John Burkardt

  Input:

    int M, N, the number of rows and columns of the matrix.
    M and N must be positive.

  Output:

    double r8utp_random[*], the R8UTP matrix.
*/
{
  double *a;
  int i;
  int mn;

  mn = r8utp_size ( m, n );
  a = ( double * ) malloc ( mn * sizeof ( double ) );

  for ( i = 0; i < mn; i++ )
  {
    a[i] = 0.0;
  }

  return a;
}