#! /usr/bin/env python3
#
def r8pbl_test ( ):

#*****************************************************************************80
#
## r8pbl_test() tests r8pbl().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    25 August 2022
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'r8pbl_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  Test r8pbl().' )

  r8pbl_dif2_test ( )
  r8pbl_indicator_test ( )
  r8pbl_mv_test ( )
  r8pbl_print_test ( )
  r8pbl_print_some_test ( )
  r8pbl_random_test ( )
  r8pbl_to_r8ge_test ( )
  r8pbl_zeros_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'r8pbl_test():' )
  print ( '  Normal end of execution.' )
  return

def i4_log_10 ( i ):

#*****************************************************************************80
#
## i4_log_10() returns the integer part of the logarithm base 10 of ABS(X).
#
#  Example:
#
#        I  VALUE
#    -----  --------
#        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:
#
#    08 May 2013
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer I, the number whose logarithm base 10 is desired.
#
#  Output:
#
#    integer VALUE, the integer part of the logarithm base 10 of
#    the absolute value of X.
#
  import numpy as np

  i = np.floor ( i )

  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

def r8pbl_dif2 ( n, ml ):

#*****************************************************************************80
#
## R8PBL_DIF2 returns the DIF2 matrix in R8PBL format.
#
#  Example:
#
#    N = 5
#
#    2 -1  .  .  .
#   -1  2 -1  .  .
#    . -1  2 -1  .
#    .  . -1  2 -1
#    .  .  . -1  2
#
#  Properties:
#
#    A is banded, with bandwidth 3.
#
#    A is tridiagonal.
#
#    Because A is tridiagonal, it has property A (bipartite).
#
#    A is a special case of the TRIS or tridiagonal scalar matrix.
#
#    A is integral, therefore det ( A ) is integral, and 
#    det ( A ) * inverse ( A ) is integral.
#
#    A is Toeplitz: constant along diagonals.
#
#    A is symmetric: A' = A.
#
#    Because A is symmetric, it is normal.
#
#    Because A is normal, it is diagonalizable.
#
#    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
#
#    A is positive definite.
#
#    A is an M matrix.
#
#    A is weakly diagonally dominant, but not strictly diagonally dominant.
#
#    A has an LU factorization A = L * U, without pivoting.
#
#      The matrix L is lower bidiagonal with subdiagonal elements:
#
#        L(I+1,I) = -I/(I+1)
#
#      The matrix U is upper bidiagonal, with diagonal elements
#
#        U(I,I) = (I+1)/I
#
#      and superdiagonal elements which are all -1.
#
#    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
#
#      L(I,I) =    sqrt ( (I+1) / I )
#      L(I,I-1) = -sqrt ( (I-1) / I )
#
#    The eigenvalues are
#
#      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
#                = 4 SIN^2(I*PI/(2*N+2))
#
#    The corresponding eigenvector X(I) has entries
#
#       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
#
#    Simple linear systems:
#
#      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
#
#      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
#
#    det ( A ) = N + 1.
#
#    The value of the determinant can be seen by induction,
#    and expanding the determinant across the first row:
#
#      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
#                = 2 * N - (N-1)
#                = N + 1
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 July 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Robert Gregory, David Karney,
#    A Collection of Matrices for Testing Computational Algorithms,
#    Wiley, 1969,
#    ISBN: 0882756494,
#    LC: QA263.68
#
#    Morris Newman, John Todd,
#    Example A8,
#    The evaluation of matrix inversion programs,
#    Journal of the Society for Industrial and Applied Mathematics,
#    Volume 6, Number 4, pages 466-476, 1958.
#
#    John Todd,
#    Basic Numerical Mathematics,
#    Volume 2: Numerical Algebra,
#    Birkhauser, 1980,
#    ISBN: 0817608117,
#    LC: QA297.T58.
#
#    Joan Westlake,
#    A Handbook of Numerical Matrix Inversion and Solution of 
#    Linear Equations,
#    John Wiley, 1968,
#    ISBN13: 978-0471936756,
#    LC: QA263.W47.
#
#  Parameters:
#
#    Input, integer N, the number of rows and columns.
#
#    Input, integer ML, the number of subdiagonals.
#    ML must be at least 0, and no more than N-1.
#
#    Output, real A(ML+1,N), the matrix.
#
  import numpy as np

  a = np.zeros ( [ ml + 1, n ] )

  a[0,0:n]   = +2.0
  a[1,0:n-1] = -1.0
 
  return a

