#! /usr/bin/env python3
#
def r83s_test ( ):
#*****************************************************************************80
#
## r83s_test() tests r83s().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 17 August 2022
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r83s_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' Test r83s().' )
r83s_cg_test ( )
r83s_dif2_test ( )
r83s_gs_sl_test ( )
r83s_indicator_test ( )
r83s_jac_sl_test ( )
r83s_mtv_test ( )
r83s_mv_test ( )
r83s_print_test ( )
r83s_print_some_test ( )
r83s_random_test ( )
r83s_res_test ( )
r83s_to_r8ge_test ( )
r83s_zeros_test ( )
#
# Terminate.
#
print ( '' )
print ( 'r83s_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 r83s_cg ( n, a, b, x_init ):
#*****************************************************************************80
#
## r83s_cg() uses the conjugate gradient method on an R83S system.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# The matrix A must be a positive definite symmetric band matrix.
#
# The method is designed to reach the solution after N computational
# steps. However, roundoff may introduce unacceptably large errors for
# some problems. In such a case, calling the routine again, using
# the computed solution as the new starting estimate, should improve
# the results.
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 July 2015
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Frank Beckman,
# The Solution of Linear Equations by the Conjugate Gradient Method,
# in Mathematical Methods for Digital Computers,
# edited by John Ralston, Herbert Wilf,
# Wiley, 1967,
# ISBN: 0471706892,
# LC: QA76.5.R3.
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be positive.
#
# Input, real A(3), the matrix.
#
# Input, real B(N), the right hand side vector.
#
# Input, real X_INIT(N), an estimate for the solution, which may be 0.
#
# Output, real X(N), the approximate solution vector.
#
import numpy as np
x = np.zeros ( n )
for i in range ( 0, n ):
x[i] = x_init[i]
#
# Initialize
# AP = A * x,
# R = b - A * x,
# P = b - A * x.
#
ap = r83s_mv ( n, n, a, x )
r = np.zeros ( n )
for i in range ( 0, n ):
r[i] = b[i] - ap[i]
p = np.zeros ( n )
for i in range ( 0, n ):
p[i] = b[i] - ap[i]
#
# Do the N steps of the conjugate gradient method.
#
for it in range ( 0, n ):
#
# Compute the matrix*vector product AP=A*P.
#
ap = r83s_mv ( n, n, a, p )
#
# Compute the dot products
# PAP = P*AP,
# PR = P*R
# Set
# ALPHA = PR / PAP.
#
pap = np.dot ( p, ap )
pr = np.dot ( p, r )
if ( pap == 0.0 ):
return x
alpha = pr / pap
#
# Set
# X = X + ALPHA * P
# R = R - ALPHA * AP.
#
for i in range ( 0, n ):
x[i] = x[i] + alpha * p[i]
for i in range ( 0, n ):
r[i] = r[i] - alpha * ap[i]
#
# Compute the vector dot product
# RAP = R*AP
# Set
# BETA = - RAP / PAP.
#
rap = np.dot ( r, ap )
beta = - rap / pap
#
# Update the perturbation vector
# P = R + BETA * P.
#
for i in range ( 0, n ):
p[i] = r[i] + beta * p[i]
return x
def r83s_cg_test ( ):
#*****************************************************************************80
#
## r83s_cg_test() tests r83s_cg().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 July 2015
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
import numpy as np
print ( '' )
print ( 'r83s_cg_test():' )
print ( ' r83s_cg() applies CG to an R83S matrix.' )
rng = default_rng ( )
n = 10
#
# Let A be the -1 2 -1 matrix.
#
a = r83s_dif2 ( n, n )
#
# Choose a random solution.
#
x1 = rng.random ( size = n )
#
# Compute the corresponding right hand side.
#
b = r83s_mv ( n, n, a, x1 )
#
# Call the CG routine.
#
x2 = np.ones ( n )
x3 = r83s_cg ( n, a, b, x2 )
#
# Compute the residual.
#
r = r83s_res ( n, n, a, x3, b )
r_norm = np.linalg.norm ( r )
#
# Compute the error.
#
e_norm = np.linalg.norm ( x1 - x3 )
#
# Report.
#
print ( '' )
print ( ' Number of variables N = ', n )
print ( ' Norm of residual ||Ax-b|| = ', r_norm )
print ( ' Norm of error ||x1-x2|| = ', e_norm )
return
def r83s_dif2 ( m, n ):
#*****************************************************************************80
#
## r83s_dif2() returns the DIF2 matrix in R83S format.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# 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), the matrix.
#
import numpy as np
a = np.array ( [ -1.0, 2.0, -1.0 ] )
return a
def r83s_dif2_test ( ):
#*****************************************************************************80
#
## r83s_dif2_test() tests r83s_dif2().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_DIF2_TEST' )
print ( ' R83S_DIF2 sets an R83S matrix to the second difference.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a = r83s_dif2 ( m, n )
r83s_print ( m, n, a, ' Second difference in R83S format:' )
return
def r83s_gs_sl ( n, a, b, x, it_max ):
#*****************************************************************************80
#
## r83s_gs_sl() solves an R83S system using Gauss-Seidel iteration.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A(I-J+2).
#
# This routine simply applies a given number of steps of the
# iteration to an input approximate solution. On first call, you can
# simply pass in the zero vector as an approximate solution. If
# the returned value is not acceptable, you may call again, using
# it as the starting point for additional iterations.
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
#
# Input, real A(3), the R83S matrix.
#
# Input, real B(N), the right hand side of the linear system.
#
# Input, real X(N), an approximate solution to the system.
#
# Input, integer IT_MAX, the maximum number of iterations.
#
# Output, real X(N), an approximate solution to the system.
#
#
# No diagonal matrix entry can be zero.
#
if ( a[1] == 0.0 ):
print ( '' )
print ( 'R83S_GS_SL - Fatal error!' )
print ( ' Zero diagonal entry.' )
raise Exception ( 'R83S_GS_SL - Fatal error!' )
for it_num in range ( 0, it_max ):
x[0] = ( b[0] - a[2] * x[1] ) / a[1]
for i in range ( 1, n - 1 ):
x[i] = ( b[i] - a[0] * x[i-1] - a[2] * x[i+1] ) / a[1]
x[n-1] = ( b[n-1] - a[0] * x[n-2] ) / a[1]
return x
def r83s_gs_sl_test ( ):
#*****************************************************************************80
#
## r83s_gs_sl_test() tests r83s_gs_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 September 2015
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
import numpy as np
print ( '' )
print ( 'r83s_gs_sl_test():' )
print ( ' r83s_gs_sl() applies Gauss-Seidel iteration with an R83S matrix' )
print ( ' to solve a linear system A*x=b.' )
rng = default_rng ( )
#
# Set A to the second difference matrix.
#
n = 10
a = r83s_dif2 ( n, n )
#
# Choose a random solution.
#
x1 = rng.random ( size = n )
#
# Compute the corresponding right hand side.
#
b = r83s_mv ( n, n, a, x1 )
#
# Call the Gauss-Seidel routine.
#
x2 = np.ones ( n )
it_max = 100
x3 = r83s_gs_sl ( n, a, b, x2, it_max )
#
# Compute the residual.
#
r = r83s_res ( n, n, a, x3, b )
r_norm = np.linalg.norm ( r )
#
# Compute the error.
#
e_norm = np.linalg.norm ( x1 - x3 )
#
# Report.
#
print ( '' )
print ( ' Number of variables N = ', n )
print ( ' Norm of residual ||Ax-b|| = ', r_norm )
print ( ' Norm of error ||x1-x2|| = ', e_norm )
return
def r83s_indicator ( m, n ):
#*****************************************************************************80
#
## r83s_indicator() sets an R83S indicator matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A(I-J+2).
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Output, real A(3), the R83S matrix.
#
import numpy as np
a = np.array ( [ 3.0, 2.0, 1.0 ] )
return a
def r83s_indicator_test ( ):
#*****************************************************************************80
#
## r83s_indicator_test() tests r83s_indicator().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_INDICATOR_TEST' )
print ( ' R83S_INDICATOR sets an R83S indicator matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a = r83s_indicator ( m, n )
r83s_print ( m, n, a, ' R83S indicator matrix:' )
return
def r83s_jac_sl ( n, a, b, x, it_max ):
#*****************************************************************************80
#
## r83s_jac_sl() solves an R83S system using Jacobi iteration.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A(I-J+2).
#
# This routine simply applies a given number of steps of the
# iteration to an input approximate solution. On first call, you can
# simply pass in the zero vector as an approximate solution. If
# the returned value is not acceptable, you may call again, using
# it as the starting point for additional iterations.
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
#
# Input, real A(3), the R83S matrix.
#
# Input, real B(N), the right hand side of the linear system.
#
# Input, real X(N), an approximate solution to the system.
#
# Input, integer IT_MAX, the maximum number of iterations.
#
# Output, real X(N), an approximate solution
# to the system.
#
import numpy as np
#
# No diagonal matrix entry can be zero.
#
if ( a[1] == 0.0 ):
print ( '' )
print ( 'R83S_JAC_SL - Fatal error!' )
print ( ' Zero diagonal entry.' )
raise Exception ( 'R83S_JAC_SL - Fatal error!' )
x_new = np.zeros ( n )
for it_num in range ( 0, it_max ):
x_new[0] = b[0] - a[2] * x[1]
for i in range ( 1, n - 1 ):
x_new[i] = b[i] - a[0] * x[i-1] - a[2] * x[i+1]
x_new[n-1] = b[n-1] - a[0] * x[n-2]
#
# Divide by diagonal terms.
#
for i in range ( 0, n ):
x[i] = x_new[i] / a[1]
return x
def r83s_jac_sl_test ( ):
#*****************************************************************************80
#
## r83s_jac_sl_test() tests r83s_jac_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
import numpy as np
print ( '' )
print ( 'R83S_JAC_SL_TEST' )
print ( ' R83S_JAC_SL applies Jacobi iteration with an R83S matrix' )
print ( ' to solve a linear system A*x=b.' )
rng = default_rng ( )
#
# Set A to the second difference matrix.
#
n = 10
a = r83s_dif2 ( n, n )
#
# Choose a random solution.
#
x1 = rng.random ( size = n )
#
# Compute the corresponding right hand side.
#
b = r83s_mv ( n, n, a, x1 )
#
# Call the Jacobi routine.
#
x2 = np.ones ( n )
it_max = 100
x3 = r83s_jac_sl ( n, a, b, x2, it_max )
#
# Compute the residual.
#
r = r83s_res ( n, n, a, x3, b )
r_norm = np.linalg.norm ( r )
#
# Compute the error.
#
e_norm = np.linalg.norm ( x1 - x3 )
#
# Report.
#
print ( '' )
print ( ' Number of variables N = ', n )
print ( ' Norm of residual ||Ax-b|| = ', r_norm )
print ( ' Norm of error ||x1-x2|| = ', e_norm )
return
def r83s_mtv ( m, n, a, x ):
#*****************************************************************************80
#
## r83s_mtv() computes b=A'*x, where A is an R83S matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A(I-J+2).
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3), the R83S matrix.
#
# Input, real X(M), 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] + a[i-j+1] * x[i]
return b
def r83s_mtv_test ( ):
#*****************************************************************************80
#
## r83s_mtv_test() tests r83s_mtv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r83s_mtv_test():' )
print ( ' r83s_mtv() computes b=A\'*x, where A is an R83S matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a_83s = r83s_random ( m, n )
x = r8vec_indicator1 ( m )
ax_83s = r83s_mtv ( m, n, a_83s, x )
a_ge = r83s_to_r8ge ( m, n, a_83s )
ax_ge = r8ge_mtv ( m, n, a_ge, x )
r8vec2_print ( ax_83s, ax_ge, ' Product comparison:' )
return
def r83s_mv ( m, n, a, x ):
#*****************************************************************************80
#
## r83s_mv() multiplies an R83S matrix times an R8VEC.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 July 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the number of rows and columns.
#
# Input, real A(3), the R83S 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] * x[j]
return b
def r83s_mv_test ( ):
#*****************************************************************************80
#
## r83s_mv_test() tests r83s_mv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r83s_mv_test():' )
print ( ' r83s_mv() computes b=A*x, where A is an R83S matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a_83s = r83s_random ( m, n )
x = r8vec_indicator1 ( n )
ax_83s = r83s_mv ( m, n, a_83s, x )
a_ge = r83s_to_r8ge ( m, n, a_83s )
ax_ge = r8ge_mv ( m, n, a_ge, x )
r8vec2_print ( ax_83s, ax_ge, ' Product comparison:' )
return
def r83s_print ( m, n, a, title ):
#*****************************************************************************80
#
## r83s_print() prints an R83S matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3), the R83S matrix.
#
# Input, string TITLE, a title.
#
r83s_print_some ( m, n, a, 1, 1, m, n, title )
return
def r83s_print_test ( ):
#*****************************************************************************80
#
## r83s_print_test() tests r83s_print().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_PRINT_TEST' )
print ( ' R83S_PRINT prints an R83S matrix.' )
m = 5
n = 5
a = r83s_indicator ( m, n )
r83s_print ( m, n, a, ' R83S matrix:' )
return
def r83s_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r83s_print_some() prints some of a R83 matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3), the R83S 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 + 2 < 1 or 3 < i - j + 2 ):
print ( ' ', end = '' )
else:
print ( '%14g' % ( a[i-j+1] ), end = '' )
print ( '' )
return
def r83s_print_some_test ( ):
#*****************************************************************************80
#
## r83s_print_some_test() tests r83s_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_PRINT_SOME_TEST' )
print ( ' R83S_PRINT_SOME prints some of an R83S matrix.' )
m = 5
n = 5
a = r83s_indicator ( m, n )
r83s_print_some ( m, n, a, 2, 2, 5, 4, ' Rows 2-5, Cols 2-4:' )
return
def r83s_random ( m, n ):
#*****************************************************************************80
#
## r83s_random() randomizes an R83S matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# Output:
#
# real A(3), the R83 matrix.
#
from numpy.random import default_rng
rng = default_rng ( )
a = rng.random ( size = 3 )
return a
def r83s_random_test ( ):
#*****************************************************************************80
#
## r83s_random_test() tests r83s_random().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 Septemer 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_RANDOM_TEST' )
print ( ' R83S_RANDOM randomizes an R83S matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a = r83s_random ( m, n )
r83s_print ( m, n, a, ' Random R83S matrix:' )
return
def r83s_res ( m, n, a, x, b ):
#*****************************************************************************80
#
## r83s_res() computes the residual R = B-A*X for R83S matrices.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A(3), the matrix.
#
# Input, real X(N), the vector to be multiplied by A.
#
# Input, real B(M), the desired result A * x.
#
# Output, real R(M), the residual R = B - A * X.
#
r = r83s_mv ( m, n, a, x )
for i in range ( 0, m ):
r[i] = b[i] - r[i]
return r
def r83s_res_test ( ):
#*****************************************************************************80
#
## r83s_res_test() tests r83s_res().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r83s_res_test():' )
print ( ' r83s_res() computes b-A*x, where A is an R83S matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a = r83s_random ( m, n )
x = r8vec_indicator1 ( n )
b = r83s_mv ( m, n, a, x )
r = r83s_res ( m, n, a, x, b )
r8vec_print ( m, r, ' Residual A*x-b:' )
return
def r83s_to_r8ge ( m, n, a_83s ):
#*****************************************************************************80
#
## r83s_to_r8ge() copies an R83S matrix to a R8GE matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Input, real A_83S(3), the R83S matrix.
#
# Output, real A_GE(M,N), the R8GE matrix.
#
import numpy as np
a_ge = np.zeros ( [ m, n ] )
for j in range ( 0, n ):
for i in range ( max ( 0, j - 1 ), min ( m, j + 2 ) ):
a_ge[i,j] = a_83s[i-j+1]
return a_ge
def r83s_to_r8ge_test ( ):
#*****************************************************************************80
#
## r83s_to_r8ge_test() tests r83s_to_r8ge().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'R83S_TO_R8GE_TEST' )
print ( ' R83S_TO_R8GE converts an R83S matrix to R8GE format.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a_83s = r83s_random ( m, n )
r83s_print ( m, n, a_83s, ' R83s matrix:' )
a_ge = r83s_to_r8ge ( m, n, a_83s )
print ( '' )
print ( ' The R8GE matrix:' )
print ( a_ge )
return
def r83s_zeros ( m, n ):
#*****************************************************************************80
#
## r83s_zeros() zeros an R83S matrix.
#
# Discussion:
#
# The R83S storage format is used for a tridiagonal scalar matrix.
# The vector A(3) contains the subdiagonal, diagonal, and superdiagonal
# values that occur on every row.
# RGE A(I,J) = R83S A[I-J+1].
#
# Example:
#
# Here is how an R83S matrix of order 5, stored as (A1,A2,A3), would
# be interpreted:
#
# A2 A1 0 0 0
# A3 A2 A1 0 0
# 0 A3 A2 A1 0
# 0 0 A3 A2 A1
# 0 0 0 A3 A2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the order of the matrix.
#
# Output, real A(3,1), the R83 matrix.
#
import numpy as np
a = np.zeros ( 3 )
return a
def r83s_zeros_test ( ):
#*****************************************************************************80
#
## r83s_zeros_test() tests r83s_zeros().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r83s_zeros_test():' )
print ( ' r83s_zeros() zeros an R83S matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 1, 4 ):
if ( i == 1 ):
m = 3
n = 5
elif ( i == 2 ):
m = 5
n = 5
elif ( i == 3 ):
m = 5
n = 3
a = r83s_zeros ( m, n )
r83s_print ( m, n, a, ' Zeroed R83S matrix:' )
return
def r8ge_mtv ( m, n, a, x ):
#*****************************************************************************80
#
## r8ge_mtv() multiplies a vector by a R8GE matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A 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:
#
# 06 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, the number of rows of the matrix.
# M must be positive.
#
# integer N, the number of columns of the matrix.
# N must be positive.
#
# real A(M,N), the R8GE matrix.
#
# real X(M), 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 ( 0, m ):
b[j] = b[j] + x[i] * a[i,j]
return b
def r8ge_mv ( m, n, a, x ):
#*****************************************************************************80
#
## r8ge_mv() multiplies an R8GE matrix times a vector.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A 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:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, the number of rows of the matrix.
# M must be positive.
#
# integer N, the number of columns of the matrix.
# N must be positive.
#
# real A(M,N), the R8GE matrix.
#
# 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 i in range ( 0, m ):
for j in range ( 0, n ):
b[i] = b[i] + a[i,j] * x[j]
return b
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 ( a1, a2, title ):
#*****************************************************************************80
#
## r8vec2_print() prints an R8VEC2.
#
# Discussion:
#
# An R8VEC2 is a dataset consisting of N pairs of real values, stored
# as two separate vectors A1 and A2.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 June 2020
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the number of components of the vector.
#
# real A1(N), A2(N), the vectors to be printed.
#
# string TITLE, a title.
#
n = len ( a1 )
print ( '' )
print ( title )
print ( '' )
for i in range ( 0, n ):
print ( ' %6d: %12g %12g' % ( i, a1[i], a2[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 ( )
r83s_test ( )
timestamp ( )