#! /usr/bin/env python3
#
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 i4_log_10_test ( ) :
#*****************************************************************************80
#
## i4_log_10_test() tests i4_log_10().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 09 May 2013
#
# Author:
#
# John Burkardt
#
import platform
n = 13
x = [ 0, 1, 2, 3, 9, 10, 11, 99, 101, -1, -2, -3, -9 ]
print ( '' )
print ( 'i4_log_10_test' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' i4_log_10: whole part of log base 10,' )
print ( '' )
print ( ' X, i4_log_10' )
print ( '' )
for i in range ( 0, n ):
j = i4_log_10 ( x[i] )
print ( '%6d %12d' % ( x[i], j ) )
return
def r8ge_print ( m, n, a, title ):
#*****************************************************************************80
#
## r8ge_print() prints an R8GE matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, the number of rows in A.
#
# integer N, the number of columns in A.
#
# real A(M,N), the matrix.
#
# string TITLE, a title.
#
r8ge_print_some ( m, n, a, 0, 0, m - 1, n - 1, title )
return
def r8ge_print_test ( ):
#*****************************************************************************80
#
## r8ge_print_test() tests r8ge_print().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 29 July 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'r8ge_print_test():' )
print ( ' r8ge_print() prints an R8GE matrix.' )
m = 4
n = 6
v = np.array ( [ \
[ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ],
[ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ],
[ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ],
[ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 )
r8ge_print ( m, n, v, ' Here is an R8GE:' )
return
def r8ge_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r8ge_print_some() prints out a portion of an R8GE.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 10 February 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of the matrix.
#
# real A(M,N), an M by N matrix to be printed.
#
# integer ILO, JLO, the first row and column to print.
#
# integer IHI, JHI, the last row and column to print.
#
# string TITLE, a title.
#
incx = 5
print ( '' )
print ( title )
if ( m <= 0 or n <= 0 ):
print ( '' )
print ( ' (None)' )
return
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, m - 1 )
for i in range ( i2lo, i2hi + 1 ):
print ( '%7d :' % ( i ), end = '' )
for j in range ( j2lo, j2hi + 1 ):
print ( '%12g ' % ( a[i,j] ), end = '' )
print ( '' )
return
def r8ge_print_some_test ( ):
#*****************************************************************************80
#
## r8ge_print_some_test() tests r8ge_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 29 July 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'r8ge_print_some_test' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_print_some prints some of an R8GE matrix.' )
m = 4
n = 6
v = np.array ( [ \
[ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ],
[ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ],
[ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ],
[ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 )
r8ge_print_some ( m, n, v, 0, 3, 2, 5, ' Rows 0:2, Cols 3:5:' )
return
def r8ncf_print ( m, n, nz_num, rowcol, a, title ):
#*****************************************************************************80
#
## r8ncf_print() prints a R8NCF matrix.
#
# Discussion:
#
# The R8NCF storage format stores NZ_NUM, the number of nonzeros,
# a real array containing the nonzero values, a 2 by NZ_NUM integer
# array storing the row and column of each nonzero entry.
#
# The R8NCF format is used by NSPCG. NSPCG requires that the information
# for the diagonal entries of the matrix must come first.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in the matrix.
#
# integer ROWCOL(2,NZ_NUM), the row and column indices
# of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# string TITLE, a title.
#
r8ncf_print_some ( m, n, nz_num, rowcol, a, 0, 0, m - 1, n - 1, title )
return
def r8ncf_print_some ( m, n, nz_num, rowcol, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r8ncf_print_some() prints some of a R8NCF matrix.
#
# Discussion:
#
# The R8NCF storage format stores NZ_NUM, the number of nonzeros,
# a real array containing the nonzero values, a 2 by NZ_NUM integer
# array storing the row and column of each nonzero entry.
#
# The R8NCF format is used by NSPCG. NSPCG requires that the information
# for the diagonal entries of the matrix must come first.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in the matrix.
#
# integer ROWCOL(2,NZ_NUM), the row and column indices
# of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# integer ILO, JLO, IHI, JHI, the first row and
# column, and the last row and column to be printed.
#
# string TITLE, a title.
#
import numpy as np
incx = 5
index = np.zeros ( incx, dtype = np.int32 )
print ( '' )
print ( title )
#
# 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 - 1 )
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 )
i2hi = min ( ihi, m - 1 )
for i in range ( i2lo, i2hi + 1 ):
#
# Print out (up to) 5 entries in row I, that lie in the current strip.
#
nonzero = False
for j2 in range ( 0, inc ):
index[j2] = -1
for k in range ( 0, nz_num ):
if ( i == rowcol[0,k] and j2lo <= rowcol[1,k] and rowcol[1,k] <= j2hi ):
j2 = rowcol[1,k] - j2lo + 1
if ( a[k] != 0.0 ):
index[j2-1] = k
nonzero = True
if ( nonzero ):
print ( '%5d ' % ( i ), end = '' )
for j2 in range ( 0, inc ):
if ( 0 <= index[j2] ):
aij = a[index[j2]]
else:
aij = 0.0
print ( '%14g' % ( aij ), end = '' )
print ( '' )
return
def r8st_cg ( n, nz_num, row, col, a, b, x_init ):
#*****************************************************************************80
#
## r8st_cg() uses the conjugate gradient method on an r8st system.
#
# Discussion:
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix.
#
# It is possible that a pair of indices (I,J) may occur more than
# once. Presumably, in this case, the intent is that the actual value
# of A(I,J) is the sum of all such entries. This is not a good thing
# to do, but I seem to have come across this in MATLAB.
#
# The r8st format is used by CSPARSE ("sparse triplet"), SLAP
# ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
#
# The matrix A must be a positive definite symmetric 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.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 09 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.
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# integer NZ_NUM, the number of nonzero elements in the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# real B(N), the right hand side vector.
#
# 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 = r8st_mv ( n, n, nz_num, row, col, 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 = r8st_mv ( n, n, nz_num, row, col, 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 r8st_cg_test ( rng ):
#*****************************************************************************80
#
## r8st_cg_test() tests r8st_cg().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 09 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# rng: the current random number generator.
#
import numpy as np
print ( '' )
print ( 'r8st_cg_test():' )
print ( ' r8st_cg() applies CG to an r8st matrix.' )
n = 10
nz_num = 3 * n - 2
#
# Set A to the [-1 2 -1] matrix.
#
row, col, a = r8st_dif2 ( n, n, nz_num )
#
# Choose a random solution.
#
x1 = rng.random ( size = n )
#
# Compute the corresponding right hand side.
#
b = r8st_mv ( n, n, nz_num, row, col, a, x1 )
#
# Call the CG routine.
#
x2 = np.ones ( n )
x3 = r8st_cg ( n, nz_num, row, col, a, b, x2 )
#
# Compute the residual.
#
r = r8st_res ( n, n, nz_num, row, col, 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 = %d' % ( n ) )
print ( ' Norm of residual ||Ax-b|| = %g' % ( r_norm ) )
print ( ' Norm of error ||x1-x2|| = %g' % ( e_norm ) )
return
def r8st_check ( m, n, nz_num, row, col ):
#*****************************************************************************80
#
## r8st_check() checks that a r8st matrix data structure is properly sorted.
#
# Discussion:
#
# This routine assumes that the data structure has been sorted,
# so that the entries of ROW are ascending sorted, and that the
# entries of COL are ascending sorted, within the group of entries
# that have a common value of ROW.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix.
#
# The r8st format is used by CSPARSE ("sparse triplet"), SLAP
# ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of
# the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# Output:
#
# bool CHECK, is TRUE if the matrix is properly defined.
#
check = True
#
# Check 1 <= ROW(*) <= M.
#
for k in range ( 0, nz_num ):
if ( row[k] < 0 or m <= row[k] ):
check = False
return check
#
# Check 1 <= COL(*) <= N.
#
for k in range ( 0, nz_num ):
if ( col[k] < 0 or n <= col[k] ):
check = False
return check
#
# Check that ROW(K) <= ROW(K+1).
#
for k in range ( 0, nz_num - 1 ):
if ( row[k+1] < row[k] ):
check = False
return check
#
# Check that, if ROW(K) == ROW(K+1), that COL(K) < COL(K+1).
#
for k in range ( 0, nz_num - 1 ):
if ( row[k] == row[k+1] ):
if ( col[k+1] <= col[k] ):
check = False
return check
return check
def r8st_diagonal ( m, n, nz_num, row, col, a ):
#*****************************************************************************80
#
## r8st_diagonal() reorders an r8st matrix so diagonal entries are first.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# This routine reorders the entries of A so that the first N entries
# are exactly the diagonal entries of the matrix, in order.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements
# of the matrix.
#
# Output:
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the updated row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the updated nonzero elements
# of the matrix.
#
found = 0
for k in range ( 0, nz_num ):
while ( row[k] == col[k] ):
if ( row[k] == k ):
found = found + 1
break
i = row[k]
j = row[i]
row[i] = row[k]
row[k] = j
j = col[i]
col[i] = col[k]
col[k] = j
t = a[i]
a[i] = a[k]
a[k] = t
found = found + 1
if ( min ( m, n ) <= found ):
break
if ( min ( m, n ) <= found ):
break
if ( found < min ( m, n ) ):
print ( '' )
print ( 'r8st_diagonal - Warning!' )
print ( ' Number of diagonal entries expected was %d' % ( min ( m, n ) ) )
print ( ' Number found was %d' % ( found ) )
return row, col, a
def r8st_diagonal_test ( ):
#*****************************************************************************80
#
## r8st_diagonal_test() tests r8st_diagonal().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
nz_num = 20
col = np.array ( [ \
4, 5, 1, 1, 2, 3, 3, 4, 0, 5, \
3, 5, 4, 0, 5, 2, 0, 1, 0, 2 ], dtype = np.int32 )
row = np.array ( [ \
0, 2, 3, 5, 4, 1, 5, 2, 0, 1, \
3, 5, 4, 3, 3, 2, 5, 1, 2, 3 ], dtype = np.int32 )
print ( '' )
print ( 'r8st_diagonal_test' )
print ( ' r8st_diagonal rearranges an r8st matrix' )
print ( ' so that the diagonal is listed first.' )
m = 6
n = 6
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
a = r8st_indicator ( m, n, nz_num, row, col )
print ( '' )
print ( ' Before rearrangement:' )
print ( ' K ROW(K) COL(K) A(K)' )
print ( '' )
for k in range ( 0, nz_num ):
print ( ' %8d %8d %8d %14.6g' % ( k, row[k], col[k], a[k] ) )
row, col, a = r8st_diagonal ( m, n, nz_num, row, col, a )
print ( '' )
print ( ' After rearrangement:' )
print ( ' K ROW(K) COL(K) A(K)' )
print ( '' )
for k in range ( 0, nz_num ):
print ( ' %8d %8d %8d %14.6g' % ( k, row[k], col[k], a[k] ) )
return
def r8st_dif2 ( m, n, nz_num ):
#*****************************************************************************80
#
## r8st_dif2() returns the DIF2 matrix in r8st format.
#
# Example:
#
# N = 5
#
# 2 -1 . . .
# -1 2 -1 . .
# . -1 2 -1 .
# . . -1 2 -1
# . . . -1 2
#
# Properties:
#
# A is banded, with bandwidth 3.
#
# A is tridiagonal.
#
# Because A is tridiagonal, it has property A (bipartite).
#
# A is a special case of the TRIS or tridiagonal scalar matrix.
#
# A is integral, therefore det ( A ) is integral, and
# det ( A ) * inverse ( A ) is integral.
#
# A is Toeplitz: constant along diagonals.
#
# A is symmetric: A' = A.
#
# Because A is symmetric, it is normal.
#
# Because A is normal, it is diagonalizable.
#
# A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
#
# A is positive definite.
#
# A is an M matrix.
#
# A is weakly diagonally dominant, but not strictly diagonally dominant.
#
# A has an LU factorization A = L * U, without pivoting.
#
# The matrix L is lower bidiagonal with subdiagonal elements:
#
# L(I+1,I) = -I/(I+1)
#
# The matrix U is upper bidiagonal, with diagonal elements
#
# U(I,I) = (I+1)/I
#
# and superdiagonal elements which are all -1.
#
# A has a Cholesky factorization A = L * L', with L lower bidiagonal.
#
# L(I,I) = sqrt ( (I+1) / I )
# L(I,I-1) = -sqrt ( (I-1) / I )
#
# The eigenvalues are
#
# LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
# = 4 SIN^2(I*PI/(2*N+2))
#
# The corresponding eigenvector X(I) has entries
#
# X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
#
# Simple linear systems:
#
# x = (1,1,1,...,1,1), A*x=(1,0,0,...,0,1)
#
# x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
#
# det ( A ) = N + 1.
#
# The value of the determinant can be seen by induction,
# and expanding the determinant across the first row:
#
# det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
# = 2 * N - (N-1)
# = N + 1
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 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.
#
# Input:
#
# integer M, N, the number of rows and columns.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# Output:
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
import numpy as np
row = np.zeros ( nz_num, dtype = np.int32 )
col = np.zeros ( nz_num, dtype = np.int32 )
a = np.zeros ( nz_num, dtype = np.float64 )
k = 0
for i in range ( 0, m ):
j = i - 1
if ( 0 <= j and j < n ):
row[k] = i
col[k] = j
a[k] = -1.0
k = k + 1
j = i
if ( 0 <= j and j < n ):
row[k] = i
col[k] = j
a[k] = 2.0
k = k + 1
j = i + 1
if ( 0 <= j and j < n ):
row[k] = i
col[k] = j
a[k] = -1.0
k = k + 1
return row, col, a
def r8st_dif2_test ( ):
#*****************************************************************************80
#
## r8st_dif2_test() tests r8st_dif2().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 16 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_dif2_test()' )
print ( ' r8st_dif2() sets an r8st matrix to the second difference.' )
m = 5
n = 5
nz_num = 3 * n - 2
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
r8st_print ( m, n, nz_num, row, col, a, ' r8st matrix A:' )
return
def r8st_ij_to_k ( nz_num, row, col, i, j ):
#*****************************************************************************80
#
## r8st_ij_to_k() seeks the compressed index of the (I,J) entry of A.
#
# Discussion:
#
# If A(I,J) is nonzero, then its value is stored in location K.
#
# This routine searches the r8st storage structure for the index K
# corresponding to (I,J), returning -1 if no such entry was found.
#
# This routine assumes that the data structure has been sorted,
# so that the entries of ROW are ascending sorted, and that the
# entries of COL are ascending sorted, within the group of entries
# that have a common value of ROW.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix.
#
# The r8st format is used by CSPARSE ("sparse triplet"), SLAP
# ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# integer I, J, the row and column indices of the
# matrix entry.
#
# Output:
#
# integer K, the r8st index of the (I,J) entry.
#
lo = 0
hi = nz_num - 1
while ( True ):
if ( hi < lo ):
k = -1
break
md = ( ( lo + hi ) // 2 )
if ( row[md] < i or ( row[md] == i and col[md] < j ) ):
lo = md + 1
elif ( i < row[md] or ( row[md] == i and j < col[md] ) ):
hi = md - 1
else:
k = md
break
return k
def r8st_ij_to_k_test ( ):
#*****************************************************************************80
#
## r8st_ij_to_k_test() tests r8st_ij_to_k().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
m = 7
n = 5
nz_num = 10
row = np.array ( [ 0, 0, 1, 1, 3, 3, 3, 4, 5, 6 ] )
col = np.array ( [ 1, 4, 0, 4, 0, 1, 2, 3, 3, 0 ] )
print ( '' )
print ( 'r8st_ij_to_k_test()' )
print ( ' r8st_ij_to_k() returns the r8st index of (I,J).' )
print ( '' )
print ( ' Matrix rows M = %d' % ( m ) )
print ( ' Matrix columns N = %d' % ( n ) )
print ( ' Matrix nonzeros = %d' % ( nz_num ) )
check = r8st_check ( m, n, nz_num, row, col )
if ( not check ):
print ( '' )
print ( 'r8st_check(): Warning!' )
print ( ' The matrix is not in the proper sorted format.' )
return
print ( '' )
print ( ' I J K' )
print ( '' )
for i in range ( 0, m ):
for j in range ( 0, n ):
k = r8st_ij_to_k ( nz_num, row, col, i, j )
print ( ' %8d %8d %8d' % ( i, j, k ) )
return
def r8st_indicator ( m, n, nz_num, row, col ):
#*****************************************************************************80
#
## r8st_indicator() sets up a r8st indicator matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero entries.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
# of the nonzero elements.
#
# Output:
#
# real A(NZ_NUM), the indicator matrix.
#
import numpy as np
fac = 10 ** ( i4_log_10 ( n ) + 1 )
a = np.zeros ( nz_num, dtype = np.float64 )
for k in range ( 0, nz_num ):
i = row[k]
j = col[k]
a[k] = float ( fac * ( i + 1 ) + ( j + 1 ) )
return a
def r8st_indicator_test ( ):
#*****************************************************************************80
#
## r8st_indicator_test() tests r8st_indicator().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
m = 6
n = 6
nz_num = 20
row = np.array ( [0,2,3,5,4,1,5,2,0,1,3,5,4,3,3,2,5,1,2,3], dtype = np.int32 )
col = np.array ( [4,5,1,1,2,3,3,4,0,5,3,5,4,0,5,2,0,1,0,2], dtype = np.int32 )
print ( '' )
print ( 'r8st_indicator_test' )
print ( ' r8st_indicator sets an r8st indicator matrix.' )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
a = r8st_indicator ( m, n, nz_num, row, col )
r8st_print ( m, n, nz_num, row, col, a, ' Matrix A:' )
return
def r8st_jac_sl ( n, nz_num, row, col, a, b, x, it_max ):
#*****************************************************************************80
#
## r8st_jac_sl() solves an r8st system using Jacobi iteration.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# This routine REQUIRES that the matrix be square, that the matrix
# have nonzero diagonal entries, and that the first N entries of
# the array A be exactly the diagonal entries of the matrix, in order.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column
# indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# real B(N), the right hand side of the linear system.
#
# real X(N), an approximate solution
# to the system.
#
# integer IT_MAX, the maximum number of iterations.
#
# Output:
#
# real X(N), an improved approximate solution
# to the system.
#
import numpy as np
x_new = np.zeros ( n )
#
# Ensure that the matrix lists diagonal entries first.
#
row, col, a = r8st_diagonal ( n, n, nz_num, row, col, a )
for it_num in range ( 0, it_max ):
#
# Initialize to right hand side.
#
for i in range ( 0, n ):
x_new[i] = b[i]
#
# Subtract off-diagonal terms.
#
for k in range ( n, nz_num ):
i = row[k]
j = col[k]
x_new[i] = x_new[i] - a[k] * x[j]
#
# Update.
#
for i in range ( 0, n ):
x[i] = x_new[i] / a[i]
return x
def r8st_jac_sl_test ( ):
#*****************************************************************************80
#
## r8st_jac_sl_test() tests r8st_jac_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
print ( '' )
print ( 'r8st_jac_sl_test()' )
print ( ' r8st_jac_sl() uses Jacobi iteration to solve a linear system' )
print ( ' with an r8st matrix.' )
m = 10
n = 10
nz_num = 3 * n - 2
it_max = 25
print ( '' )
print ( ' Matrix order M = %8d' % ( m ) )
print ( ' Matrix order N = %8d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %8d' % ( nz_num ) )
print ( ' Iterations per = %8d' % ( it_max ) )
#
# Set the matrix values.
#
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Solving A * x = b.' )
print ( '' )
title = 'Matrix A:'
r8st_print ( m, n, nz_num, row, col, a, title )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
print ( 'Solution x:' )
print ( x )
#
# Compute the corresponding right hand side.
#
b = r8st_mv ( n, n, nz_num, row, col, a, x )
print ( 'Right hand size b:' )
print ( b )
#
# Set the starting solution.
#
x = np.zeros ( n )
#
# Solve the linear system.
#
for i in range ( 0, 3 ):
x = r8st_jac_sl ( n, nz_num, row, col, a, b, x, it_max )
r8vec_print ( n, x, ' Current solution estimate:' )
return
def r8st_mtv ( m, n, nz_num, row, col, a, x ):
#*****************************************************************************80
#
## r8st_mtv() multiplies an R8VEC times an r8st matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column
# indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the 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 k in range ( 0, nz_num ):
i = col[k]
j = row[k]
b[i] = b[i] + a[k] * x[j]
return b
def r8st_mtv_test ( ):
#*****************************************************************************80
#
## r8st_mtv_test() tests r8st_mtv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_mtv_test' )
print ( ' r8st_mtv computes b=A\'*x, where A is an r8st matrix.' )
m = 5
n = 4
if ( m == n ):
nz_num = 3 * n - 2
else:
nz_num = 3 * n - 1
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
x = r8vec_indicator1 ( m )
r8vec_print ( m, x, ' x:' )
b = r8st_mtv ( m, n, nz_num, row, col, a, x )
r8vec_print ( n, b, ' b=A\'*x:' )
return
def r8st_mv ( m, n, nz_num, row, col, a, x ):
#*****************************************************************************80
#
## r8st_mv() multiplies an r8st matrix by an R8VEC.
#
# Discussion:
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix.
#
# It is possible that a pair of indices (I,J) may occur more than
# once. Presumably, in this case, the intent is that the actual value
# of A(I,J) is the sum of all such entries. This is not a good thing
# to do, but I seem to have come across this in MATLAB.
#
# The r8st format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
# ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 09 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of
# the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# real X(N), the vector to be multiplied by A.
#
# Output:
#
# real B(M), the product vector A*X.
#
import numpy as np
b = np.zeros ( m, dtype = np.float64 )
for k in range ( 0, nz_num ):
i = row[k]
j = col[k]
b[i] = b[i] + a[k] * x[j]
return b
def r8st_mv_test ( ):
#*****************************************************************************80
#
## r8st_mv_test() tests r8st_mv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_mv_test' )
print ( ' r8st_mv computes b=A*x, where A is an r8st matrix.' )
m = 5
n = 4
if ( m == n ):
nz_num = 3 * n - 2
else:
nz_num = 3 * n - 1
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
x = r8vec_indicator1 ( n )
r8vec_print ( n, x, ' x:' )
b = r8st_mv ( m, n, nz_num, row, col, a, x )
r8vec_print ( m, b, ' b=A*x:' )
return
def r8st_print ( m, n, nz_num, row, col, a, title ):
#*****************************************************************************80
#
## r8st_print() prints a r8st matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
# of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# string TITLE, a title.
#
r8st_print_some ( m, n, nz_num, row, col, a, 0, 0, m - 1, n - 1, title )
return
def r8st_print_test ( ):
#*****************************************************************************80
#
## r8st_print_test() tests r8st_print().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_print_test' )
print ( ' r8st_print prints an r8st matrix.' )
m = 5
n = 5
nz_num = 3 * n - 2
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
r8st_print ( m, n, nz_num, row, col, a, ' r8st matrix A:' )
return
def r8st_print_some ( m, n, nz_num, row, col, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r8st_print_some() prints some of an r8st matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in the matrix.
#
# integer ROW[NZ_NUM], COL[NZ_NUM], the row and column indices
# of the nonzero elements.
#
# real A[NZ_NUM], the nonzero elements
# of the matrix.
#
# integer ILO, JLO, IHI, JHI, the first row and
# column, and the last row and column to be printed.
#
# string TITLE, a title.
#
import numpy as np
incx = 5
index = np.zeros ( incx, dtype = np.int32 )
print ( '' )
print ( title )
#
# 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 - 1 )
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 )
i2hi = min ( ihi, m - 1 )
for i in range ( i2lo, i2hi + 1 ):
#
# Print out (up to) 5 entries in row I, that lie in the current strip.
#
nonzero = False
for j2 in range ( 0, inc ):
index[j2] = -1
for k in range ( 0, nz_num ):
if ( i == row[k] and j2lo <= col[k] and col[k] <= j2hi ):
j2 = col[k] - j2lo + 1
if ( a[k] != 0.0 ):
index[j2-1] = k
nonzero = True
if ( nonzero ):
print ( '%5d ' % ( i ), end = '' )
for j2 in range ( 0, inc ):
if ( 0 <= index[j2] ):
aij = a[index[j2]]
else:
aij = 0.0
print ( '%14g' % ( aij ), end = '' )
print ( '' )
return
def r8st_print_some_test ( ):
#*****************************************************************************80
#
## r8st_print_some_test() tests r8st_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_print_some_test' )
print ( ' r8st_print_some prints some of an r8st matrix.' )
m = 5
n = 5
nz_num = 3 * n - 2
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
r8st_print_some ( m, n, nz_num, row, col, a, 2, 3, 4, 5, ' Rows 2:4, Cols 3:5:' )
return
def r8st_random ( m, n, nz_num, row, col, rng ):
#*****************************************************************************80
#
## r8st_random() randomizes an r8st matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero entries.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
# of the nonzero elements.
#
# rng: the current random number generator.
#
# Output:
#
# real A(NZ_NUM), the indicator matrix.
#
import numpy as np
a = np.zeros ( nz_num, dtype = np.float64 )
for k in range ( 0, nz_num ):
i = row[k]
j = col[k]
a[k] = rng.random ( )
return a
def r8st_random_test ( rng ):
#*****************************************************************************80
#
## r8st_random_test() tests r8st_random().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# rng: the current random number generator.
#
import numpy as np
m = 6
n = 6
nz_num = 20
row = np.array ( [0,2,3,5,4,1,5,2,0,1,3,5,4,3,3,2,5,1,2,3], dtype = np.int32 )
col = np.array ( [4,5,1,1,2,3,3,4,0,5,3,5,4,0,5,2,0,1,0,2], dtype = np.int32 )
print ( '' )
print ( 'r8st_random_test()' )
print ( ' r8st_random() randomizes an r8st indicator matrix.' )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
a = r8st_random ( m, n, nz_num, row, col, rng )
r8st_print ( m, n, nz_num, row, col, a, ' Matrix A:' )
return
def r8st_read ( input_file, m, n, nz_num ):
#*****************************************************************************80
#
## r8st_read() reads an r8st matrix from a file.
#
# Discussion:
#
# This routine needs the value of NZ_NUM, which can be determined
# by a call to r8st_read_size.
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# string INPUT_FILE, the name of the file to be read.
#
# Unused, integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# Output:
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements
# of the matrix.
#
import numpy as np
row = np.zeros ( nz_num, dtype = np.int32 )
col = np.zeros ( nz_num, dtype = np.int32 )
a = np.zeros ( nz_num, dtype = np.float64 )
input_unit = open ( input_file, 'r' )
k = 0
for line in input_unit:
words = line.split ( )
row[k] = int ( words[0] )
col[k] = int ( words[1] )
a[k] = float ( words[2] )
k = k + 1
input_unit.close ( )
return row, col, a
def r8st_read_size ( input_file ):
#*****************************************************************************80
#
## r8st_read_size() reads the size of an r8st matrix from a file.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# string INPUT_FILE, the name of the file to be read.
#
# Output:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
import numpy as np
input_unit = open ( input_file, 'r' )
m = -1
n = -1
nz_num = 0
for line in input_unit:
words = line.split ( )
i = int ( words[0] )
m = max ( m, i )
j = int ( words[1] )
n = max ( n, j )
nz_num = nz_num + 1
input_unit.close ( )
return m, n, nz_num
def r8st_read_test ( ):
#*****************************************************************************80
#
## r8st_read_test() tests r8st_read().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_read_test' )
print ( ' r8st_read_size reads the size of an r8st matrix.' )
print ( ' r8st_read reads the r8st matrix from a file.' )
input_file = 'r8st_matrix.txt'
m, n, nz_num = r8st_read_size ( input_file )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
row, col, a = r8st_read ( input_file, m, n, nz_num )
r8st_print_some ( m, n, nz_num, row, col, a, 1, 1, \
10, 10, ' Initial 10x10 block of recovered r8st matrix:' )
return
def r8st_res ( m, n, nz_num, row, col, a, x, b ):
#*****************************************************************************80
#
## r8st_res() computes the residual R = B-A*X for r8st matrices.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 09 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.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and
# column indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements of the matrix.
#
# real X(N), the vector to be multiplied by A.
#
# real B(M), the desired result A * x.
#
# Output:
#
# real R(M), the residual R = B - A * X.
#
r = r8st_mv ( m, n, nz_num, row, col, a, x )
for i in range ( 0, m ):
r[i] = b[i] - r[i]
return r
def r8st_res_test ( ):
#*****************************************************************************80
#
## r8st_res_test() tests r8st_res().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_res_test' )
print ( ' r8st_res computes r=b-A*x, where A is an r8st matrix.' )
m = 5
n = 4
if ( m == n ):
nz_num = 3 * n - 2
else:
nz_num = 3 * n - 1
row, col, a = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
x = r8vec_indicator1 ( n )
r8vec_print ( n, x, ' x:' )
b = r8st_mv ( m, n, nz_num, row, col, a, x )
r = r8st_res ( m, n, nz_num, row, col, a, x, b )
r8vec_print ( m, r, ' r=b-A*x:' )
return
def r8st_to_r8ge ( m, n, nz_num, row, col, a_r8st ):
#*****************************************************************************80
#
## r8st_to_r8ge() copies an r8st matrix to an R8GE matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional bool
# array can be thought of as a vector of M*N entries, starting with
# the M entries in the column 1, then the M entries in column 2
# and so on. Considered as a vector, the entry A(I,J) is then stored
# in vector location I+(J-1)*M.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column
# indices of the nonzero elements.
#
# real A_r8st(NZ_NUM), the nonzero elements of the matrix.
#
# Output:
#
# real A_r8ge(M,N), the R8GE matrix.
#
import numpy as np
a_r8ge = np.zeros ( [ m, n ] )
for k in range ( 0, nz_num ):
i = row[k]
j = col[k]
a_r8ge[i,j] = a_r8ge[i,j] + a_r8st[k]
return a_r8ge
def r8st_to_r8ge_test ( ):
#*****************************************************************************80
#
## r8st_to_r8ge_test() tests r8st_to_r8ge().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_to_r8ge_test' )
print ( ' r8st_to_r8ge converts an r8st matrix to R8GE format.' )
m = 5
n = 5
nz_num = 3 * n - 2
row, col, a_r8st = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
r8st_print ( m, n, nz_num, row, col, a_r8st, ' r8st matrix A:' )
a_r8ge = r8st_to_r8ge ( m, n, nz_num, row, col, a_r8st )
r8ge_print ( m, n, a_r8ge, ' R8GE matrix A:' )
return
def r8st_to_r8ncf ( m, n, nz_num, row, col, a_r8st ):
#*****************************************************************************80
#
## r8st_to_r8ncf() copies an r8st matrix to an R8NCF matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# The R8NCF storage format stores NZ_NUM, the number of nonzeros,
# a real array containing the nonzero values, a 2 by NZ_NUM integer
# array storing the row and column of each nonzero entry.
#
# The R8NCF format is used by NSPCG. NSPCG requires that the information
# for the diagonal entries of the matrix must come first.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column
# indices of the nonzero elements.
#
# real A_r8st(NZ_NUM), the nonzero elements of the matrix.
#
# Output:
#
# integer NZ_NUM, the number of nonzero elements in the DNCF matrix.
#
# integer ROWCOL(2,NZ_NUM), the row and column indices
# of the nonzero elements.
#
# real A_r8ncf(M,N), the R8NCF matrix.
#
import numpy as np
rowcol = np.zeros ( [ 2, nz_num ], dtype = np.int32 )
a_r8ncf = np.zeros ( nz_num )
for k in range ( 0, nz_num ):
rowcol[0,k] = row[k]
rowcol[1,k] = col[k]
a_r8ncf[k] = a_r8st[k]
return nz_num, rowcol, a_r8ncf
def r8st_to_r8ncf_test ( ):
#*****************************************************************************80
#
## r8st_to_r8ncf_test() tests r8st_to_r8ncf().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_to_r8ncf_test' )
print ( ' r8st_to_r8ncf converts an r8st matrix to R8NCF format.' )
m = 5
n = 5
nz_num = 3 * n - 2
row, col, a_r8st = r8st_dif2 ( m, n, nz_num )
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
r8st_print ( m, n, nz_num, row, col, a_r8st, ' r8st matrix A:' )
nz_num, rowcol, a_r8ncf = r8st_to_r8ncf ( m, n, nz_num, row, col, a_r8st )
r8ncf_print ( m, n, nz_num, rowcol, a_r8ncf, ' R8NCF matrix A:' )
return
def r8st_write ( m, n, nz_num, row, col, a, output_file ):
#*****************************************************************************80
#
## r8st_write() writes an r8st matrix to a file.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2006
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero elements in
# the matrix.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column
# indices of the nonzero elements.
#
# real A(NZ_NUM), the nonzero elements
# of the matrix.
#
# string OUTPUT_FILE, the name of the file to which
# the information is to be written.
#
output_unit = open ( output_file, 'w' )
for k in range ( 0, nz_num ):
s = '%d %d %g\n' % ( row[k], col[k], a[k] )
output_unit.write ( s )
output_unit.close ( )
return
def r8st_write_test ( ):
#*****************************************************************************80
#
## r8st_write_test() tests r8st_write().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
print ( '' )
print ( 'r8st_write_test' )
print ( ' r8st_write writes an r8st matrix to a file.' )
m = 100
n = 100
nz_num = 3 * n - 2
print ( '' )
print ( ' Matrix order M = %d' % ( m ) )
print ( ' Matrix order N = %d' % ( n ) )
print ( ' Matrix nonzeros NZ_NUM = %d' % ( nz_num ) )
#
# Set the matrix values.
#
row, col, a = r8st_dif2 ( m, n, nz_num )
#
# Print a bit of the matrix.
#
r8st_print_some ( m, n, nz_num, row, col, a, 1, 1, \
10, 10, ' Initial 10x10 block of r8st matrix:' )
#
# Write the matrix to a file.
#
output_file = 'r8st_matrix.txt'
r8st_write ( m, n, nz_num, row, col, a, output_file )
print ( '' )
print ( ' Matrix data written to "%s".' % ( output_file ) )
return
def r8st_zeros ( m, n, nz_num, row, col ):
#*****************************************************************************80
#
## r8st_zeros() zeros an r8st matrix.
#
# Discussion:
#
# The r8st storage format corresponds to the SLAP Triad format.
#
# The r8st storage format stores the row, column and value of each nonzero
# entry of a sparse matrix. The entries may be given in any order. No
# check is made for the erroneous case in which a given matrix entry is
# specified more than once.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# integer NZ_NUM, the number of nonzero entries.
#
# integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
# of the nonzero elements.
#
# Output:
#
# real A(NZ_NUM), the indicator matrix.
#
import numpy as np
a = np.zeros ( nz_num )
return a
def r8st_zeros_test ( ):
#*****************************************************************************80
#
## r8st_zeros_test() tests r8st_zeros().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 September 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
m = 5
n = 4
nz_num = 11
row = np.array ( [0,1,2,3,0,1,1,2,3,4,4], dtype = np.int32 )
col = np.array ( [0,1,2,3,3,0,2,3,1,0,1], dtype = np.int32 )
print ( '' )
print ( 'r8st_zeros_test' )
print ( ' r8st_zeros zeros out space for an r8st matrix.' )
print ( '' )
print ( ' Matrix order M, N = %d, %d' % ( m, n ) )
a = r8st_zeros ( m, n, nz_num, row, col )
r8st_print ( m, n, nz_num, row, col, a, ' Matrix A:' )
return
def r8st_test ( ):
#*****************************************************************************80
#
## r8st_test() tests r8st().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 January 2024
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
import platform
print ( '' )
print ( 'r8st_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' Test r8st().' )
rng = default_rng ( )
r8st_cg_test ( rng )
r8st_diagonal_test ( )
r8st_dif2_test ( )
r8st_ij_to_k_test ( )
r8st_indicator_test ( )
r8st_jac_sl_test ( )
r8st_mtv_test ( )
r8st_mv_test ( )
r8st_print_test ( )
r8st_print_some_test ( )
r8st_random_test ( rng )
r8st_read_test ( )
r8st_res_test ( )
r8st_to_r8ge_test ( )
r8st_to_r8ncf_test ( )
r8st_write_test ( )
r8st_zeros_test ( )
#
# Terminate.
#
print ( '' )
print ( 'r8st_test():' )
print ( ' Normal end of execution.' )
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
a = numpy.zeros ( n );
for i in range ( 0, n ):
a[i] = i + 1
return a
def r8vec_indicator1_test ( ):
#*****************************************************************************80
#
## r8vec_indicator1_test() tests r8vec_indicator1().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 February 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'r8vec_indicator1_test' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8vec_indicator1 returns the 1-based indicator matrix.' )
n = 10
a = r8vec_indicator1 ( n )
r8vec_print ( n, a, ' The 1-based indicator vector:' )
return
def r8vec_print ( n, a, title ):
#*****************************************************************************80
#
## r8vec_print() prints an R8VEC.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 August 2014
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the dimension of the vector.
#
# real A(N), the vector to be printed.
#
# string TITLE, a title.
#
print ( '' )
print ( title )
print ( '' )
for i in range ( 0, n ):
print ( '%6d: %12g' % ( i, a[i] ) )
return
def timestamp ( ):
#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 August 2019
#
# Author:
#
# John Burkardt
#
import time
t = time.time ( )
print ( time.ctime ( t ) )
return
if ( __name__ == '__main__' ):
timestamp ( )
r8st_test ( )
timestamp ( )