def r8pbl_dif2_test ( ):

#*****************************************************************************80
#
## R8PBL_DIF2_TEST tests R8PBL_DIF2.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 5
  ml = 1

  print ( '' )
  print ( 'R8PBL_DIF2_TEST' )
  print ( '  R8PBL_DIF2 sets an R8PBL second difference matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_dif2 ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL second difference matrix:' )

  return

def r8pbl_indicator ( n, ml ):

#*****************************************************************************80
#
## R8PBL_INDICATOR sets up a R8PBL indicator matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 February 2004
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer ML, the number of subdiagonals in the matrix.
#    ML must be at least 0 and no more than N-1.
#
#    Output, real A(ML+1,N), the R8PBL matrix.
#
  import numpy as np

  a = np.zeros ( [ ml + 1, n ] )

  fac = 10 ** ( i4_log_10 ( n ) + 1 )

  for i in range ( 0, n ):
    for j in range ( max ( 0, i - ml ), i + 1 ):
      a[i-j,j] = float ( fac * ( i + 1 ) + ( j + 1 ) )

  return a

def r8pbl_indicator_test ( ):

#*****************************************************************************80
#
## R8PBL_INDICATOR_TEST tests R8PBL_INDICATOR.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    14 February 2003
#
#  Author:
#
#    John Burkardt
#
  n = 9
  ml = 3

  print ( '' )
  print ( 'R8PBL_INDICATOR_TEST' )
  print ( '  R8PBL_INDICATOR sets up a R8PBL indicator matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_indicator ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL indicator matrix:' )

  return

def r8pbl_mv ( n, ml, a, x ):

#*****************************************************************************80
#
## R8PBL_MV multiplies an R8PBL matrix by an R8VEC.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 July 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#
#    Input, integer ML, the number of subdiagonals in the matrix.
#    ML must be at least 0 and no more than N-1.
#
#    Input, real A(ML+1,N), the R8PBL matrix.
#
#    Input, real X(N), the vector to be multiplied by A.
#
#    Output, real B(N), the result vector A * x.
#
  import numpy as np

  b = np.zeros ( n )
#
#  Multiply X by the diagonal of the matrix.
#
  for i in range ( 0, n ):
    b[i] = a[0,i] * x[i]
#
#  Multiply X by the subdiagonals of the matrix.
#
  for k in range ( 0, ml ):
    for j in range ( 0, n - k ):
      i = j + k
      aij = a[k+1,j]
      b[i] = b[i] + aij * x[j]
      b[j] = b[j] + aij * x[i]

  return b

def r8pbl_mv_test ( ):

#*****************************************************************************80
#
## R8PBL_MV_TEST tests R8PBL_MV.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 5
  ml = 2

  print ( '' )
  print ( 'R8PBL_MV_TEST' )
  print ( '  R8PBL_MV computes A*x, where A is an R8PBL matrix.' )
  print ( '' )
  print ( '  Matrix order N =     ', n )
  print ( '  Lower bandwidth ML = ', ml )
#
#  Set the matrix.
#
  a = r8pbl_random ( n, ml )

  r8pbl_print ( n, ml, a, '  Matrix A:' )
#
#  Set the desired solution.
#
  x = r8vec_indicator1 ( n )

  r8vec_print ( n, x, '  Vector x:' )
#
#  Compute the corresponding right hand side.
#
  b = r8pbl_mv ( n, ml, a, x )

  r8vec_print ( n, b, '  Product b=A*x' )

  return

def r8pbl_print ( n, ml, a, title ):

#*****************************************************************************80
#
## R8PBL_PRINT prints a R8PBL matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 July 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer ML, the upper (and lower) bandwidth.
#    ML must be nonnegative, and no greater than N-1.
#
#    Input, real A(ML+1,N), the R8PBL matrix.
#
#    Input, string TITLE, a title to be printed.
#
  r8pbl_print_some ( n, ml, a, 0, 0, n - 1, n - 1, title )

  return

