#! /usr/bin/env python3
#
def r83_np_test ( ):
#*****************************************************************************80
#
## r83_np_test() tests r83_np().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 24 May 2016
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r83_np_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' Test r83_np().' )
r83_np_det_test ( )
r83_np_fa_test ( )
r83_np_fs_test ( )
r83_np_ml_test ( )
r83_np_sl_test ( )
#
# Terminate.
#
print ( '' )
print ( 'r83_np_test():' )
print ( ' Normal end of execution.' )
return
def r83_dif2 ( m, n ):
#*****************************************************************************80
#
## r83_dif2() returns the DIF2 matrix in R83 format.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# 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:
#
# 07 July 2015
#
# 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 M, N, the order of the matrix.
#
# Output, real A(3,N), the matrix.
#
import numpy as np
a = np.zeros( [ 3, n ] )
for j in range ( 0, n):
for i in range ( max ( 0, j - 1 ), min ( m, j + 2 ) ):
if ( i == j - 1 ):
a[i-j+1,j] = -1.0
elif ( i == j ):
a[i-j+1,j] = +2.0
elif ( i == j + 1 ):
a[i-j+1,j] = -1.0
return a
def r83_mtv ( m, n, a, x ):
#*****************************************************************************80
#
## r83_mtv() multiplies a vector by an R83 matrix.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3,N), the R83 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 j in range ( 0, n ):
for i in range ( max ( 0, j - 1 ), min ( m, j + 2 ) ):
b[j] = b[j] + x[i] * a[i-j+1,j]
return b
def r83_mv ( m, n, a, x ):
#*****************************************************************************80
#
## r83_mv() multiplies a R83 matrix times a vector.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the linear system.
#
# Input, real A(3,N), the R83 matrix.
#
# Input, real X(N), the vector to be multiplied by A.
#
# Output, real B(M), the product A * x.
#
import numpy as np
b = np.zeros ( m )
for j in range ( 0, n ):
for i in range ( max ( 0, j - 1 ), min ( m, j + 2 ) ):
b[i] = b[i] + a[i-j+1,j] * x[j]
return b
def r83_np_det ( n, a ):
#*****************************************************************************80
#
## r83_np_det() returns the determinant of a R83 system factored by r83_np_fa.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:N), the diagonal in
# entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
# original matrix is "collapsed" vertically into the array.
#
# Example:
#
# Here is how a R83 matrix of order 5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be at least 2.
#
# Input, real A(3,N), the tridiagonal factor information computed
# by r83_np_fa.
#
# Output, real DET, the determinant of the matrix.
#
det = 1.0
for j in range ( 0, n ):
det = det * a[1,j]
return det
def r83_np_det_test ( ):
#*****************************************************************************80
#
## r83_np_det_test() tests r83_np_det().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 May 2016
#
# Author:
#
# John Burkardt
#
n = 10
print ( '' )
print ( 'r83_np_det_test():' )
print ( ' r83_np_det() computes the determinant of a tridiagonal' )
print ( ' matrix factored by r83_np_fa.' )
print ( '' )
print ( ' Matrix order N = ', n )
#
# Set the matrix.
#
a = r83_dif2 ( n, n )
#
# Factor the matrix.
#
a_lu, info = r83_np_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r83_np_det_test(): Fatal error!' )
print ( ' r83_np_fa() returns INFO = ', info )
raise Exception ( 'r83_np_det_test(): Fatal error!' )
r83_print ( n, n, a_lu, ' The factored R83 matrix:' )
#
# Compute the determinant.
#
det = r83_np_det ( n, a_lu )
print ( '' )
print ( ' r83_np_det() computes determinant = ', det )
print ( ' Exact determinant is = ', n + 1 )
return
def r83_np_fa ( n, a ):
#*****************************************************************************80
#
## r83_np_fa() factors a R83 matrix without pivoting.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:N), the diagonal in
# entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
# original matrix is "collapsed" vertically into the array.
#
# Because this routine does not use pivoting, it can fail even when
# the matrix is not singular, and it is liable to make larger
# errors.
#
# r83_np_fa and R83_NP_SL may be preferable to the corresponding
# LINPACK routine SGTSL for tridiagonal systems, which factors and solves
# in one step, and does not save the factorization.
#
# Example:
#
# Here is how a R83 matrix of order 5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be at least 2.
#
# Input, real A(3,N), the tridiagonal matrix.
#
# Output, integer INFO, singularity flag.
# 0, no singularity detected.
# nonzero, the factorization failed on the INFO-th step.
#
# Output, real A_LU(3,N), factorization information.
#
info = 0
a_lu = a.copy ( )
for j in range ( 0, n - 1 ):
if ( a_lu[1,j] == 0.0 ):
info = i
print ( '' )
print ( 'r83_np_fa(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'r83_np_fa(): Fatal error!' )
#
# Store the multiplier in L.
#
a_lu[2,j] = a_lu[2,j] / a_lu[1,j]
#
# Modify the diagonal entry in the next column.
#
a_lu[1,j+1] = a_lu[1,j+1] - a_lu[2,j] * a_lu[0,j+1]
if ( a_lu[1,n-1] == 0.0 ):
info = n - 1
print ( '' )
print ( 'r83_np_fa(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'r83_np_fa(): Fatal error!' )
return a_lu, info
def r83_np_fa_test ( ):
#*****************************************************************************80
#
## r83_np_fa_test() tests r83_np_fa().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
n = 10
print ( '' )
print ( 'r83_np_fa_test():' )
print ( ' r83_np_fa() factors a tridiagonal matrix with no pivoting' )
print ( '' )
print ( ' Matrix order N = ', n )
#
# Set the matrix.
#
a = r83_random ( n, n )
r83_print ( n, n, a, ' The tridiagonal matrix:' )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r83_mv ( n, n, a, x )
x = np.zeros ( n )
#
# Factor the matrix.
#
a_lu, info = r83_np_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r83_np_fa_test(): Fatal error!' )
print ( ' The test matrix is singular.' )
raise Exception ( 'r83_np_fa_test(): Fatal error!' )
#
# Solve the linear system.
#
job = 0
x = r83_np_sl ( n, a_lu, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side, using the factored matrix.
#
job = 1
b = r83_np_ml ( n, a_lu, x, job )
#
# Solve the linear system.
#
job = 1
x = r83_np_sl ( n, a_lu, b, job )
r8vec_print ( n, x, ' Solution to tranposed system:' )
return
def r83_np_fs ( n, a, b ):
#*****************************************************************************80
#
## r83_np_fs() factors and solves a R83 system with no pivoting.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:N), the diagonal in
# entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
# original matrix is "collapsed" vertically into the array.
#
# This algorithm requires that each diagonal entry be nonzero.
# It does not use pivoting, and so can fail on systems that
# are actually nonsingular.
#
# Example:
#
# Here is how a R83 matrix of order 5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the linear system.
#
# Input/output, real A(3,N).
# On input, the tridiagonal matrix.
# On output, the data in these vectors has been overwritten
# by factorization information.
#
# Input, real B(N), the right hand side of the linear system.
#
# Output, real X(N), the solution of the linear system.
#
#
# The diagonal entries can't be zero.
#
for j in range ( 0, n ):
if ( a[1,j] == 0.0 ):
print ( '' )
print ( 'r83_np_fs(): Fatal error!' )
print ( ' A(2,%d) = 0.', j )
raise Exception ( 'r83_np_fs(): Fatal error!' )
x = b.copy ( )
for i in range ( 1, n ):
xmult = a[2,i-1] / a[1,i-1]
a[1,i] = a[1,i] - xmult * a[0,i]
x[i] = x[i] - xmult * x[i-1]
x[n-1] = x[n-1] / a[1,n-1]
for i in range ( n - 2, -1, -1 ):
x[i] = ( x[i] - a[0,i+1] * x[i+1] ) / a[1,i]
return x
def r83_np_fs_test ( ):
#*****************************************************************************80
#
## r83_np_fs_test() tests r83_np_fs().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
n = 10
print ( '' )
print ( 'r83_np_fs_test():' )
print ( ' r83_np_fs() factors and solves a tridiagonal linear system.' )
print ( '' )
print ( ' Matrix order N = ', n )
#
# Set the matrix elements.
#
a = r83_random ( n, n )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute b = A * x.
#
b = r83_mv ( n, n, a, x )
#
# Solve the system.
#
x = r83_np_fs ( n, a, b )
r8vec_print ( n, x, ' Solution:' )
return
def r83_np_ml ( n, a_lu, x, job ):
#*****************************************************************************80
#
## r83_np_ml() computes A * x or x * A, where A has been factored by r83_np_fa().
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:N), the diagonal in
# entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
# original matrix is "collapsed" vertically into the array.
#
# Example:
#
# Here is how a R83 matrix of order 5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 May 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be at least 2.
#
# Input, real A_LU(3,N), the LU factors from R83_FA.
#
# Input, real X(N), the vector to be multiplied by A.
#
# Input, integer JOB, specifies the product to find.
# 0, compute A * x.
# nonzero, compute A' * x.
#
# Output, real B(N), the product A*x or A'*x.
#
b = x.copy ( )
if ( job == 0 ):
#
# Compute X := U * X
#
for i in range ( 0, n):
b[i] = a_lu[1,i] * b[i]
if ( i < n - 1 ):
b[i] = b[i] + a_lu[0,i+1] * b[i+1]
#
# Compute X: = L * X.
#
for i in range ( n - 1, 0, -1 ):
b[i] = b[i] + a_lu[2,i-1] * b[i-1]
else:
#
# Compute X: = L' * X.
#
for i in range ( 0, n - 1 ):
b[i] = b[i] + a_lu[2,i] * b[i+1]
#
# Compute X: = U' * X.
#
for i in range ( n - 1, 0, -1 ):
b[i] = a_lu[1,i] * b[i]
b[i] = b[i] + a_lu[0,i] * b[i-1]
b[0] = a_lu[1,0] * b[0]
return b
def r83_np_ml_test ( ):
#*****************************************************************************80
#
## r83_np_ml_test() tests r83_np_ml().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
n = 10
print ( '' )
print ( 'r83_np_ml_test():' )
print ( ' r83_np_ml() computes A*x or A\'*x' )
print ( ' where A has been factored by r83_fa().' )
print ( '' )
print ( ' Matrix order N = ', n )
for job in range ( 0, 2 ):
#
# Set the matrix.
#
a = r83_random ( n, n )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
if ( job == 0 ):
b = r83_mv ( n, n, a, x )
else:
b = r83_mtv ( n, n, a, x )
#
# Factor the matrix.
#
a_lu, info = r83_np_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r83_np_ml_test(): Fatal error!' )
print ( ' r83_np_fa() declares the matrix is singular!' )
print ( ' The value of INFO is ', info )
raise Exception ( 'r83_np_ml_test(): Fatal error!' )
#
# Now multiply factored matrix times solution to get right hand side again.
#
b2 = r83_np_ml ( n, a_lu, x, job )
if ( job == 0 ):
r8vec2_print_some ( n, b, b2, 10, ' A*x and PLU*x:' )
else:
r8vec2_print_some ( n, b, b2, 10, ' A\'*x and (PLU)\'*x' )
return
def r83_np_sl ( n, a_lu, b, job ):
#*****************************************************************************80
#
## r83_np_sl() solves a R83 system factored by r83_np_fa().
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:N), the diagonal in
# entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
# original matrix is "collapsed" vertically into the array.
#
# Example:
#
# Here is how a R83 matrix of order 5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 May 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be at least 2.
#
# Input, real A_LU(3,N), the LU factor information
# returned by r83_np_fa.
#
# Input, real B(N), the right hand side of the linear system.
#
# 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.
#
x = b.copy ( )
if ( job == 0 ):
#
# Solve L * Y = B.
#
for i in range ( 1, n ):
x[i] = x[i] - a_lu[2,i-1] * x[i-1]
#
# Solve U * X = Y.
#
for i in range ( n - 1, -1, -1 ):
x[i] = x[i] / a_lu[1,i]
if ( 0 < i ):
x[i-1] = x[i-1] - a_lu[0,i] * x[i]
else:
#
# Solve U' * Y = B
#
for i in range ( 0, n ):
x[i] = x[i] / a_lu[1,i]
if ( i < n - 1 ):
x[i+1] = x[i+1] - a_lu[0,i+1] * x[i]
#
# Solve L' * X = Y.
#
for i in range ( n - 2, -1, -1 ):
x[i] = x[i] - a_lu[2,i] * x[i+1]
return x
def r83_np_sl_test ( ):
#*****************************************************************************80
#
## r83_np_sl_test() tests r83_np_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 May 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
n = 10
print ( '' )
print ( 'r83_np_sl_test():' )
print ( ' r83_np_sl() solves a linear system after the tridiagonal' )
print ( ' matrix has been factored by r83_np_fa().' )
print ( '' )
print ( ' Matrix order N = ', n )
#
# Set the matrix.
#
a = r83_random ( n, n )
r83_print ( n, n, a, ' The tridiagonal matrix:' )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r83_mv ( n, n, a, x )
x = np.zeros ( n )
#
# Factor the matrix.
#
[ a_lu, info ] = r83_np_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r83_np_sl_test(): Fatal error!' )
print ( ' The test matrix is singular.' )
raise Exception ( 'r83_np_sl_test(): Fatal error!' )
#
# Solve the linear system.
#
job = 0
x = r83_np_sl ( n, a_lu, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side, using the factored matrix.
#
job = 1
b = r83_np_ml ( n, a_lu, x, job )
#
# Solve the linear system.
#
job = 1
x = r83_np_sl ( n, a_lu, b, job )
r8vec_print ( n, x, ' Solution to tranposed system:' )
return
def r83_print ( m, n, a, title ):
#*****************************************************************************80
#
## r83_print() prints a R83 matrix.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3,N), the R83 matrix.
#
# Input, string TITLE, a title.
#
r83_print_some ( m, n, a, 0, 0, m - 1, n - 1, title )
return
def r83_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r83_print_some() prints some of a R83 matrix.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3,N), the R83 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
#
# Print the columns of the matrix, in strips of 5.
#
for j2lo in range ( jlo, jhi + 1, incx ):
j2hi = j2lo + incx - 1
j2hi = min ( j2hi, n )
j2hi = min ( j2hi, jhi )
inc = j2hi + 1 - j2lo
print ( '' )
print ( ' Col: ', end = '' )
for j in range ( j2lo, j2hi + 1 ):
print ( '%7d ' % ( j ), end = '' )
print ( '' )
print ( ' Row' )
print ( ' ---' )
#
# Determine the range of the rows in this strip.
#
i2lo = max ( ilo, 0 )
i2lo = max ( i2lo, j2lo - 1 )
i2hi = min ( ihi, m - 1 )
i2hi = min ( i2hi, j2hi + 1 )
for i in range ( i2lo, i2hi + 1 ):
#
# Print out (up to) 5 entries in row I, that lie in the current strip.
#
print ( '%5d:' % ( i ), end = '' )
for j2 in range ( 1, inc + 1 ):
j = j2lo - 1 + j2
if ( i - j + 1 < 0 or 2 < i - j + 1 ):
print ( ' ', end = '' )
else:
print ( '%14g' % ( a[i-j+1,j] ), end = '' )
print ( '' )
return
def r83_random ( m, n ):
#*****************************************************************************80
#
## r83_random() randomizes a R83 matrix.
#
# Discussion:
#
# The R83 storage format is used for a tridiagonal matrix.
# The superdiagonal is stored in entries (1,2:min(M+1,N)).
# The diagonal in entries (2,1:min(M,N)).
# The subdiagonal in (3,min(M-1,N)).
# R8GE A(I,J) = R83 A[I-J+1+J*3] (0 based indexing).
#
# Example:
#
# An R83 matrix of order 3x5 would be stored:
#
# * A12 A23 A34 *
# A11 A22 A33 * *
# A21 A32 * * *
#
# An R83 matrix of order 5x5 would be stored:
#
# * A12 A23 A34 A45
# A11 A22 A33 A44 A55
# A21 A32 A43 A54 *
#
# An R83 matrix of order 5x3 would be stored:
#
# * A12 A23
# A11 A22 A33
# A21 A32 A43
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the linear system.
#
# Output:
#
# real A(3,N), the R83 matrix.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
a = np.zeros ( [ 3, n ] )
for j in range ( 0, n ):
for i in range ( max ( 0, j - 1 ), min ( m, j + 2 ) ):
a[i-j+1,j] = rng.random ( size = 1 )
return a
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 r8vec2_print_some ( n, x1, x2, max_print, title ):
#*****************************************************************************80
#
## r8vec2_print_some() prints "some" of an R8VEC2.
#
# Discussion:
#
# An R8VEC2 is two R8VEC's.
#
# An R8VEC is a vector of R8 values.
#
# The user specifies MAX_print, the maximum number of lines to print.
#
# If N, the size of the vectors, is no more than MAX_print, then
# the entire vectors are printed, one entry of each per line.
#
# Otherwise, if possible, the first MAX_print-2 entries are printed,
# followed by a line of periods suggesting an omission,
# and the last entry.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the number of entries of the vectors.
#
# real X1(N), X2(N), the vector to be printed.
#
# integer MAX_print, the maximum number of lines to print.
#
# string TITLE, a title.
#
if ( max_print <= 0 ):
return
if ( n <= 0 ):
return
print ( '' )
print ( title )
print ( '' )
if ( n <= max_print ):
for i in range ( 0, n ):
print ( '%6d: %14g %14g' % ( i, x1[i], x2[i] ) )
elif ( 3 <= max_print ):
for i in range ( 0, max_print - 2 ):
print ( '%6d: %14g %14g' % ( i, x1[i], x2[i] ) )
print ( '...... .............. ..............' )
i = n - 1
print ( '%6d: %14g %14g' % ( i, x1[i], x2[i] ) )
else:
for i in range ( 0, max_print - 1 ):
print ( '%6d: %14g %14g' % ( i, x1[i], x2[i] ) )
i = max_print - 1
print ( '%6d: %14g %14g ...more entries...' % ( i, x1[i], x2[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 ( )
r83_np_test ( )
timestamp ( )