#! /usr/bin/env python3 # def r8ci_test ( ): #*****************************************************************************80 # ## r8ci_test() tests r8ci(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 August 2022 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'r8ci_test():' ) print ( ' Python version: ' + platform.python_version ( ) ) print ( ' Test r8ci().' ) r8ci_det_test ( ) r8ci_dif2_test ( ) r8ci_eval_test ( ) r8ci_indicator_test ( ) r8ci_mtv_test ( ) r8ci_mv_test ( ) r8ci_print_test ( ) r8ci_print_some_test ( ) r8ci_random_test ( ) r8ci_sl_test ( ) r8ci_to_r8ge_test ( ) r8ci_zeros_test ( ) # # Terminate. # print ( '' ) print ( 'r8ci_test():' ) print ( ' Normal end of execution.' ) return def c8_le_l2 ( c1, c2 ): #*****************************************************************************80 # ## c8_le_l2() := C1 <= C1 for C8's, and the L2 norm. # # Discussion: # # The L2 norm can be defined here as: # # c8_norm_l2(C) = sqrt ( ( real (C) ) ^ 2 + abs ( imag (C) ) ^ 2 ) # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 19 February 2015 # # Author: # # John Burkardt # # Input: # # complex C1, C2, the values to be compared. # # Output: # # bool VALUE, is TRUE if C1 <= C2. # if ( c1.real ** 2 + c1.imag ** 2 <= c2.real ** 2 + c2.imag ** 2 ) : value = True else: value = False return value def c8vec_print ( n, a, title ): #*****************************************************************************80 # ## c8vec_print() prints a C8VEC. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 August 2022 # # Author: # # John Burkardt # # Input: # # integer N, the dimension of the vector. # # complex A(N), the vector to be printed. # # string TITLE, a title. # print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( '%6d %12g %12g' % ( i, a.real[i], a.imag[i] ) ) return def c8vec_sort_a_l2 ( n, x ): #*****************************************************************************80 # ## c8vec_sort_a_l2() ascending sorts a C8VEC by L2 norm. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 28 March 2015 # # Author: # # John Burkardt # # Input: # # integer N, the number of entries in the array. # # complex X(N), an unsorted array. # # Output: # # complex X(N), the sorted array. # if ( 1 < n ): i = 0 indx = 0 isgn = 0 j = 0 i_save = 0 j_save = 0 k_save = 0 l_save = 0 n_save = 0 while ( True ): [ indx, i, j, i_save, j_save, k_save, l_save, n_save ] = sort_safe_rc ( \ n, indx, isgn, i_save, j_save, k_save, l_save, n_save ) if ( 0 < indx ): temp = x[i-1] x[i-1] = x[j-1] x[j-1] = temp elif ( indx < 0 ): if ( c8_le_l2 ( x[i-1], x[j-1] ) ): isgn = -1 else: isgn = +1 elif ( indx == 0 ): break return x def c8vec_unity ( n ): #*****************************************************************************80 # ## c8vec_unity() returns the N roots of unity. # # Discussion: # # X(1:N) = exp ( 2 * PI * (0:N-1) / N ) # # X(1:N)^N = ( (1,0), (1,0), ..., (1,0) ). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 13 February 2015 # # Author: # # John Burkardt # # Input: # # integer N, the number of elements. # # Output: # # complex A(N), the array. # import numpy as np a = np.zeros ( n, 'complex' ) for i in range ( 0, n ): t = 2.0 * np.pi * float ( i ) / float ( n ) a[i] = np.cos ( t ) + 1j * np.sin ( t ) return a 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 r8ci_det ( n, a ): #*****************************************************************************80 # ## R8CI_DET returns the determinant of an R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI matrix. # # Output, real DET, the determinant. # import numpy as np w = c8vec_unity ( n ) lam = np.zeros ( n, 'complex' ) for i in range ( 0, n ): lam[i] = a[n-1] for i in range ( n - 2, -1, -1 ): for j in range ( 0, n ): lam[j] = lam[j] * w[j] + a[i] detc = 1.0 + 0j for i in range ( 0, n ): detc = detc * lam[i] # # DETC should be real, but we need to copy it into a real variable. # det = detc.real return det def r8ci_det_test ( ): #*****************************************************************************80 # ## R8CI_DET_TEST tests R8CI_DET. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_DET_TEST' ) print ( ' R8CI_DET finds the determinant of a real circulant system.' ) print ( '' ) print ( ' Matrix order N = ', n ) # # Set the matrix. # a = r8ci_random ( n ) r8ci_print ( n, a, ' The circulant matrix:' ) det = r8ci_det ( n, a ) print ( '' ) print ( ' Computed determinant = ', det ) return def r8ci_dif2 ( n ): #*****************************************************************************80 # ## R8CI_DIF2 sets up a R8CI second difference matrix. # # Discussion: # # This is actually a periodic second difference matrix. # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Output, real A(N), the R8CI matrix. # import numpy as np a = np.zeros ( n ) a[0] = 2.0 a[1] = -1.0 a[n-1] = -1.0 return a def r8ci_dif2_test ( ): #*****************************************************************************80 # ## R8CI_DIF2_TEST tests R8CI_DIF2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_DIF2_TEST' ) print ( ' R8CI_DIF2 sets up an R8CI periodic second difference matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_dif2 ( n ) r8ci_print ( n, a, ' The R8CI periodic second difference matrix:' ) return def r8ci_eval ( n, a ): #*****************************************************************************80 # ## R8CI_EVAL returns the eigenvalues of a R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Reference: # # Philip Davis, # Circulant Matrices, # Wiley, 1979. # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI matrix. # # Output, complex LAMBDA(N), the eigenvalues. # import numpy as np w = c8vec_unity ( n ) lam = np.zeros ( n, 'complex' ) for i in range ( 0, n ): lam[i] = a[n-1] for i in range ( n - 2, -1, -1 ): for j in range ( 0, n ): lam[j] = lam[j] * w[j] + a[i] lam = c8vec_sort_a_l2 ( n, lam ) return lam def r8ci_eval_test ( ): #*****************************************************************************80 # ## R8CI_EVAL_TEST tests R8CI_EVAL. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_EVAL_TEST' ) print ( ' R8CI_EVAL finds the eigenvalues of' ) print ( ' a real circulant system.' ) print ( '' ) print ( ' Matrix order N = ', n ) # # Set the matrix. # a = r8ci_random ( n ) r8ci_print ( n, a, ' The circulant matrix:' ) lam = r8ci_eval ( n, a ) c8vec_print ( n, lam, ' The eigenvalues:' ) return def r8ci_indicator ( n ): #*****************************************************************************80 # ## R8CI_INDICATOR sets up a R8CI indicator matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Output, real A(N), the R8CI matrix. # import numpy as np a = np.zeros ( n ) fac = 10 ** ( i4_log_10 ( n ) + 1 ) for j in range ( 0, n ): a[j] = float ( fac * 1 + j + 1 ) return a def r8ci_indicator_test ( ): #*****************************************************************************80 # ## R8CI_INDICATOR_TEST tests R8CI_INDICATOR. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_INDICATOR_TEST' ) print ( ' R8CI_INDICATOR sets up a R8CI indicator matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print ( n, a, ' The R8CI indicator matrix:' ) return def r8ci_mtv ( n, a, x ): #*****************************************************************************80 # ## R8CI_MTV multiplies a vector by a R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI matrix. # # Input, real X(N), the vector to be multiplied by A. # # Output, real B(N), the product A' * X. # import numpy as np b = np.zeros ( n ) for i in range ( 1, n + 1 ): for j in range ( 1, i + 1 ): b[i-1] = b[i-1] + a[i+1-j-1] * x[j-1] for j in range ( i + 1, n + 1 ): b[i-1] = b[i-1] + a[n+i+1-j-1] * x[j-1] return b def r8ci_mtv_test ( ): #*****************************************************************************80 # ## R8CI_MTV_TEST tests R8CI_MTV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_MTV_TEST' ) print ( ' R8CI_MTV computes b=A''*x where A is an R8CI indicator matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print ( n, a, ' The R8CI matrix A:' ) x = r8vec_indicator1 ( n ) r8vec_print ( n, x, ' The vector x:' ) b = r8ci_mtv ( n, a, x ) r8vec_print ( n, b, ' The product b=A''*x:' ) return def r8ci_mv ( n, a, x ): #*****************************************************************************80 # ## R8CI_MV multiplies a R8CI matrix times a vector. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI matrix. # # Input, real X(N), the vector to be multiplied by A. # # Output, real B(N), the product A * x. # import numpy as np b = np.zeros ( n ) for i in range ( 1, n + 1 ): for j in range ( 1, i ): b[i-1] = b[i-1] + a[n+j+1-i-1] * x[j-1] for j in range ( i, n + 1 ): b[i-1] = b[i-1] + a[j+1-i-1] * x[j-1] return b def r8ci_mv_test ( ): #*****************************************************************************80 # ## R8CI_MV_TEST tests R8CI_MV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_MV_TEST' ) print ( ' R8CI_MV computes b=A*x where A is an R8CI indicator matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print ( n, a, ' The R8CI matrix A:' ) x = r8vec_indicator1 ( n ) r8vec_print ( n, x, ' The vector x:' ) b = r8ci_mv ( n, a, x ) r8vec_print ( n, b, ' The product b=A*x:' ) return def r8ci_print ( n, a, title ): #*****************************************************************************80 # ## R8CI_PRINT prints a R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Input, real A(N), the R8CI matrix. # # Input, string TITLE, a title to be printed. # r8ci_print_some ( n, a, 0, 0, n - 1, n - 1, title ) return def r8ci_print_test ( ): #*****************************************************************************80 # ## R8CI_PRINT_TEST tests R8CI_PRINT. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_PRINT_TEST' ) print ( ' R8CI_PRINT prints an R8CI matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print ( n, a, ' The R8CI matrix:' ) return def r8ci_print_some ( n, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## R8CI_PRINT_SOME prints some of a R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Input, real A(N), the R8CI 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. # print ( '' ) print ( title ) incx = 5 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 ) i2hi = min ( ihi, n - 1 ) for i in range ( i2lo, i2hi + 1 ): print ( '%7d :' % ( i ), end = '' ) for j in range ( j2lo, j2hi + 1 ): if ( i <= j ): aij = a[j-i] else: aij = a[n+j-i] print ( '%12g ' % ( aij ), end = '' ) print ( '' ) return def r8ci_print_some_test ( ): #*****************************************************************************80 # ## R8CI_PRINTS_SOME_TEST tests R8CI_PRINT_SOME. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 10 print ( '' ) print ( 'R8CI_PRINT_SOME_TEST' ) print ( ' R8CI_PRINT_SOME prints some of an R8CI matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print_some ( n, a, 2, 3, 6, 5, ' Rows 2-6, Cols 3-5:' ) return def r8ci_random ( n ): #*****************************************************************************80 # ## r8ci_random() randomizes a R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Output, real A(N), the R8CI matrix. # from numpy.random import default_rng import numpy as np rng = default_rng ( ) a = np.zeros ( n ) for i in range ( 0, n ): a[i] = rng.uniform ( size = 1 ) return a def r8ci_random_test ( ): #*****************************************************************************80 # ## R8CI_RANDOM_TEST tests R8CI_RANDOM. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_RANDOM_TEST' ) print ( ' R8CI_RANDOM computes a random R8CI matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_random ( n ) r8ci_print ( n, a, ' The random R8CI matrix:' ) return def r8ci_sl ( n, a, b, job ): #*****************************************************************************80 # ## R8CI_SL solves a R8CI system. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI matrix. # # Input, real B(N), the right hand side. # # Input, integer JOB, specifies the system to solve. # 0, solve A * x = b. # nonzero, solve A' * x = b. # # Output, real X(N), the solution of the linear system. # import numpy as np x = np.zeros ( n ) work = np.zeros ( 2 * n - 2 ) if ( job == 0 ): # # Solve the system with the principal minor of order 1. # r1 = a[0] x[0] = b[0] / r1 r2 = 0.0 # # Recurrent process for solving the system. # for nsub in range ( 2, n + 1 ): # # Compute multiples of the first and last columns of # the inverse of the principal minor of order N. # r5 = a[n+2-nsub-1] r6 = a[nsub-1] if ( 2 < nsub ): work[nsub-2] = r2 for i in range ( 1, nsub - 1 ): r5 = r5 + a[n-i] * work[nsub-i-1] r6 = r6 + a[i] * work[n-2+i] r2 = - r5 / r1 r3 = - r6 / r1 r1 = r1 + r5 * r3 if ( 2 < nsub ): r6 = work[n-1] work[n+nsub-3] = 0.0 for i in range ( 2, nsub ): r5 = work[n-2+i] work[n-2+i] = work[i-1] * r3 + r6 work[i-1] = work[i-1] + r6 * r2 r6 = r5 work[n-1] = r3 # # Compute the solution of the system with the principal minor of order NSUB. # r5 = 0.0 for i in range ( 1, nsub ): r5 = r5 + a[n-i] * x[nsub-i-1] r6 = ( b[nsub-1] - r5 ) / r1 for i in range ( 1, nsub ): x[i-1] = x[i-1] + work[n+i-2] * r6 x[nsub-1] = r6 else: # # Solve the system with the principal minor of order 1. # r1 = a[0] x[0] = b[0] / r1 r2 = 0.0 # # Recurrent process for solving the system. # for nsub in range ( 2, n + 1 ): # # Compute multiples of the first and last columns of # the inverse of the principal minor of order N. # r5 = a[nsub-1] r6 = a[n+1-nsub] if ( 2 < nsub ): work[nsub-2] = r2 for i in range ( 1, nsub - 1 ): r5 = r5 + a[i] * work[nsub-i-1] r6 = r6 + a[n-i] * work[n-2+i] r2 = - r5 / r1 r3 = - r6 / r1 r1 = r1 + r5 * r3 if ( 2 < nsub ): r6 = work[n-1] work[n+nsub-3] = 0.0 for i in range ( 2, nsub ): r5 = work[n-2+i] work[n-2+i] = work[i-1] * r3 + r6 work[i-1] = work[i-1] + r6 * r2 r6 = r5 work[n-1] = r3 # # Compute the solution of the system with the principal minor of order NSUB. # r5 = 0.0 for i in range ( 1, nsub ): r5 = r5 + a[i] * x[nsub-i-1] r6 = ( b[nsub-1] - r5 ) / r1 for i in range ( 1, nsub ): x[i-1] = x[i-1] + work[n-2+i] * r6 x[nsub-1] = r6 return x def r8ci_sl_test ( ): #*****************************************************************************80 # ## R8CI_SL_TEST tests R8CI_SL. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 10 print ( '' ) print ( 'R8CI_SL_TEST' ) print ( ' R8CI_SL solves a circulant system.' ) print ( '' ) print ( ' Matrix order N = ', n ) # # Set the matrix. # a = r8ci_random ( n ) r8ci_print ( n, a, ' The circulant matrix:' ) for job in range ( 0, 2 ): # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # if ( job == 0 ): b = r8ci_mv ( n, a, x ) else: b = r8ci_mtv ( n, a, x ) # # Solve the linear system. # x = r8ci_sl ( n, a, b, job ) if ( job == 0 ): r8vec_print ( n, x, ' Solution:' ) else: r8vec_print ( n, x, ' Solution to transposed system:' ) return def r8ci_to_r8ge ( n, a ): #*****************************************************************************80 # ## R8CI_TO_R8GE copies a R8CI matrix into a R8GE matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # # Input, real A(N), the R8CI 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, i ): b[i,j] = a[n+j-i] for j in range ( i, n ): b[i,j] = a[j-i] return b def r8ci_to_r8ge_test ( ): #*****************************************************************************80 # ## R8CI_TO_R8GE_TEST tests R8CI_TO_R8GE. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_TO_R8GE_TEST' ) print ( ' R8CI_TO_R8GE converts an R8CI matrix to R8GE format.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_indicator ( n ) r8ci_print ( n, a, ' The R8CI matrix:' ) a_r8ge = r8ci_to_r8ge ( n, a ) print ( '' ) print ( ' The R8GE matrix:' ) print ( a_r8ge ) return def r8ci_zeros ( n ): #*****************************************************************************80 # ## R8CI_ZEROS zeros an R8CI matrix. # # Discussion: # # The R8CI storage format is used for an N by N circulant matrix. # An N by N circulant matrix A has the property that the entries on # row I appear again on row I+1, shifted one position to the right, # with the final entry of row I appearing as the first of row I+1. # The R8CI format simply records the first row of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the order of the matrix. # N must be positive. # # Output, real A(N), the R8CI matrix. # import numpy as np a = np.zeros ( n ) return a def r8ci_zeros_test ( ): #*****************************************************************************80 # ## R8CI_ZEROS_TEST tests R8CI_ZEROS. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 June 2016 # # Author: # # John Burkardt # n = 5 print ( '' ) print ( 'R8CI_ZEROS_TEST' ) print ( ' R8CI_ZEROS zeros an R8CI matrix.' ) print ( '' ) print ( ' Matrix order N = ', n ) a = r8ci_zeros ( n ) r8ci_print ( n, a, ' The zero R8CI matrix:' ) return 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 sort_safe_rc ( n, indx, isgn, i_save, j_save, k_save, l_save, n_save ): #*****************************************************************************80 # ## sort_safe_rc() externally sorts a list of items into ascending order. # # Discussion: # # This is a version of SORT_RC which does not rely on # storing certain work variables internally to the function. This makes # the function somewhat more awkward to call, but easier to program # in a variety of languages, and safe to use in a parallel programming # environment, or in cases where the sorting of several vectors is to # be carried out at more or less the same time. # # The actual list of data is not passed to the routine. Hence this # routine may be used to sort integers, reals, numbers, names, # dates, shoe sizes, and so on. After each call, the routine asks # the user to compare or interchange two items, until a special # return value signals that the sorting is completed. # # Example: # # n = 100 # indx = 0 # isgn = 0 # i_save = 0 # j_save = 0 # k_save = 0 # l_save = 0 # n_save = 0 # # while ( 1 ) # # indx, i, j, i_save, j_save, k_save, l_save, n_save = # sort_safe_rc ( n, indx, isgn, i_save, j_save, k_save, l_save, n_save ) # # if ( indx < 0 ) # # isgn = 1 # if ( a(i) <= a(j) ) # isgn = -1 # end # # elseif ( 0 < indx ) # # k = a(i) # a(i) = a(j) # a(j) = k # # else # # break # # end # # end # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 10 March 2015 # # Author: # # Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. # MATLAB version by John Burkardt # # Reference: # # Albert Nijenhuis, Herbert Wilf. # Combinatorial Algorithms, # Academic Press, 1978, second edition, # ISBN 0-12-519260-6. # # Parameters: # # Input, integer N, the number of items to be sorted. # # Input, integer INDX, the main communication signal. # The user must set INDX to 0 before the first call. # Thereafter, the user should set the input value of INDX # to the output value from the previous call. # # Input, integer ISGN, results of comparison of elements I and J. # (Used only when the previous call returned INDX less than 0). # ISGN <= 0 means I is less than or equal to J # 0 <= ISGN means I is greater than or equal to J. # # Output, integer INDX, the main communication signal. # If INDX is # * greater than 0, the user should: # interchange items I and J # call again. # * less than 0, the user should: # compare items I and J # set ISGN = -1 if I < J, ISGN = +1 if J < I # call again. # * equal to 0, the sorting is done. # # Output, integer I, J, the indices of two items. # On return with INDX positive, elements I and J should be interchanged. # On return with INDX negative, elements I and J should be compared, and # the result reported in ISGN on the next call. # # Input/output, integer I_SAVE, J_SAVE, K_SAVE, L_SAVE, N_SAVE, workspace # needed by the routine. Before calling the function, # the user should declare variables to hold these values, but should # not change them, and need not ever examine them. # # # INDX = 0: This is the first call. # if ( indx == 0 ): k_save = ( n // 2 ) l_save = k_save n_save = n # # INDX < 0: The user is returning the results of a comparison. # elif ( indx < 0 ): if ( indx == -2 ): if ( isgn < 0 ): i_save = i_save + 1 j_save = l_save l_save = i_save indx = -1 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save if ( 0 < isgn ): indx = 2 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save if ( k_save <= 1 ): if ( n_save == 1 ): i_save = 0 j_save = 0 indx = 0 else: i_save = n_save n_save = n_save - 1 j_save = 1 indx = 1 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save k_save = k_save - 1 l_save = k_save # # 0 < INDX, the user was asked to make an interchange. # elif ( indx == 1 ): l_save = k_save while ( True ): i_save = 2 * l_save if ( i_save == n_save ): j_save = l_save l_save = i_save indx = -1 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save elif ( i_save <= n_save ): j_save = i_save + 1 indx = -2 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save if ( k_save <= 1 ): break k_save = k_save - 1 l_save = k_save if ( n_save == 1 ): i_save = 0 j_save = 0 indx = 0 i = i_save j = j_save else: i_save = n_save n_save = n_save - 1 j_save = 1 indx = 1 i = i_save j = j_save return indx, i, j, i_save, j_save, k_save, l_save, n_save 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 ( ) r8ci_test ( ) timestamp ( )