def r8pbl_print_test ( ):

#*****************************************************************************80
#
## R8PBL_PRINT_TEST tests R8PBL_PRINT.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 9
  ml = 3

  print ( '' )
  print ( 'R8PBL_PRINT_TEST' )
  print ( '  R8PBL_PRINT prints an R8PBL matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_indicator ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL matrix:' )

  return

def r8pbl_print_some ( n, ml, a, ilo, jlo, ihi, jhi, title ):

#*****************************************************************************80
#
## R8PBL_PRINT_SOME prints some of a R8PBL matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 July 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer ML, the upper (and lower) bandwidth.
#    ML must be nonnegative, and no greater than N-1.
#
#    Input, real A(ML+1,N), the R8PBL matrix.
#
#    Input, integer ILO, JLO, IHI, JHI, the first row and
#    column, and the last row and column to be printed.
#
#    Input, string TITLE, a title.
#
  incx = 5

  print ( '' )
  print ( title )

  for j2lo in range ( max ( jlo, 0 ), min ( jhi + 1, n ), incx ):

    j2hi = j2lo + incx - 1
    j2hi = min ( j2hi, n - 1 )
    j2hi = min ( j2hi, jhi )
    
    print ( '' )
    print ( '  Col: ', end = '' )

    for j in range ( j2lo, j2hi + 1 ):
      print ( '%7d       ' % ( j ), end = '' )

    print ( '' )
    print ( '  Row' )

    i2lo = max ( ilo, 0 )
    i2lo = max ( i2lo, j2lo - ml )
    i2hi = min ( ihi, n - 1 )
    i2hi = min ( i2hi, j2hi + ml )

    for i in range ( i2lo, i2hi + 1 ):

      print ( '%4d :' % ( i ), end = '' )
      
      for j in range ( j2lo, j2hi + 1 ):

        if ( j < i - ml or i + ml < j ):
          print ( '              ', end = '' )
        elif ( i <= j and j <= i + ml ):
          print ( '%12g  ' % ( a[j-i,i] ), end = '' )
        elif ( j <= i and i <= j + ml ):
          print ( '%12g  ' % ( a[i-j,j] ), end = '' )
 
      print ( '' )

  return

def r8pbl_print_some_test ( ):

#*****************************************************************************80
#
## R8PBL_PRINT_SOME_TEST tests R8PBL_PRINT_SOME.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 9
  ml = 4

  print ( '' )
  print ( 'R8PBL_PRINT_SOME_TEST' )
  print ( '  R8PBL_PRINT_SOME prints some of an R8PBL matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_indicator ( n, ml )

  r8pbl_print_some ( n, ml, a, 3, 4, 7, 8, '  Row(3:7), Col(4:8):' )

  return

def r8pbl_random ( n, ml ):

#*****************************************************************************80
#
## R8PBL_RANDOM randomizes a R8PBL matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#    The matrix returned will be positive definite, but of limited
#    randomness.  The off diagonal elements are random values between
#    0 and 1, and the diagonal element of each row is selected to
#    ensure strict diagonal dominance.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    16 February 2005
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer ML, the number of subdiagonals in the matrix.
#    ML must be at least 0 and no more than N-1.
#
#    Output, real A(ML+1,N), the R8PBL matrix.
#
  from numpy.random import default_rng
  import numpy as np

  rng = default_rng ( )

  a = np.zeros ( [ ml + 1, n ] )
#
#  Set the off diagonal values.
#
  for i in range ( 0, n ):
    for j in range ( max ( 0, i - ml ), i ):
      a[i-j,j] = rng.random ( size = 1 )
#
#  Set the diagonal values.
#
  for i in range ( 0, n ):

    sum2 = 0.0

    jlo = max ( 0, i - ml )
    for j in range ( jlo, i ):
      sum2 = sum2 + abs ( a[i-j,j] )

    jhi = min ( i + ml, n - 1 )
    for j in range ( i + 1, jhi + 1 ):
      sum2 = sum2 + abs ( a[j-i,i] )

    r = rng.random ( size = 1 )

    a[0,i] = ( 1.0 + r ) * ( sum2 + 0.01 )

  return a

def r8pbl_random_test ( ):

#*****************************************************************************80
#
## R8PBL_RANDOM_TEST tests R8PBL_RANDOM.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 5
  ml = 2

  print ( '' )
  print ( 'R8PBL_RANDOM_TEST' )
  print ( '  R8PBL_RANDOM randomizes an R8PBL matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_random ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL random matrix:' )

  return

def r8pbl_to_r8ge ( n, ml, a ):

#*****************************************************************************80
#
## R8PBL_TO_R8GE copies a R8PBL matrix to a R8GE matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 July 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrices.
#    N must be positive.
#
#    Input, integer ML, the upper bandwidth of A.
#    ML must be nonnegative, and no greater than N-1.
#
#    Input, real A(ML+1,N), the R8PBL matrix.
#
#    Output, real B(N,N), the R8GE matrix.
#
  import numpy as np

  b = np.zeros ( [ n, n ] )

  for i in range ( 0, n ):
    for j in range ( 0, n ):
      if ( i <= j and j <= i + ml ):
        b[i,j] = a[j-i,i]
      elif ( i - ml <= j and j < i ):
        b[i,j] = a[i-j,j]

  return b

def r8pbl_to_r8ge_test ( ):

#*****************************************************************************80
#
## R8PBL_TO_R8GE_TEST tests R8PBL_TO_R8GE.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 5
  ml = 2

  print ( '' )
  print ( 'R8PBL_TO_R8GE_TEST' )
  print ( '  R8PBL_TO_R8GE converts an R8PBL matrix to R8GE format.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_random ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL matrix:' )

  a_r8ge = r8pbl_to_r8ge ( n, ml, a )

  print ( '' )
  print ( '  The R8GE matrix:' )
  print ( a_r8ge )

  return

def r8pbl_zeros ( n, ml ):

#*****************************************************************************80
#
## R8PBL_ZEROS zeros an R8PBL matrix.
#
#  Discussion:
#
#    The R8PBL storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and lower triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row 1 of the array.
#    The first subdiagonal in row 2, columns 1 through ML.
#    The second subdiagonal in row 3, columns 1 through ML-1.
#    The ML-th subdiagonal in row ML+1, columns 1 through 1.
#
#    The matrix returned will be positive definite, but of limited
#    randomness.  The off diagonal elements are random values between
#    0 and 1, and the diagonal element of each row is selected to
#    ensure strict diagonal dominance.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 August 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer ML, the number of subdiagonals in the matrix.
#    ML must be at least 0 and no more than N-1.
#
#    Output, real A(ML+1,N), the R8PBL matrix.
#
  import numpy as np

  a = np.zeros ( [ ml + 1, n ] )

  return a

def r8pbl_zeros_test ( ):

#*****************************************************************************80
#
## R8PBL_ZEROS_TEST tests R8PBL_ZEROS.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 July 2016
#
#  Author:
#
#    John Burkardt
#
  n = 5
  ml = 2

  print ( '' )
  print ( 'R8PBL_ZEROS_TEST' )
  print ( '  R8PBL_ZEROS zeros an R8PBL matrix.' )
  print ( '' )
  print ( '  Matrix order N = ', n )
  print ( '  Bandwidth ML = ', ml )

  a = r8pbl_zeros ( n, ml )

  r8pbl_print ( n, ml, a, '  The R8PBL zero matrix:' )

def r8vec_indicator1 ( n ):

#*****************************************************************************80
#
## r8vec_indicator1() sets an R8VEC to the indicator vector (1,2,3,...).
#
#  Discussion:
#
#    An R8VEC is a vector of R8's.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 September 2014
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of elements of the vector.
#
#  Output:
#
#    real A(N), the indicator array.
#
  import numpy as np

  a = np.zeros ( n );

  for i in range ( 0, n ):
    a[i] = i + 1

  return a

def r8vec_print ( n, a, title ):

#*****************************************************************************80
#
## r8vec_print() prints an R8VEC.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2014
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the dimension of the vector.
#
#    real A(N), the vector to be printed.
#
#    string TITLE, a title.
#
  print ( '' )
  print ( title )
  print ( '' )
  for i in range ( 0, n ):
    print ( '%6d:  %12g' % ( i, a[i] ) )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  r8pbl_test ( )
  timestamp ( )