#! /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.
#
i = int ( 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 numpy as np
n = 13
x = np.array ( [ 0, 1, 2, 3, 9, 10, 11, 99, 101, -1, -2, -3, -9 ] )
print ( '' )
print ( 'i4_log_10_test():' )
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 i4vec_print ( n, a, title ):
#*****************************************************************************80
#
## i4vec_print() prints an I4VEC.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 August 2014
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the dimension of the vector.
#
# integer A(N), the vector to be printed.
#
# string TITLE, a title.
#
print ( '' )
print ( title )
print ( '' )
for i in range ( 0, n ):
print ( '%6d %6d' % ( i, a[i] ) )
return
def i4vec_print_test ( ):
#*****************************************************************************80
#
## i4vec_print_test() tests i4vec_print().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 October 2014
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'i4vec_print_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' i4vec_print prints an I4VEC.' )
n = 4
v = np.array ( [ 91, 92, 93, 94 ], dtype = np.int32 )
i4vec_print ( n, v, ' Here is an I4VEC:' )
return
def r8_sign ( x ):
#*****************************************************************************80
#
## r8_sign() returns the sign of an R8.
#
# Discussion:
#
# The value is +1 if the number is positive or zero, and it is -1 otherwise.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 June 2013
#
# Author:
#
# John Burkardt
#
# Input:
#
# real X, the number whose sign is desired.
#
# Output:
#
# real VALUE, the sign of X.
#
if ( x < 0.0 ):
value = -1.0
else:
value = +1.0
return value
def r8_sign_test ( ):
#*****************************************************************************80
#
## r8_sign_test() tests r8_sign().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 September 2014
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
test_num = 5
r8_test = np.array ( [ -1.25, -0.25, 0.0, +0.5, +9.0 ] )
print ( '' )
print ( 'r8_sign_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8_sign returns the sign of an R8.' )
print ( '' )
print ( ' R8 r8_sign(R8)' )
print ( '' )
for test in range ( 0, test_num ):
r8 = r8_test[test]
s = r8_sign ( r8 )
print ( ' %8.4f %8.0f' % ( r8, s ) )
return
def r8ge_cg ( n, a, b, x ):
#*****************************************************************************80
#
## r8ge_cg() uses the conjugate gradient method on an R8GE system.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# 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.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 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.
#
# real A(N,N), the matrix.
#
# real B(N), the right hand side vector.
#
# real X(N), an estimate for the solution, which may be 0.
#
# Output:
#
# real X(N), the approximate solution.
#
import numpy as np
#
# Initialize
# AP = A * x,
# R = b - A * x,
# P = b - A * x.
#
ap = r8ge_mv ( n, n, a, x )
r = np.zeros ( n, dtype = np.float64 )
for i in range ( 0, n ):
r[i] = b[i] - ap[i]
p = np.zeros ( n, dtype = np.float64 )
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 = r8ge_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 r8ge_cg_test ( ):
#*****************************************************************************80
#
## r8ge_cg_test() tests r8ge_cg() for a full storage matrix.
#
# 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
rng = default_rng ( )
print ( '' )
print ( 'r8ge_cg_test():' )
print ( ' r8ge_cg applies CG to an R8GE matrix.' )
#
# Choose a random positive definite symmetric matrix A.
#
n = 10
key = 123456789
a = r8ge_spd_random ( n, key )
#
# Choose a random solution.
#
x1 = rng.random ( size = n )
#
# Compute the corresponding right hand side.
#
b = r8ge_mv ( n, n, a, x1 )
#
# Call the CG routine.
#
x2 = np.ones ( n )
x3 = r8ge_cg ( n, a, b, x2 )
#
# Compute the residual.
#
r = r8ge_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 = %d' % ( n ) )
print ( ' Norm of residual ||Ax-b|| = %g' % ( r_norm ) )
print ( ' Norm of error ||x1-x2|| = %g' % ( e_norm ) )
return
def r8ge_co ( n, a ):
#*****************************************************************************80
#
## r8ge_co() factors a R8GE matrix and estimates its condition number.
#
# 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.
#
# For the system A * X = B, relative perturbations in A and B
# of size EPSILON may cause relative perturbations in X of size
# EPSILON/RCOND.
#
# If RCOND is so small that the logical expression
# 1.0E+00 + rcond == 1.0E+00
# is true, then A may be singular to working precision. In particular,
# RCOND is zero if exact singularity is detected or the estimate
# underflows.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt.
#
# Reference:
#
# Dongarra, Bunch, Moler, Stewart,
# LINPACK User's Guide,
# SIAM, 1979
#
# Input:
#
# integer N, the order of the matrix A.
#
# real A(N,N), a matrix to be factored.
#
# Output:
#
# real A_LU(N,N), the LU factorization of the matrix.
#
# integer PIVOT(N), the pivot indices.
#
# real RCOND, an estimate of the reciprocal condition number of A.
#
# real Z(N), a work vector whose contents are usually unimportant.
# If A is close to a singular matrix, then Z is an approximate null vector
# in the sense that
# norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ).
#
import numpy as np
#
# Compute the L1 norm of A.
#
anorm = 0.0
for j in range ( 0, n ):
t = 0.0
for i in range ( 0, n ):
t = t + abs ( a[i,j] )
anorm = max ( anorm, t )
#
# Compute the LU factorization.
#
a_lu, pivot, info = r8ge_fa ( n, a )
#
# RCOND = 1 / ( norm(A) * (estimate of norm(inverse(A))) )
#
# estimate of norm(inverse(A)) = norm(Z) / norm(Y)
#
# where
# A * Z = Y
# and
# A' * Y = E
#
# The components of E are chosen to cause maximum local growth in the
# elements of W, where U'*W = E. The vectors are frequently rescaled
# to avoid overflow.
#
# Solve U' * W = E.
#
ek = 1.0
z = np.zeros ( n )
for k in range ( 0, n ):
if ( z[k] != 0.0 ):
ek = - r8_sign ( z[k] ) * abs ( ek )
if ( abs ( a_lu[k,k] ) < abs ( ek - z[k] ) ):
s = abs ( a_lu[k,k] ) / abs ( ek - z[k] )
for i in range ( 0, n ):
z[i] = s * z[i]
ek = s * ek
wk = ek - z[k]
wkm = - ek - z[k]
s = abs ( wk )
sm = abs ( wkm )
if ( a_lu[k,k] != 0.0 ):
wk = wk / a_lu[k,k]
wkm = wkm / a_lu[k,k]
else:
wk = 1.0
wkm = 1.0
if ( k + 1 <= n - 1 ):
for j in range ( k + 1, n ):
sm = sm + abs ( z[j] + wkm * a_lu[k,j] )
z[j] = z[j] + wk * a_lu[k,j]
s = s + abs ( z[j] )
if ( s < sm ):
t = wkm - wk
wk = wkm
for i in range ( k + 1, n ):
z[i] = z[i] + t * a_lu[k,i]
z[k] = wk
t = 0.0
for i in range ( 0, n ):
t = t + abs ( z[i] )
for i in range ( 0, n ):
z[i] = z[i] / t
#
# Solve L' * Y = W
#
for k in range ( n - 1, -1, -1 ):
for i in range ( k + 1, n ):
z[k] = z[k] + a_lu[i,k] * z[k]
t = abs ( z[k] )
if ( 1.0 < t ):
for i in range ( 0, n ):
z[i] = z[i] / t
l = pivot[k]
t = z[l]
z[l] = z[k]
z[k] = t
t = 0.0
for i in range ( 0, n ):
t = t + abs ( z[i] )
for i in range ( 0, n ):
z[i] = z[i] / t
ynorm = 1.0
#
# Solve L * V = Y.
#
for k in range ( 0, n ):
l = pivot[k]
t = z[l]
z[l] = z[k]
z[k] = t
for i in range ( k + 1, n ):
z[i] = z[i] + t * a_lu[i,k]
if ( 1.0 < abs ( z[k] ) ):
ynorm = ynorm / abs ( z[k] )
for i in range ( 0, n ):
z[i] = z[i] / abs ( z[k] )
s = 0.0
for i in range ( 0, n ):
s = s + abs ( z[i] )
for i in range ( 0, n ):
z[i] = z[i] / s
ynorm = ynorm / s
#
# Solve U * Z = V.
#
for k in range ( n - 1, -1, -1 ):
if ( abs ( a_lu[k,k] ) < abs ( z[k] ) ):
s = abs ( a_lu[k,k] ) / abs ( z[k] )
for i in range ( 0, n ):
z[i] = s * z[i]
ynorm = s * ynorm
if ( a_lu[k,k] != 0.0 ):
z[k] = z[k] / a_lu[k,k]
else:
z[k] = 1.0
for i in range ( 0, k ):
z[i] = z[i] - z[k] * a_lu[i,k]
#
# Normalize Z in the L1 norm.
#
t = 0.0
for i in range ( 0, n ):
t = t + abs ( z[i] )
for i in range ( 0, n ):
z[i] = z[i] / t
ynorm = ynorm / t
if ( anorm != 0.0E+00 ):
rcond = ynorm / anorm
else:
rcond = 0.0
return a_lu, pivot, rcond, z
def r8ge_co_test ( ):
#*****************************************************************************80
#
## r8ge_co_test() tests r8ge_co().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 4
x = 2.0
y = 3.0
print ( '' )
print ( 'r8ge_co_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_co estimates the condition number of an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
if ( i == j ):
a[i,j] = x + y
else:
a[i,j] = y
a_norm_l1 = 0.0
for j in range ( 0, n ):
t = 0.0
for i in range ( 0, n ):
t = t + abs ( a[i,j] )
a_norm_l1 = max ( a_norm_l1, t )
a_lu, pivot, info = r8ge_fa ( n, a )
a_inverse = r8ge_inverse ( n, a_lu, pivot )
a_inverse_norm_l1 = 0.0
for j in range ( 0, n ):
t = 0.0
for i in range ( 0, n ):
t = t + abs ( a_inverse[i,j] )
a_inverse_norm_l1 = max ( a_inverse_norm_l1, t )
cond_l1 = a_norm_l1 * a_inverse_norm_l1
print ( '' )
print ( ' The L1 condition number is %g' % ( cond_l1 ) )
#
# Factor the matrix.
#
a_lu, pivot, rcond, z = r8ge_co ( n, a )
print ( '' )
print ( ' The r8ge_co estimate is %g' % ( 1.0 / rcond ) )
return
def r8ge_det ( n, a_lu, pivot ):
#*****************************************************************************80
#
## r8ge_det(): determinant of a matrix factored by r8ge_fa() or r8ge_trf().
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
# LINPACK User's Guide,
# SIAM, 1979,
# ISBN13: 978-0-898711-72-1,
# LC: QA214.L56.
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A_LU(N,N), the LU factors from r8ge_fa
# or r8ge_trf.
#
# integer PIVOT(N), as computed by r8ge_fa or r8ge_trf.
#
# Output:
#
# real DET, the determinant of the matrix.
#
value = 1.0
for i in range ( 0, n ):
value = value * a_lu[i,i]
if ( pivot[i] != i ):
value = - value
return value
def r8ge_det_test ( ):
#*****************************************************************************80
#
## r8ge_det_test() tests r8ge_det().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 24 August 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 4
x = 2.0
y = 3.0
print ( '' )
print ( 'r8ge_det_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_det computes the determinant of an R8GE matrix.' )
#
# Set the matrix.
#
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
if ( i == j ):
a[i,j] = x + y
else:
a[i,j] = y
#
# Factor the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
#
# Compute the determinant.
#
det = r8ge_det ( n, a_lu, pivot )
exact = x ** ( n - 1 ) * ( x + n * y )
print ( '' )
print ( ' r8ge_det computes the determinant = %g' % ( det ) )
print ( ' Exact determinant = %g' % ( exact ) )
return
def r8ge_dif2 ( m, n ):
#*****************************************************************************80
#
## r8ge_dif2() returns the DIF2 matrix in R8GE 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:
#
# 06 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.
#
# Input:
#
# integer M, N, the order of the matrix.
#
# Output:
#
# real A(M,N), the matrix.
#
import numpy as np
a = r8ge_zeros ( m, n )
for j in range ( 0, n ):
for i in range ( 0, m ):
if ( j == i - 1 ):
a[i,j] = -1.0
elif ( j == i ):
a[i,j] = 2.0
elif ( j == i + 1 ):
a[i,j] = -1.0
return a
def r8ge_dif2_test ( ):
#*****************************************************************************80
#
## r8ge_dif2_test() tests r8ge_dif2().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_dif2_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_dif2 returns the second difference matrix.' )
a = r8ge_dif2 ( m, n )
r8ge_print ( m, n, a, ' DIF2 matrix:' )
return
def r8ge_dilu ( m, n, a ):
#*****************************************************************************80
#
## r8ge_dilu() produces the diagonal incomplete LU factor of an R8GE matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# 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 R8GE matrix.
#
# Output:
#
# real D(M), the DILU factor.
#
import numpy as np
d = np.zeros ( m, dtype = np.float64 )
for i in range ( 0, m ):
if ( i < n ):
d[i] = a[i,i]
mn = min ( m, n )
for i in range ( 0, mn ):
d[i] = 1.0 / d[i]
for j in range ( i + 1, mn ):
d[j] = d[j] - a[j,i] * d[i] * a[i,j]
return d
def r8ge_dilu_test ( ):
#*****************************************************************************80
#
## r8ge_dilu_test() tests r8ge_dilu().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
ncol = 3
nrow = 3
n = nrow * ncol
m = n
print ( '' )
print ( 'r8ge_dilu_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_dilu returns the DILU factors of an R8GE matrix.' )
print ( '' )
print ( ' Matrix rows M = %d' % ( m ) )
print ( ' Matrix columns N = %d' % ( n ) )
a = np.zeros ( [ m, n ] )
for i in range ( 0, nrow * ncol ):
for j in range ( 0, nrow * ncol ):
if ( i == j ):
a[i,j] = 4.0
elif ( i == j + 1 or i == j - 1 or i == j + nrow or i == j - nrow ):
a[i,j] = -1.0
else:
a[i,j] = 0.0
r8ge_print ( m, n, a, ' Matrix A:' )
#
# Compute the incomplete LU factorization.
#
d = r8ge_dilu ( m, n, a )
r8vec_print ( m, d, ' DILU factor:' )
return
def r8ge_fa ( n, a ):
#*****************************************************************************80
#
## r8ge_fa() performs a LINPACK style PLU factorization of 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.
#
# r8ge_fa is a simplified version of the LINPACK routine R8GEFA.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Dongarra, Bunch, Moler, Stewart,
# LINPACK User's Guide,
# SIAM, 1979
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A(N,N), the matrix to be factored.
#
# Output:
#
# real A_LU(N,N), an upper triangular matrix and
# the multipliers used to obtain it. The factorization
# can be written A = L * U, where L is a product of
# permutation and unit lower triangular matrices and U
# is upper triangular.
#
# integer PIVOT(N), a vector of pivot indices.
#
# integer INFO, singularity flag.
# 0, no singularity detected.
# nonzero, the factorization failed on the INFO-th step.
#
import numpy as np
a_lu = r8ge_zeros ( n, n )
for j in range ( 0, n ):
for i in range ( 0, n ):
a_lu[i,j] = a[i,j]
info = 0
pivot = np.zeros ( n, dtype = np.int32 )
for k in range ( 0, n - 1 ):
#
# Find L, the index of the pivot row.
#
l = k
for i in range ( k + 1, n ):
if ( abs ( a_lu[l,k] ) < abs ( a_lu[i,k] ) ):
l = i
pivot[k] = l
#
# If the pivot index is zero, the algorithm has failed.
#
if ( a_lu[l,k] == 0.0 ):
info = k
return a_lu, pivot, info
#
# Interchange rows L and K if necessary.
#
if ( l != k ):
t = a_lu[l,k]
a_lu[l,k] = a_lu[k,k]
a_lu[k,k] = t
#
# Normalize the values that lie below the pivot entry A(K,K).
#
for i in range ( k + 1, n ):
a_lu[i,k] = - a_lu[i,k] / a_lu[k,k]
#
# Row elimination with column indexing.
#
for j in range ( k + 1, n ):
if ( l != k ):
t = a_lu[l,j]
a_lu[l,j] = a_lu[k,j]
a_lu[k,j] = t
for i in range ( k + 1, n ):
a_lu[i,j] = a_lu[i,j] + a_lu[i,k] * a_lu[k,j]
pivot[n-1] = n - 1
if ( a_lu[n-1,n-1] == 0.0 ):
info = n - 1
return a_lu, pivot, info
def r8ge_fa_test01 ( ):
#*****************************************************************************80
#
## r8ge_fa_test01() tests r8ge_fa(), r8ge_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 10
print ( '' )
print ( 'r8ge_fa_test01' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_fa computes the LU factors of an R8GE matrix,' )
print ( ' r8ge_sl solves a factored R8GE system.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = r8ge_random ( n, n )
#
# Set the desired solution.
#
x = np.zeros ( n )
for i in range ( 0, n ):
x[i] = float ( i + 1 )
#
# Compute the corresponding right hand side.
#
b = r8ge_mv ( n, n, a, x )
#
# Factor the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_fa_test01(): Warning!' )
print ( ' r8ge_fa declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
return
#
# Solve the linear system.
#
job = 0
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution.
#
for i in range ( 0, n ):
x[i] = 1.0
#
# Compute the corresponding right hand side.
#
job = 0
b = r8ge_ml ( n, a_lu, pivot, x, job )
#
# Solve the system
#
job = 0
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution.
#
x = np.zeros ( n )
for i in range ( 0, n ):
x[i] = float ( i + 1 )
#
# Compute the corresponding right hand side.
#
job = 1
b = r8ge_ml ( n, a_lu, pivot, x, job )
#
# Solve the system
#
job = 1
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution of transposed system:' )
return
def r8ge_fa_test02 ( ):
#*****************************************************************************80
#
## r8ge_fa_test02() tests r8ge_fa(), r8ge_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 5
print ( '' )
print ( 'r8ge_fa_test02' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_fa computes the LU factors of an R8GE system,' )
print ( ' r8ge_sl solves a factored R8GE system.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = r8ge_random ( n, n )
r8ge_print ( n, n, a, ' The matrix:' )
#
# Set the desired solution.
#
x = np.zeros ( n )
for i in range ( 0, n ):
x[i] = float ( i + 1 )
#
# Compute the corresponding right hand side.
#
b = r8ge_mv ( n, n, a, x )
#
# Factor the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_fa_test02 - Warning!' )
print ( ' r8ge_fa declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
#
# Display the gory details.
#
r8ge_print ( n, n, a_lu, ' The compressed LU factors:' )
i4vec_print ( n, pivot, ' The pivot vector P:' )
#
# Solve the linear system.
#
job = 0
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution:' )
return
def r8ge_fs ( n, a, b ):
#*****************************************************************************80
#
## r8ge_fs() factors and solves a R8GE system.
#
# 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.
#
# r8ge_fs does not save the LU factors of the matrix, and hence cannot
# be used to efficiently solve multiple linear systems, or even to
# factor A at one time, and solve a single linear system at a later time.
#
# r8ge_fs uses partial pivoting, but no pivot vector is required.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A(N,N), the coefficient matrix of the linear system.
#
# real B(N), the right hand side of the linear system.
#
# Output:
#
# real X(N), the solution of the linear system.
#
import numpy as np
info = 0
x = b.copy ( )
for jcol in range ( 0, n ):
#
# Find the maximum element in column I.
#
piv = abs ( a[jcol,jcol] )
ipiv = jcol
for i in range ( jcol + 1, n ):
if ( piv < abs ( a[i,jcol] ) ):
piv = abs ( a[i,jcol] )
ipiv = i
if ( piv == 0.0 ):
info = jcol
return
#
# Switch rows JCOL and IPIV, and B.
#
if ( jcol != ipiv ):
for j in range ( 0, n ):
t = a[jcol,j]
a[jcol,j] = a[ipiv,j]
a[ipiv,j] = t
t = x[jcol]
x[jcol] = x[ipiv]
x[ipiv] = t
#
# Scale the pivot row.
#
t = a[jcol,jcol]
a[jcol,jcol] = 1.0
for k in range ( jcol + 1, n ):
a[jcol,k] = a[jcol,k] / t
x[jcol] = x[jcol] / t
#
# Use the pivot row to eliminate lower entries in that column.
#
for i in range ( jcol + 1, n ):
if ( a[i,jcol] != 0.0 ):
t = - a[i,jcol]
a[i,jcol] = 0.0
for k in range ( jcol + 1, n ):
a[i,k] = a[i,k] + t * a[jcol,k]
x[i] = x[i] + t * x[jcol]
#
# Back solve.
#
for jcol in range ( n - 1, 0, -1 ):
for k in range ( 0, jcol ):
x[k] = x[k] - a[k,jcol] * x[jcol]
return x
def r8ge_fs_test ( ):
#*****************************************************************************80
#
## r8ge_fs_test() tests r8ge_fs().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import platform
n = 10
print ( '' )
print ( 'r8ge_fs_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_fs factors and solves a linear system involving' )
print ( ' an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = r8ge_random ( n, n )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r8ge_mv ( n, n, a, x )
#
# Factor and solve the system.
#
x = r8ge_fs ( n, a, b )
r8vec_print ( n, x, ' Solution:' )
return
def r8ge_fss ( n, a, nb, b ):
#*****************************************************************************80
#
## r8ge_fss() factors and solves a R8GE system.
#
# 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.
#
# This routine does not save the LU factors of the matrix, and hence cannot
# be used to efficiently solve multiple linear systems, or even to
# factor A at one time, and solve a single linear system at a later time.
#
# This routine uses partial pivoting, but no pivot vector is required.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 23 June 2009
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A(N,N), the coefficient matrix of the linear system.
#
# integer NB, the number of right hand sides.
#
# real B(N,NB), the right hand side of the linear system.
#
# Output:
#
# real B(N,NB), the solution of the linear system.
#
import numpy as np
info = 0
for jcol in range ( 0, n ):
#
# Find the maximum element in column I.
#
piv = abs ( a[jcol,jcol] )
ipiv = jcol
for i in range ( jcol + 1, n ):
if ( piv < abs ( a[i,jcol] ) ):
piv = abs ( a[i,jcol] )
ipiv = i
if ( piv == 0.0 ):
info = jcol
return
#
# Switch rows JCOL and IPIV, and B.
#
if ( jcol != ipiv ):
for j in range ( 0, n ):
t = a[jcol,j]
a[jcol,j] = a[ipiv,j]
a[ipiv,j] = t
for j in range ( 0, nb ):
t = b[jcol,j]
b[jcol,j] = b[ipiv,j]
b[ipiv,j] = t
#
# Scale the pivot row.
#
t = a[jcol,jcol]
a[jcol,jcol] = 1.0
for k in range ( jcol + 1, n ):
a[jcol,k] = a[jcol,k] / t
for k in range ( 0, nb ):
b[jcol,k] = b[jcol,k] / t
#
# Use the pivot row to eliminate lower entries in that column.
#
for i in range ( jcol + 1, n ):
if ( a[i,jcol] != 0.0 ):
t = - a[i,jcol]
a[i,jcol] = 0.0
for k in range ( jcol + 1, n ):
a[i,k] = a[i,k] + t * a[jcol,k]
for k in range ( 0, nb ):
b[i,k] = b[i,k] + t * b[jcol,k]
#
# Back solve.
#
for j in range ( 0, nb ):
for jcol in range ( n - 1, 0, -1 ):
for k in range ( 0, jcol ):
b[k,j] = b[k,j] - a[k,jcol] * b[jcol,j]
return b
def r8ge_fss_test ( ):
#*****************************************************************************80
#
## r8ge_fss_test() tests r8ge_fss().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 10
nb = 3
print ( '' )
print ( 'r8ge_fss_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_fss factors and solves multiple linear systems' )
print ( ' associated with an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = r8ge_random ( n, n )
#
# Set the desired solutions.
#
x = np.zeros ( [ n, 3 ] )
for i in range ( 0, n ):
x[i,0] = 1.0
x[i,1] = i + 1
x[i,2] = ( i % 3 ) + 1
b = r8ge_mm ( n, n, nb, a, x )
#
# Factor and solve the system.
#
x = r8ge_fss ( n, a, nb, b )
r8ge_print ( n, nb, x, ' Solution:' )
return
def r8ge_hilbert ( m, n ):
#*****************************************************************************80
#
## r8ge_hilbert() returns the Hilbert matrix.
#
# Formula:
#
# A(I,J) = 1 / ( I + J - 1 )
#
# Example:
#
# N = 5
#
# 1/1 1/2 1/3 1/4 1/5
# 1/2 1/3 1/4 1/5 1/6
# 1/3 1/4 1/5 1/6 1/7
# 1/4 1/5 1/6 1/7 1/8
# 1/5 1/6 1/7 1/8 1/9
#
# Rectangular Properties:
#
# A is a Hankel matrix: constant along anti-diagonals.
#
# Square Properties:
#
# A is positive definite.
#
# A is symmetric: A' = A.
#
# Because A is symmetric, it is normal.
#
# Because A is normal, it is diagonalizable.
#
# A is totally positive.
#
# A is a Cauchy matrix.
#
# A is nonsingular.
#
# A is very ill-conditioned.
#
# The entries of the inverse of A are all integers.
#
# The sum of the entries of the inverse of A is N*N.
#
# The ratio of the absolute values of the maximum and minimum
# eigenvalues is roughly EXP(3.5*N).
#
# The determinant of the Hilbert matrix of order 10 is
# 2.16417... * 10^(-53).
#
# If the (1,1) entry of the 5 by 5 Hilbert matrix is changed
# from 1 to 24/25, the matrix is exactly singular. And there
# is a similar rule for larger Hilbert matrices.
#
# The family of matrices is nested as a function of N.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
# Reference:
#
# MD Choi,
# Tricks or treats with the Hilbert matrix,
# American Mathematical Monthly,
# Volume 90, 1983, pages 301-312.
#
# Robert Gregory, David Karney,
# Example 3.8,
# A Collection of Matrices for Testing Computational Algorithms,
# Wiley, New York, 1969, page 33,
# LC: QA263.G68.
#
# Nicholas Higham,
# Accuracy and Stability of Numerical Algorithms,
# Society for Industrial and Applied Mathematics, Philadelphia, PA,
# USA, 1996; section 26.1.
#
# Donald Knuth,
# The Art of Computer Programming,
# Volume 1, Fundamental Algorithms, Second Edition
# Addison-Wesley, Reading, Massachusetts, 1973, page 37.
#
# Morris Newman and John Todd,
# Example A13,
# The evaluation of matrix inversion programs,
# Journal of the Society for Industrial and Applied Mathematics,
# Volume 6, 1958, pages 466-476.
#
# Joan Westlake,
# Test Matrix A12,
# A Handbook of Numerical Matrix Inversion and Solution of Linear Equations,
# John Wiley, 1968.
#
# Input:
#
# integer M, N, the number of rows and columns of A.
#
# Output:
#
# real A(M,N), the matrix.
#
import numpy as np
a = np.zeros ( [ m, n ] )
for i in range ( 0, m ):
for j in range ( 0, n ):
a[i,j] = 1.0 / ( i + j + 1 )
return a
def r8ge_hilbert_test ( ):
#*****************************************************************************80
#
## r8ge_hilbert_test() tests r8ge_hilbert().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_hilbert_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_hilbert returns the Hilbert matrix.' )
a = r8ge_hilbert ( m, n )
r8ge_print ( m, n, a, ' Hilbert matrix:' )
return
def r8ge_hilbert_inverse ( n ):
#*****************************************************************************80
#
## r8ge_hilbert_inverse() returns the inverse of the Hilbert matrix.
#
# Formula:
#
# A(I,J) = (-1)^(I+J) * (N+I-1)! * (N+J-1)! /
# [ (I+J-1) * ((I-1)!*(J-1)!)^2 * (N-I)! * (N-J)! ]
#
# Example:
#
# N = 5
#
# 25 -300 1050 -1400 630
# -300 4800 -18900 26880 -12600
# 1050 -18900 79380 -117600 56700
# -1400 26880 -117600 179200 -88200
# 630 -12600 56700 -88200 44100
#
# Properties:
#
# A is symmetric.
#
# Because A is symmetric, it is normal, so diagonalizable.
#
# A is almost impossible to compute accurately by general routines
# that compute the inverse.
#
# A is integral.
#
# The sum of the entries of A is N^2.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of A.
#
# Output:
#
# real A(N,N), the inverse Hilbert matrix.
#
import numpy as np
a = np.zeros ( [ n, n ] )
#
# Set the (1,1) entry.
#
a[0,0] = n * n
#
# Define Row 1, Column J by recursion on Row 1 Column J-1
#
i = 0
for j in range ( 1, n ):
a[i,j] = - a[i,j-1] * float ( ( n + j ) * ( i + j ) * ( n - j ) ) \
/ float ( ( i + j + 1 ) * j * j )
#
# Define Row I by recursion on row I-1
#
for i in range ( 1, n ):
for j in range ( 0, n ):
a[i,j] = - a[i-1,j] * float ( ( n + i ) * ( i + j ) * ( n - i ) ) \
/ float ( ( i + j + 1 ) * i * i )
return a
def r8ge_hilbert_inverse_test ( ):
#*****************************************************************************80
#
## r8ge_hilbert_inverse_test() tests r8ge_hilbert_inverse().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
import platform
n = 4
print ( '' )
print ( 'r8ge_hilbert_inverse_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_hilbert_inverse computes the inverse of the' )
print ( ' Hilbert matrix, stored as an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
a = r8ge_hilbert ( n, n )
r8ge_print ( n, n, a, ' Matrix A:' )
b = r8ge_hilbert_inverse ( n )
r8ge_print ( n, n, b, ' Inverse matrix B:' )
#
# Check.
#
c = r8ge_mm ( n, n, n, a, b )
r8ge_print ( n, n, c, ' Product A * B:' )
return
def r8ge_house_axh ( n, a, v ):
#*****************************************************************************80
#
## r8ge_house_axh() computes A*H where H is a compact Householder matrix.
#
# Discussion:
#
# The Householder matrix H(V) is defined by
#
# H(V) = I - 2 * v * v' / ( v' * v )
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 24 April 2013
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of A.
#
# real A(N,N), the matrix to be postmultiplied.
#
# real V(N), a vector defining a Householder matrix.
#
# Output:
#
# real AH(N,N), the product A*H.
#
import numpy as np
vtv = 0.0
for i in range ( 0, n ):
vtv = vtv + v[i] ** 2
ah = np.zeros ( ( n, n ) )
for j in range ( 0, n ):
for i in range ( 0, n ):
ah[i,j] = a[i,j]
for k in range ( 0, n ):
ah[i,j] = ah[i,j] - 2.0 * a[i,k] * v[k] * v[j] / vtv
return ah
def r8ge_house_axh_test ( ):
#*****************************************************************************80
#
## r8ge_house_axh_test() tests r8ge_house_axh().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 February 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 5
print ( '' )
print ( 'r8ge_house_axh_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_house_axh multiplies a matrix A times a' )
print ( ' compact Householder matrix.' )
r8_lo = -5.0
r8_hi = +5.0
a = r8ge_random_ab ( n, n, r8_lo, r8_hi )
r8ge_print ( n, n, a, ' Matrix A:' )
#
# Request V, the compact form of the Householder matrix H
# such that H*A packs column 3 of A.
#
k = 3
km1 = k - 1
a_col = np.zeros ( n )
for i in range ( 0, n ):
a_col[i] = a[i,km1]
v = r8vec_house_column ( n, a_col, km1 )
r8vec_print ( n, v, ' Compact vector V so column 3 of H*A is packed:' )
h = r8ge_house_form ( n, v )
r8ge_print ( n, n, h, ' Householder matrix H:' )
#
# Compute A*H.
#
ah = r8ge_house_axh ( n, a, v )
r8ge_print ( n, n, ah, ' Indirect product A*H:' )
#
# Compare with a direct calculation.
#
ah = r8ge_mm ( n, n, n, a, h )
r8ge_print ( n, n, ah, ' Direct product A*H:' )
#
# Compute H*A to demonstrate packed column 3:
#
ha = r8ge_mm ( n, n, n, h, a )
r8ge_print ( n, n, ha, ' H*A should pack column 3:' )
return
def r8ge_house_form ( n, v ):
#*****************************************************************************80
#
## r8ge_house_form() constructs a Householder matrix from its compact form.
#
# Discussion:
#
# H(v) = I - 2 * v * v' / ( v' * v )
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 April 2013
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
#
# real V(N,1), the vector defining the Householder matrix.
#
# Output:
#
# real H(N,N), the Householder matrix.
#
import numpy as np
v_dot_v = 0.0
for i in range ( 0, n ):
v_dot_v = v_dot_v + v[i] * v[i]
c = - 2.0 / v_dot_v
h = np.zeros ( ( n, n ) )
for j in range ( 0, n ):
h[j,j] = 1.0
for i in range ( 0, n ):
h[i,j] = h[i,j] + c * v[i] * v[j]
return h
def r8ge_house_form_test ( ):
#*****************************************************************************80
#
## r8ge_house_form_test() tests r8ge_house_form().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 16 February 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 5
v = np.array ( ( 0.0, 0.0, 1.0, 2.0, 3.0 ) )
print ( '' )
print ( 'r8ge_house_form_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_house_form forms a Householder' )
print ( ' matrix from its compact form.' )
r8vec_print ( n, v, ' Compact vector form V:' )
h = r8ge_house_form ( n, v )
r8ge_print ( n, n, h, ' Householder matrix H:' )
return
def r8ge_identity ( m, n ):
#*****************************************************************************80
#
## r8ge_identity() copies the identity matrix to an R8GE matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of A.
#
# Output:
#
# real A(M,N), the identity matrix.
#
import numpy as np
a = np.zeros ( [ m, n ] )
for i in range ( 0, min ( m, n ) ):
a[i,i] = 1.0
return a
def r8ge_identity_test ( ):
#*****************************************************************************80
#
## r8ge_identity_test() tests r8ge_identity().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_identity_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_identity returns the identity matrix.' )
a = r8ge_identity ( m, n )
r8ge_print ( m, n, a, ' Identity matrix:' )
return
def r8ge_ilu ( m, n, a ):
#*****************************************************************************80
#
## r8ge_ilu() produces the incomplete LU factors of 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.
#
# The incomplete LU factors of the M by N matrix A are:
#
# L, an M by M unit lower triangular matrix,
# U, an M by N upper triangular matrix
#
# with the property that L and U are computed in the same way as
# the usual LU factors, except that, whenever an off diagonal element
# of the original matrix is zero, then the corresponding value of
# U is forced to be zero.
#
# This condition means that it is no longer the case that A = L*U.
#
# On the other hand, L and U will have a simple sparsity structure
# related to that of A. The incomplete LU factorization is generally
# used as a preconditioner in iterative schemes applied to sparse
# matrices. It is presented here merely for illustration.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# 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 R8GE matrix.
#
# Output:
#
# real L(M,M), the M by M unit lower triangular factor.
#
# real U(M,N), the M by N upper triangular factor.
#
import numpy as np
l = np.zeros ( [ m, m ] )
for i in range ( 0, m ):
l[i,i] = 1.0
u = a.copy ( )
for j in range ( 0, min ( m - 1, n ) ):
#
# Zero out the entries in column J, from row J+1 to M.
#
for i in range ( j + 1, m ):
if ( u[i,j] != 0.0 ):
l[i,j] = u[i,j] / u[j,j]
u[i,j] = 0.0
for k in range ( j + 1, n ):
if ( u[i,k] != 0.0 ):
u[i,k] = u[i,k] - l[i,j] * u[j,k]
return l, u
def r8ge_ilu_test ( ):
#*****************************************************************************80
#
## r8ge_ilu_test() tests r8ge_ilu().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
ncol = 3
nrow = 3
n = nrow * ncol
m = n
print ( '' )
print ( 'r8ge_ilu_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_ilu returns the ILU factors of an R8GE matrix.' )
print ( '' )
print ( ' Matrix rows M = %d' % ( m ) )
print ( ' Matrix columns N = %d' % ( n ) )
a = np.zeros ( [ m, n ] )
for i in range ( 0, nrow * ncol ):
for j in range ( 0, nrow * ncol ):
if ( i == j ):
a[i,j] = 4.0
elif ( i == j + 1 or i == j - 1 or i == j + nrow or i == j - nrow ):
a[i,j] = -1.0
else:
a[i,j] = 0.0
r8ge_print ( m, n, a, ' Matrix A:' )
#
# Compute the incomplete LU factorization.
#
l, u = r8ge_ilu ( m, n, a )
r8ge_print ( m, m, l, ' Factor L:' )
r8ge_print ( m, n, u, ' Factor U:' )
lu = r8ge_mm ( m, m, n, l, u )
r8ge_print ( m, n, lu, ' Product L*U:' )
return
def r8ge_indicator ( m, n ):
#*****************************************************************************80
#
## r8ge_indicator() sets an R8GE indicator matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of A.
#
# Output:
#
# real A(M,N), the matrix.
#
import numpy as np
fac = 10 ** ( i4_log_10 ( n ) + 1 )
a = r8ge_zeros ( m, n )
for j in range ( 0, n ):
for i in range ( 0, m ):
a[i,j] = float ( fac * ( i + 1 ) + ( j + 1 ) )
return a
def r8ge_indicator_test ( ):
#*****************************************************************************80
#
## r8ge_indicator_test() tests r8ge_indicator().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_indicator_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_indicator returns the indicator matrix.' )
a = r8ge_indicator ( m, n )
r8ge_print ( m, n, a, ' Indicator matrix:' )
return
def r8ge_inverse ( n, a_lu, pivot ):
#*****************************************************************************80
#
## r8ge_inverse() computes the inverse of a matrix factored by r8ge_fa.
#
# 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.
#
# r8ge_inverse is a simplified standalone version of the LINPACK routine
# R8GEDI.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix A.
#
# real A_LU(N,N), the factor information computed by r8ge_fa.
#
# integer PIVOT(N), the pivot vector from r8ge_fa.
#
# Output:
#
# real A_inverse(N,N), the inverse matrix.
#
import numpy as np
a_inverse = a_lu.copy ( )
#
# Compute Inverse(U).
#
for k in range ( 0, n ):
a_inverse[k,k] = 1.0 / a_inverse[k,k]
for i in range ( 0, k ):
a_inverse[i,k] = -a_inverse[i,k] * a_inverse[k,k]
for j in range ( k + 1, n ):
temp = a_inverse[k,j]
a_inverse[k,j] = 0.0
for i in range ( 0, k + 1 ):
a_inverse[i,j] = a_inverse[i,j] + a_inverse[i,k] * temp
#
# Form Inverse(U) * Inverse(L).
#
work = np.zeros ( n )
for k in range ( n - 2, -1, -1 ):
for i in range ( k + 1, n ):
work[i] = a_inverse[i,k]
a_inverse[i,k] = 0.0
for j in range ( k + 1, n ):
for i in range ( 0, n ):
a_inverse[i,k] = a_inverse[i,k] + a_inverse[i,j] * work[j]
if ( pivot[k] != k ):
for i in range ( 0, n ):
t = a_inverse[i,k]
a_inverse[i,k] = a_inverse[i,pivot[k]]
a_inverse[i,pivot[k]] = t
return a_inverse
def r8ge_inverse_test ( ):
#*****************************************************************************80
#
## r8ge_inverse_test() tests r8ge_inverse().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 14 March 2009
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 4
x = 2.0
y = 3.0
print ( '' )
print ( 'r8ge_inverse_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_inverse computes the inverse of an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = np.zeros ( [ n, n ] )
for j in range ( 0, n ):
for i in range ( 0, n ):
if ( i == j ):
a[i,i] = x + y
else:
a[i,j] = y
r8ge_print ( n, n, a, ' Matrix A:' )
#
# Factor and invert the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_inverse_test(): Fatal error!' )
print ( ' r8ge_fa reports the matrix is singular.' )
return
a_inverse = r8ge_inverse ( n, a_lu, pivot )
r8ge_print ( n, n, a_inverse, ' Inverse matrix B:' )
#
# Check.
#
c = r8ge_mm ( n, n, n, a, a_inverse )
r8ge_print ( n, n, c, ' Product matrix:' )
return
def r8ge_ml ( n, a_lu, pivot, x, job ):
#*****************************************************************************80
#
## r8ge_ml() computes A * x or A' * x, using r8ge_fa factors.
#
# 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.
#
# It is assumed that r8ge_fa has overwritten the original matrix
# information by LU factors. r8ge_ml is able to reconstruct the
# original matrix from the LU factor data.
#
# r8ge_ml allows the user to check that the solution of a linear
# system is correct, without having to save an unfactored copy
# of the matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A_LU(N,N), the LU factors from r8ge_fa.
#
# integer PIVOT(N), the pivot vector computed by r8ge_fa.
#
# real X(N), the vector to be multiplied.
#
# integer JOB, specifies the operation to be done:
# JOB = 0, compute A * x.
# JOB nonzero, compute A' * X.
#
# Output:
#
# real B(N), the result of the multiplication.
#
import numpy as np
b = x.copy ( )
if ( job == 0 ):
#
# Y = U * X.
#
for j in range ( 0, n ):
for i in range ( 0, j ):
b[i] = b[i] + a_lu[i,j] * b[j]
b[j] = a_lu[j,j] * b[j]
#
# B = PL * Y = PL * U * X = A * x.
#
for j in range ( n - 2, -1, -1 ):
for i in range ( j + 1, n ):
b[i] = b[i] - a_lu[i,j] * b[j]
k = pivot[j]
if ( k != j ):
t = b[k]
b[k] = b[j]
b[j] = t
else:
#
# Y = (PL)' * X:
#
for j in range ( 0, n - 1 ):
k = pivot[j]
if ( k != j ):
t = b[k]
b[k] = b[j]
b[j] = t
for i in range ( j + 1, n ):
b[j] = b[j] - b[i] * a_lu[i,j]
#
# B = U' * Y = ( PL * U )' * X = A' * X.
#
for i in range ( n - 1, -1, -1 ):
for j in range ( i + 1, n ):
b[j] = b[j] + b[i] * a_lu[i,j]
b[i] = b[i] * a_lu[i,i]
return b
def r8ge_ml_test ( ):
#*****************************************************************************80
#
## r8ge_ml_test() tests r8ge_ml().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
import platform
n = 10
print ( '' )
print ( 'r8ge_ml_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_ml computes A*x or A\'*X' )
print ( ' where A has been factored by r8ge_fa.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
for job in range ( 0, 2 ):
#
# Set the matrix.
#
a = r8ge_random ( n, n )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
if ( job == 0 ):
b = r8ge_mv ( n, n, a, x )
else:
b = r8ge_mtv ( n, n, a, x )
#
# Factor the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_ml_test(): Fatal error!' )
print ( ' r8ge_fa declares the matrix is singular!' )
print ( ' The value of INFO is %d,' % ( info ) )
return
#
# Now multiply factored matrix times solution to get right hand side again.
#
b2 = r8ge_ml ( n, a_lu, pivot, 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 r8ge_mm ( n1, n2, n3, a, b ):
#*****************************************************************************80
#
## r8ge_mm() multiplies two R8GE's.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 14 February 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N1, N2, N3, the order of the matrices.
#
# real A(N1,N2), B(N2,N3), the matrices to multiply.
#
# Output:
#
# real C(N1,N3), the product matrix C = A * B.
#
import numpy as np
c = r8ge_zeros ( n1, n3 )
for j in range ( 0, n3 ):
for i in range ( 0, n1 ):
for k in range ( 0, n2 ):
c[i,j] = c[i,j] + a[i,k] * b[k,j]
return c
def r8ge_mm_test ( ):
#*****************************************************************************80
#
## r8ge_mm_test() tests r8ge_mm().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 14 February 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n1 = 4
n2 = 3
n3 = n1
print ( '' )
print ( 'r8ge_mm_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_mm computes a matrix-matrix product C = A * B;' )
a = r8ge_zeros ( n1, n2 )
for i in range ( 0, n1 ):
for j in range ( 0, n2 ):
if ( j == 0 ):
a[i,j] = 1.0
elif ( i == 0 ):
a[i,j] = 0.0
else:
a[i,j] = a[i-1,j-1] + a[i-1,j]
b = np.transpose ( a )
c = r8ge_mm ( n1, n2, n3, a, b )
r8ge_print ( n1, n2, a, ' A:' )
r8ge_print ( n2, n3, b, ' B:' )
r8ge_print ( n1, n3, c, ' C = A*B:' )
return
def r8ge_mtm ( n1, n2, n3, a, b ):
#*****************************************************************************80
#
## r8ge_mtm() computes A' * B for two R8GE's.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N1, N2, N3, the order of the matrices.
#
# real A(N2,N1), B(N2,N3), the matrices to multiply.
#
# Output:
#
# real C(N1,N3), the product matrix C = A' * B.
#
import numpy as np
c = r8ge_zeros ( n1, n3 )
for j in range ( 0, n3 ):
for i in range ( 0, n1 ):
for k in range ( 0, n2 ):
c[i,j] = c[i,j] + a[k,i] * b[k,j]
return c
def r8ge_mtm_test ( ):
#*****************************************************************************80
#
## r8ge_mtm_test() tests r8ge_mtm().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 Augsut 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n1 = 4
n2 = 3
n3 = n1
print ( '' )
print ( 'r8ge_mtm_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_mtm computes a matrix-matrix product C = A\' * B;' )
a = r8ge_zeros ( n2, n1 )
for i in range ( 0, n2 ):
for j in range ( 0, n1 ):
if ( j == 0 ):
a[i,j] = 1.0
elif ( i == 0 ):
a[i,j] = 0.0
else:
a[i,j] = a[i-1,j-1] + a[i-1,j]
b = a
c = r8ge_mtm ( n1, n2, n3, a, b )
r8ge_print ( n2, n1, a, ' A:' )
r8ge_print ( n2, n3, b, ' B:' )
r8ge_print ( n1, n3, c, ' C = A\'*B:' )
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_mtv_test ( ):
#*****************************************************************************80
#
## r8ge_mtv_test() tests r8ge_mtv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 September 2015
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8ge_mtv_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_mtv computes a matrix product b=A\'*x for an R8GE matrix.' )
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 = r8ge_indicator ( m, n )
r8ge_print ( m, n, a, ' The matrix A:' )
x = r8vec_indicator1 ( m )
r8vec_print ( m, x, ' The vector x:' )
b = r8ge_mtv ( m, n, a, x )
r8vec_print ( n, b, ' The vector b=A\'*x:' )
return
def r8ge_mu ( m, n, a_lu, trans, pivot, x ):
#*****************************************************************************80
#
## r8ge_mu() computes A * x or A' * x, using r8ge_trf factors.
#
# 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.
#
# It is assumed that r8ge_trf has overwritten the original matrix
# information by PLU factors. r8ge_mu is able to reconstruct the
# original matrix from the PLU factor data.
#
# r8ge_mu allows the user to check that the solution of a linear
# system is correct, without having to save an unfactored copy
# of the matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum,
# Hammarling, McKenney, Ostrouchov, Sorensen,
# LAPACK User's Guide,
# Second Edition,
# SIAM, 1995.
#
# Input:
#
# integer M, the number of rows in the matrix.
#
# integer N, the number of columns in the matrix.
#
# real A_LU(M,N), the LU factors from r8ge_trf.
#
# character TRANS, specifies the form of the system of equations:
# 'N': A * x = b (No transpose)
# 'T': A'* X = B (Transpose)
# 'C': A'* X = B (Conjugate transpose = Transpose)
#
# integer PIVOT(*), the pivot vector computed by r8ge_trf.
#
# real X(*), the vector to be multiplied.
# For the untransposed case, X should have N entries.
# For the transposed case, X should have M entries.
#
# Output:
#
# real B(*), the result of the multiplication.
# For the untransposed case, B should have M entries.
# For the transposed case, B should have N entries.
#
import numpy as np
npiv = min ( m - 1, n )
mn_max = max ( m, n )
if ( trans == 'n' or trans == 'N' ):
#
# Y[MN] = U[MNxN] * X[N].
#
b = np.zeros ( m )
y = np.zeros ( n )
for j in range ( 0, n ):
for i in range ( 0, min ( j + 1, m ) ):
y[i] = y[i] + a_lu[i,j] * x[j]
#
# Z[M] = L[MxMN] * Y[MN] = L[MxMN] * U[MNxN] * X[N].
#
for i in range ( 0, m ):
if ( i <= n - 1 ):
b[i] = y[i]
else:
b[i] = 0.0
for j in range ( min ( m - 2, n - 1 ), -1, -1 ):
for i in range ( j + 1, m ):
b[i] = b[i] + a_lu[i,j] * y[j]
#
# B = P * Z = P * L * Y = P * L * U * X = A * x.
#
for j in range ( npiv - 1, -1, -1 ):
k = pivot[j]
if ( k != j ):
t = b[k]
b[k] = b[j]
b[j] = t
elif ( trans == 't' or trans == 'T' or trans == 'c' or trans == 'C' ):
b = np.zeros ( n )
#
# Y = transpose(P) * X:
#
for i in range ( 0, npiv ):
k = pivot[i]
if ( k != i ):
t = x[k]
x[k] = x[i]
x[i] = t
for i in range ( 0, n ):
if ( i <= m - 1 ):
b[i] = x[i]
else:
b[i] = 0.0
#
# Z = tranpose(L) * Y:
#
for j in range ( 0, min ( m - 1, n ) ):
for k in range ( j + 1, m ):
b[j] = b[j] + x[k] * a_lu[k,j]
#
# B = U' * Z.
#
for i in range ( m - 1, -1, -1 ):
for j in range ( i + 1, n ):
b[j] = b[j] + b[i] * a_lu[i,j]
if ( i <= n - 1 ):
b[i] = b[i] * a_lu[i,i]
#
# Now restore X.
#
for i in range ( npiv - 1, -1, -1 ):
k = pivot[i]
if ( k != i ):
t = x[k]
x[k] = x[i]
x[i] = t
else:
print ( '' )
print ( 'r8ge_mu(): Fatal error!' )
print ( ' Illegal value of TRANS = %c' % ( trans ) )
raise Exception ( 'r8ge_mu(): Fatal error!' )
return b
def r8ge_mu_test ( ):
#*****************************************************************************80
#
## r8ge_mu_test() tests r8ge_mu().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 3
print ( '' )
print ( 'r8ge_mu_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_mu computes A*x or A\'*X' )
print ( ' where A has been factored by r8ge_trf.' )
print ( '' )
print ( ' Matrix rows M = %d' % ( m ) )
print ( ' Matrix columns N = %d' % ( n ) )
for job in range ( 0, 2 ):
if ( job == 0 ):
trans = 'N'
else:
trans = 'T'
#
# Set the matrix.
#
amn = r8ge_random ( m, n )
if ( job == 0 ):
xn = r8vec_indicator1 ( n )
cm = r8ge_mv ( m, n, amn, xn )
else:
xm = r8vec_indicator1 ( m )
cn = r8ge_mtv ( m, n, amn, xm )
#
# Factor the matrix.
#
amn_lu, pivot, info = r8ge_trf ( m, n, amn )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_mu_test(): Fatal error!' )
print ( ' r8ge_trf declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
continue
#
# Now multiply the factored matrix times solution to get right hand side again.
#
if ( job == 0 ):
bm = r8ge_mu ( m, n, amn_lu, trans, pivot, xn )
r8vec2_print_some ( m, cm, bm, 10, ' A*x and PLU*x' )
else:
bn = r8ge_mu ( m, n, amn_lu, trans, pivot, xm )
r8vec2_print_some ( n, cn, bn, 10, ' A\'*x and (PLU)\'*x' )
print ( '' )
print ( ' Matrix is %d by %d' % ( n, m ) )
for job in range ( 0, 2 ):
if ( job == 0 ):
trans = 'N'
else:
trans = 'T'
#
# Set the matrix.
#
anm = r8ge_random ( n, m )
if ( job == 0 ):
xm = r8vec_indicator1 ( m )
cn = r8ge_mv ( n, m, anm, xm )
else:
xn = r8vec_indicator1 ( n )
cm = r8ge_mtv ( n, m, anm, xn )
#
# Factor the matrix.
#
anm_lu, pivot, info = r8ge_trf ( n, m, anm )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_mu_test(): Fatal error!' )
print ( ' r8ge_trf declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
continue
#
# Now multiply factored matrix times solution to get right hand side again.
#
if ( job == 0 ):
bn = r8ge_mu ( n, m, anm_lu, trans, pivot, xm )
r8vec2_print_some ( n, cn, bn, 10, ' A*x and PLU*x' )
else:
bm = r8ge_mu ( n, m, anm_lu, trans, pivot, xn )
r8vec2_print_some ( m, cm, bm, 10, ' A\'*x and (PLU)\'*x' )
return
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 r8ge_mv_test ( ):
#*****************************************************************************80
#
## r8ge_mv_test() tests r8ge_mv().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 September 2015
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_mv_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_mv computes a matrix product b=A*x for an R8GE matrix.' )
a = r8ge_indicator ( m, n )
r8ge_print ( m, n, a, ' The matrix A:' )
x = r8vec_indicator1 ( m )
r8vec_print ( m, x, ' The vector X:' )
b = r8ge_mv ( m, n, a, x )
r8vec_print ( n, b, ' The vector b=A*x:' )
return
def r8ge_orth_random ( n, key ):
#*****************************************************************************80
#
## r8ge_orth_random() returns a random orthogonal matrix.
#
# Properties:
#
# The inverse of A is equal to A'.
#
# A is orthogonal: A * A' = A' * A = I.
#
# Because A is orthogonal, it is normal: A' * A = A * A'.
#
# Columns and rows of A have unit Euclidean norm.
#
# Distinct pairs of columns of A are orthogonal.
#
# Distinct pairs of rows of A are orthogonal.
#
# The L2 vector norm of A*x = the L2 vector norm of x for any vector x.
#
# The L2 matrix norm of A*B = the L2 matrix norm of B for any matrix B.
#
# det ( A ) = +1 or -1.
#
# A is unimodular.
#
# All the eigenvalues of A have modulus 1.
#
# All singular values of A are 1.
#
# All entries of A are between -1 and 1.
#
# Discussion:
#
# Thanks to Eugene Petrov, B I Stepanov Institute of Physics,
# National Academy of Sciences of Belarus, for convincingly
# pointing out the severe deficiencies of an earlier version of
# this routine.
#
# Essentially, the computation involves saving the Q factor of the
# QR factorization of a matrix whose entries are normally distributed.
# However, it is only necessary to generate this matrix a column at
# a time, since it can be shown that when it comes time to annihilate
# the subdiagonal elements of column K, these (transformed) elements of
# column K are still normally distributed random values. Hence, there
# is no need to generate them at the beginning of the process and
# transform them K-1 times.
#
# For computational efficiency, the individual Householder transformations
# could be saved, as recommended in the reference, instead of being
# accumulated into an explicit matrix format.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 17 March 2015
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Pete Stewart,
# Efficient Generation of Random Orthogonal Matrices With an Application
# to Condition Estimators,
# SIAM Journal on Numerical Analysis,
# Volume 17, Number 3, June 1980, pages 403-409.
#
# Input:
#
# integer N, the order of A.
#
# integer KEY, a positive value that selects the data.
#
# Output:
#
# real A(N,N), the matrix.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( key )
#
# Start with A = the identity matrix.
#
a = np.zeros ( ( n, n ) )
for i in range ( 0, n ):
a[i,i] = 1.0
#
# Now behave as though we were computing the QR factorization of
# some other random matrix. Generate the N elements of the first column,
# compute the Householder matrix H1 that annihilates the subdiagonal elements,
# and set A := A * H1' = A * H.
#
# On the second step, generate the lower N-1 elements of the second column,
# compute the Householder matrix H2 that annihilates them,
# and set A := A * H2' = A * H2 = H1 * H2.
#
# On the N-1 step, generate the lower 2 elements of column N-1,
# compute the Householder matrix HN-1 that annihilates them, and
# and set A := A * H(N-1)' = A * H(N-1) = H1 * H2 * ... * H(N-1).
# This is our random orthogonal matrix.
#
x_col = np.zeros ( n )
for j in range ( 0, n - 1 ):
for i in range ( 0, j ):
x_col[i] = 0.0
for i in range ( j, n ):
x_col[i] = rng.standard_normal ( )
#
# Compute the vector V that defines a Householder transformation matrix
# H(V) that annihilates the subdiagonal elements of X.
#
v = r8vec_house_column ( n, x_col, j )
#
# Postmultiply the matrix A by H'(V) = H(V).
#
a = r8ge_house_axh ( n, a, v )
return a
def r8ge_orth_random_test ( ):
#*****************************************************************************80
#
## r8ge_orth_random_test() tests r8ge_orth_random().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 17 March 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'r8ge_orth_random_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_orth_random computes a random orthogonal matrix.' )
m = 5
n = m
key = 123456789
a = r8ge_orth_random ( n, key )
r8ge_print ( m, n, a, ' orth_random matrix:' )
return
def r8ge_spd_random ( n, key ):
#*****************************************************************************80
#
## r8ge_spd_random() returns a random symmetric positive definite matrix.
#
# Discussion:
#
# The matrix returned will have eigenvalues in the range [0,1].
#
# Properties:
#
# A is symmetric: A' = A.
#
# A is positive definite: 0 < x'*A*x for nonzero x.
#
# The eigenvalues of A will be real.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 March 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of A.
#
# integer KEY, a positive value that selects the data.
#
# Output:
#
# real A(N,N), the matrix.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( key )
#
# Get a random set of eigenvalues.
#
lam = rng.random ( size = n )
#
# Get a random orthogonal matrix Q.
#
q = r8ge_orth_random ( n, key )
#
# Set A = Q * Lambda * Q'.
#
a = np.zeros ( ( n, n ) )
for i in range ( 0, n ):
for j in range ( 0, n ):
for k in range ( 0, n ):
a[i,j] = a[i,j] + q[i,k] * lam[k] * q[j,k]
return a
def r8ge_spd_random_test ( ):
#*****************************************************************************80
#
## r8ge_spd_random_test() tests r8ge_spd_random().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 March 2015
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8ge_spd_random_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_spd_random computes the spd_random matrix.' )
n = 5
key = 123456789
a = r8ge_spd_random ( n, key )
r8ge_print ( n, n, a, ' spd_random matrix:' )
return
def r8ge_plu ( m, n, a ):
#*****************************************************************************80
#
## r8ge_plu() produces the PLU factors of 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.
#
# The PLU factors of the M by N matrix A are:
#
# P, an M by M permutation matrix P,
# L, an M by M unit lower triangular matrix,
# U, an M by N upper triangular matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# 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 R8GE matrix.
#
# Output:
#
# real P(M,M), the M by M permutation factor.
#
# real L(M,M), the M by M unit lower triangular factor.
#
# real U(M,N), the M by N upper triangular factor.
#
l = r8ge_identity ( m, m )
p = r8ge_identity ( m, m )
u = a.copy ( )
#
# On step J, find the pivot row and the pivot value.
#
for j in range ( 0, min ( m - 1, n ) ):
pivot_value = 0.0
pivot_row = -1
for i in range ( j, m ):
if ( pivot_value < abs ( u[i,j] ) ):
pivot_value = abs ( u[i,j] )
pivot_row = i
#
# If the pivot row is nonzero, swap rows J and PIVOT_ROW.
#
if ( pivot_row != -1 ):
for k in range ( 0, n ):
t = u[j,k]
u[j,k] = u[pivot_row,k]
u[pivot_row,k] = t
for k in range ( 0, m ):
t = l[j,k]
l[j,k] = l[pivot_row,k]
l[pivot_row,k] = t
for k in range ( 0, m ):
t = l[k,j]
l[k,j] = l[k,pivot_row]
l[k,pivot_row] = t
for k in range ( 0, m ):
t = p[k,j]
p[k,j] = p[k,pivot_row]
p[k,pivot_row] = t
#
# Zero out the entries in column J, from row J+1 to M.
#
for i in range ( j + 1, m ):
if ( u[i,j] != 0.0 ):
l[i,j] = u[i,j] / u[j,j]
u[i,j] = 0.0
for k in range ( j + 1, n ):
u[i,k] = u[i,k] - l[i,j] * u[j,k]
return p, l, u
def r8ge_plu_test ( ):
#*****************************************************************************80
#
## r8ge_plu_test() tests r8ge_plu().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_plu_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_plu returns the PLU factors of an R8GE matrix.' )
print ( '' )
print ( ' Matrix rows M = %d' % ( m ) )
print ( ' Matrix columns N = %d' % ( n ) )
a = r8ge_random ( m, n )
r8ge_print ( m, n, a, ' Matrix A:' )
#
# Compute the PLU factors.
#
p, l, u = r8ge_plu ( m, n, a )
r8ge_print ( m, m, p, ' Factor P:' )
r8ge_print ( m, m, l, ' Factor L:' )
r8ge_print ( m, n, u, ' Factor U:' )
lu = r8ge_mm ( m, m, n, l, u )
plu = r8ge_mm ( m, m, n, p, lu )
r8ge_print ( m, n, plu, ' Product P*L*U:')
return
def r8ge_poly ( n, a ):
#*****************************************************************************80
#
## r8ge_poly() computes the characteristic polynomial of 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:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A(N,N), the R8GE matrix.
#
# Output:
#
# real P(1:N+1), the coefficients of the characteristic
# polynomial of A. P(I+1) contains the coefficient of X**I.
#
import numpy as np
#
# Initialize WORK1 to the identity matrix.
#
work1 = r8ge_identity ( n, n )
p = np.zeros ( n + 1 )
p[n] = 1.0
for order in range ( n - 1, -1, -1 ):
#
# Work2 = A * WORK1.
#
work2 = r8ge_mm ( n, n, n, a, work1 )
#
# Take the trace.
#
trace = 0.0
for i in range ( 0, n ):
trace = trace + work2[i,i]
#
# P(ORDER) = -Trace ( WORK2 ) / ( N - ORDER )
#
p[order] = - trace / ( n - order )
#
# WORK1 := WORK2 + P(ORDER) * Identity.
#
work1 = work2.copy ( )
for i in range ( 0, n ):
work1[i,i] = work1[i,i] + p[order]
return p
def r8ge_poly_test ( ):
#*****************************************************************************80
#
## r8ge_poly_test() tests r8ge_poly().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 12
print ( '' )
print ( 'r8ge_poly_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_poly computes the characteristic polynomial of an R8GE matrix.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
p_true = np.array ( [ \
1.0, -23.0, 231.0, -1330.0, 4845.0, \
-11628.0, 18564.0, -19448.0, 12870.0, -5005.0, \
1001.0, -78.0, 1.0 ] )
#
# Set the matrix.
#
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
a[i,j] = float ( min ( i, j ) + 1 )
#
# Get the characteristic polynomial.
#
p = r8ge_poly ( n, a )
#
# Compare.
#
r8vec2_print_some ( n+1, p, p_true, 10, 'I, P(I), True P(I)' )
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 ( ' Python version: %s' % ( platform.python_version ( ) ) )
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 r8ge_random ( m, n ):
#*****************************************************************************80
#
## r8ge_random() randomizes 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 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.
#
# Output:
#
# real A(M,N), the R8GE matrix.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
a = rng.random ( size = [ m, n ] )
return a
def r8ge_random_test ( ):
#*****************************************************************************80
#
## r8ge_random_test() tests r8ge_random().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_random_test():' )
print ( ' r8ge_random() computes a random R8GE.' )
print ( '' )
print ( ' 0 <= X <= 1' )
v = r8ge_random ( m, n )
r8ge_print ( m, n, v, ' Random R8GE:' )
return
def r8ge_random_ab ( m, n, a, b ):
#*****************************************************************************80
#
## r8ge_random_ab() returns a scaled pseudorandom R8GE.
#
# 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:
#
# 20 August 2015
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Paul Bratley, Bennett Fox, Linus Schrage,
# A Guide to Simulation,
# Second Edition,
# Springer, 1987,
# ISBN: 0387964673,
# LC: QA76.9.C65.B73.
#
# Bennett Fox,
# Algorithm 647:
# Implementation and Relative Efficiency of Quasirandom
# Sequence Generators,
# ACM Transactions on Mathematical Software,
# Volume 12, Number 4, December 1986, pages 362-376.
#
# Pierre L'Ecuyer,
# Random Number Generation,
# in Handbook of Simulation,
# edited by Jerry Banks,
# Wiley, 1998,
# ISBN: 0471134031,
# LC: T57.62.H37.
#
# Peter Lewis, Allen Goodman, James Miller,
# A Pseudo-Random Number Generator for the System/360,
# IBM Systems Journal,
# Volume 8, Number 2, 1969, pages 136-143.
#
# Input:
#
# integer M, N, the number of rows and columns in the array.
#
# real A, B, the range of the pseudorandom values.
#
# Output:
#
# real R(M,N), an array of random values between 0 and 1.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
r = a + ( b - a ) * rng.random ( size = [ m, n ] )
return r
def r8ge_random_ab_test ( ):
#*****************************************************************************80
#
## r8ge_random_ab_test() tests r8ge_random_ab().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 20 August 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
m = 5
n = 4
a = -1.0
b = +5.0
print ( '' )
print ( 'r8ge_random_ab_test():' )
print ( ' r8ge_random_ab computes a random R8GE.' )
print ( '' )
print ( ' %g <= X <= %g' % ( a, b ) )
v = r8ge_random_ab ( m, n, a, b )
r8ge_print ( m, n, v, ' Random R8GE:' )
return
def r8ge_res ( m, n, a, x, b ):
#*****************************************************************************80
#
## r8ge_res() computes the residual vector for an R8GE system.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns of
# the matrix. M and N must be positive.
#
# real A(M,N), the original, UNFACTORED R8GE matrix.
#
# real X(N), the estimated solution.
#
# real B(M), the right hand side vector.
#
# Output:
#
# real R(M), the residual vector, b - A * x.
#
import numpy as np
r = np.zeros ( m )
for i in range ( 0, m ):
r[i] = b[i]
for j in range ( 0, n ):
r[i] = r[i] - a[i,j] * x[j]
return r
def r8ge_res_test ( ):
#*****************************************************************************80
#
## r8ge_res_test() tests r8ge_res().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 04 February 2016
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8ge_res_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_res computes b-A*x, where A is an R8GE matrix.' )
print ( ' We check three cases, MN.' )
for i in range ( 0, 3 ):
if ( i == 0 ):
m = 3
n = 5
elif ( i == 1 ):
m = 5
n = 5
elif ( i == 2 ):
m = 5
n = 3
a = r8ge_random ( m, n )
x = r8vec_indicator1 ( n )
b = r8ge_mv ( m, n, a, x )
r = r8ge_res ( m, n, a, x, b )
r8vec_print ( m, r, ' Residual A*x-b:' )
return
def r8ge_sl ( n, a_lu, pivot, b, job ):
#*****************************************************************************80
#
## r8ge_sl() solves a system factored by r8ge_fa().
#
# 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.
#
# r8ge_sl is a simplified version of the LINPACK routine R8GESL.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 July 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A_LU(N,N), the LU factors from r8ge_fa.
#
# integer PIVOT(N), the pivot vector from r8ge_fa.
#
# real B(N), the right hand side vector.
#
# integer JOB, specifies the operation.
# 0, solve A * x = b.
# nonzero, solve A' * x = b.
#
# Output:
#
# real X(N), the solution vector.
#
import numpy as np
x = np.zeros ( n )
for i in range ( 0, n ):
x[i] = b[i]
#
# Solve A * x = b.
#
if ( job == 0 ):
#
# Solve PL * Y = B.
#
for k in range ( 0, n - 1 ):
l = pivot[k]
if ( l != k ):
t = x[l]
x[l] = x[k]
x[k] = t
for i in range ( k + 1, n ):
x[i] = x[i] + a_lu[i,k] * x[k]
#
# Solve U * X = Y.
#
for k in range ( n - 1, -1, -1 ):
x[k] = x[k] / a_lu[k,k]
for i in range ( 0, k ):
x[i] = x[i] - a_lu[i,k] * x[k]
#
# Solve A' * X = B.
#
else:
#
# Solve U' * Y = B.
#
for k in range ( 0, n ):
for i in range ( 0, k ):
x[k] = x[k] - x[i] * a_lu[i,k]
x[k] = x[k] / a_lu[k,k]
#
# Solve ( PL )' * X = Y.
#
for k in range ( n - 2, -1, -1 ):
for i in range ( k + 1, n ):
x[k] = x[k] + x[i] * a_lu[i,k]
l = pivot[k]
if ( l != k ):
t = x[l]
x[l] = x[k]
x[k] = t
return x
def r8ge_sl_test ( ):
#*****************************************************************************80
#
## r8ge_sl_test() tests r8ge_sl().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 March 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 10
print ( '' )
print ( 'r8ge_sl_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_sl solves a linear system A*x=b that was factored' )
print ( ' by r8ge_fa().' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the matrix.
#
a = r8ge_random ( n, n )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r8ge_mv ( n, n, a, x )
#
# Factor the matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_sl_test(): Fatal error!' )
print ( ' r8ge_fa() declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
return
#
# Solve the linear system.
#
job = 0
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution.
#
x = np.ones ( n )
#
# Compute the corresponding right hand side.
#
job = 0
b = r8ge_ml ( n, a_lu, pivot, x, job )
#
# Solve the system
#
job = 0
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution:' )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
job = 1
b = r8ge_ml ( n, a_lu, pivot, x, job )
#
# Solve the system
#
job = 1
x = r8ge_sl ( n, a_lu, pivot, b, job )
r8vec_print ( n, x, ' Solution of transposed system:' )
return
def r8ge_sl_it ( n, a, a_lu, pivot, b, job, x ):
#*****************************************************************************80
#
## r8ge_sl_it() applies one step of iterative refinement following r8ge_sl.
#
# 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.
#
# It is assumed that:
#
# * the original matrix A has been factored by r8ge_fa
# * the linear system A * x = b has been solved once by r8ge_sl.
#
# (Actually, it is not necessary to solve the system once using r8ge_sl.
# You may simply supply the initial estimated solution X = 0.)
#
# Each time this routine is called, it will compute the residual in
# the linear system, apply one step of iterative refinement, and
# add the computed correction to the current solution.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
# N must be positive.
#
# real A(N,N), the original, UNFACTORED R8GE matrix.
#
# real A_LU(N,N), the LU factors from r8ge_fa.
#
# integer PIVOT(N), the pivot vector from r8ge_fa.
#
# real B(N), the right hand side vector.
#
# integer JOB, specifies the operation.
# 0, solve A*X=B.
# nonzero, solve A'*X=B.
#
# real X(N), an estimate of the solution of A * x = b.
#
# Output:
#
# real X(N), an improved estimate of the solution.
#
# real DX(N), contains the correction terms added to X.
#
#
# Compute the residual vector.
#
r = r8ge_res ( n, n, a, x, b )
#
# Solve A * dx = r
#
dx = r8ge_sl ( n, a_lu, pivot, r, job )
#
# Add dx to x.
#
for i in range ( 0, n ):
x[i] = x[i] + dx[i]
return x, dx
def r8ge_sl_it_test ( ):
#*****************************************************************************80
#
## r8ge_sl_it_test() tests r8ge_sl_it().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 6
print ( '' )
print ( 'r8ge_sl_it_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_sl_it applies one step of iterative' )
print ( ' refinement to a r8ge_sl solution.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
#
# Set the coefficient matrix.
#
a = r8ge_hilbert_inverse ( n )
#
# Set the right hand side b.
#
b = np.zeros ( n )
b[n-1] = 1.0
#
# Compute the factored coefficient matrix.
#
a_lu, pivot, info = r8ge_fa ( n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_sl_it_test(): Fatal error!' )
print ( ' r8ge_fa() declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
return
#
# Solve the system.
#
job = 1
x = r8ge_sl ( n, a_lu, pivot, b, job )
#
# Compute and print the residual.
#
r = r8ge_res ( n, n, a, x, b )
r8vec2_print_some ( n, x, r, 10, ' i, x, b-A*x' )
#
# Take a few steps of iterative refinement.
#
for j in range ( 1, 6 ):
print ( '' )
print ( 'Iterative refinement step %d' % ( j ) )
#
# Improve the solution.
#
x_new, r = r8ge_sl_it ( n, a, a_lu, pivot, b, job, x )
r8vec_print_some ( n, r, 10, ' I, DX:' )
#
# Compute and print the residual.
#
r = r8ge_res ( n, n, a, x_new, b )
r8vec2_print_some ( n, x_new, r, 10, ' i, x, b-A*x' )
x = x_new.copy ( )
return
def r8ge_to_r8gb ( m, n, ml, mu, a ):
#*****************************************************************************80
#
## r8ge_to_r8gb() copies a R8GE matrix to a R8GB matrix.
#
# Discussion:
#
# It usually doesn't make sense to try to store a general matrix
# in a band matrix format. You can always do it, but it will take
# more space, unless the general matrix is actually banded.
#
# The purpose of this routine is to allow a user to set up a
# banded matrix in the easy-to-use general format, and have this
# routine take care of the compression of the data into general
# format. All the user has to do is specify the bandwidths.
#
# Note that this routine "believes" what the user says about the
# bandwidth. It will assume that all entries in the general matrix
# outside of the bandwidth are zero.
#
# The original M by N matrix is "collapsed" downward, so that diagonals
# become rows of the storage array, while columns are preserved. The
# collapsed array is logically 2*ML+MU+1 by N.
#
# The R8GB storage format is for an M by N banded matrix, with lower
# bandwidth ML and upper bandwidth MU. Storage includes room for ML
# extra superdiagonals, which may be required to store nonzero entries
# generated during Gaussian elimination.
#
# LINPACK and LAPACK band storage requires that an extra ML
# superdiagonals be supplied to allow for fillin during Gauss
# elimination. Even though a band matrix is described as
# having an upper bandwidth of MU, it effectively has an
# upper bandwidth of MU+ML. This routine will copy nonzero
# values it finds in these extra bands, so that both unfactored
# and factored matrices can be handled.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 02 July 2016
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum,
# Hammarling, McKenney, Ostrouchov, Sorensen,
# LAPACK User's Guide,
# Second Edition,
# SIAM, 1995.
#
# Dongarra, Bunch, Moler, Stewart,
# LINPACK User's Guide,
# SIAM, 1979.
#
# Input:
#
# integer M, the number of rows of the matrices.
# M must be positive.
#
# integer N, the number of columns of the matrices.
# N must be positive.
#
# integer ML, MU, the lower and upper bandwidths of the matrix.
# ML and MU must be nonnegative, and no greater than min(M,N)-1.
#
# real A(M,N), the R8GE matrix.
#
# Output:
#
# real B(2*ML+MU+1,N), the R8GB matrix.
#
import numpy as np
b = np.zeros ( [ 2*ml+mu+1, n ] )
for i in range ( 1, m + 1 ):
jlo = max ( i - ml, 1 )
jhi = min ( i + mu, n )
for j in range ( jlo, jhi + 1 ):
b[ml+mu+i-j,j-1] = a[i-1,j-1]
return b
def r8ge_to_r8gb_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8gb_test() tests r8ge_to_r8gb().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 02 July 2016
#
# Author:
#
# John Burkardt
#
import platform
from r8gb import r8gb_print
m = 5
n = 5
ml = 2
mu = 1
print ( '' )
print ( 'r8ge_to_r8gb_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_to_r8gb converts an R8GE matrix to R8GB format.' )
print ( '' )
print ( ' Matrix order M = %d, N = %d' % ( m, n ) )
print ( ' R8GB bands ML = %d, MU = %d' % ( ml, mu ) )
a = r8ge_random ( m, n )
r8ge_print ( m, n, a, ' The random R8GE matrix:' )
b = r8ge_to_r8gb ( m, n, ml, mu, a )
r8gb_print ( m, n, ml, mu, b, ' The R8GB matrix:' )
return
def r8ge_to_r8lt ( m, n, a_ge ):
#*****************************************************************************80
#
## r8ge_to_r8lt() copies an R8GE matrix to an R8LT matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# The R8LT storage format is used for an M by N lower triangular
# matrix. The format stores all M*N entries, even those which are zero.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# real A_GE(M,N), the R8GE matrix.
#
# Output:
#
# real A_LT(M,N), the R8LT matrix.
#
import numpy as np
a_lt = np.zeros ( [ m, n ] )
for j in range ( 0, n ):
for i in range ( j, m ):
a_lt[i,j] = a_ge[i,j]
return a_lt
def r8ge_to_r8lt_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8lt_test() tests r8ge_to_r8lt().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 22 August 2015
#
# Author:
#
# John Burkardt
#
import platform
from r8lt import r8lt_print
m = 5
n = 4
print ( '' )
print ( 'r8ge_to_r8lt_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_to_r8lt converts an R8GE matrix to R8LT format.' )
a_ge = r8ge_random ( m, n )
r8ge_print ( m, n, a_ge, ' The random R8GE matrix:' )
a_lt = r8ge_to_r8lt ( m, n, a_ge )
r8lt_print ( m, n, a_lt, ' The R8LT matrix:' )
return
def r8ge_to_r8po ( n, a ):
#*****************************************************************************80
#
## r8ge_to_r8po() copies an R8GE matrix to an R8PO matrix.
#
# Discussion:
#
# The R8PO format assumes the matrix is square and symmetric; it is also
# typically assumed that the matrix is positive definite. These are not
# required here. The copied R8PO matrix simply zeros out the lower triangle
# of the R8GE matrix.
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# The R8PO storage format is used for a symmetric positive definite
# matrix and its inverse. (The Cholesky factor of an R8PO matrix is an
# upper triangular matrix, so it will be in R8GE storage format.)
#
# Only the diagonal and upper triangle of the square array are used.
# This same storage scheme is used when the matrix is factored by
# R8PO_fa, or inverted by R8PO_inverse. For clarity, the lower triangle
# is set to zero.
#
# R8PO storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 02 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
#
# real A(N,N), the R8GE matrix.
#
# Output:
#
# real B(N,N), the R8PO matrix.
#
import numpy as np
b = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( i, n ):
b[i,j] = a[i,j]
return b
def r8ge_to_r8po_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8po_test() tests r8ge_to_r8po().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 12 August 2022
#
# Author:
#
# John Burkardt
#
n = 5
print ( '' )
print ( 'r8ge_to_r8po_test():' )
print ( ' r8ge_to_r8po() converts an R8GE matrix to R8PO format.' )
print ( '' )
print ( ' Matrix order N = %d' % ( n ) )
a = r8ge_random ( n, n )
print ( '' )
print ( ' The random R8GE matrix:' )
print ( a )
b = r8ge_to_r8po ( n, a )
print ( '' )
print ( ' The R8PO matrix:' )
print ( b )
return
def r8ge_to_r8pp ( n, a ):
#*****************************************************************************80
#
## r8ge_to_r8pp() copies an R8GE matrix to an R8PO matrix.
#
# Discussion:
#
# The R8PO format assumes the matrix is square and symmetric; it is also
# typically assumed that the matrix is positive definite. These are not
# required here. The copied R8PO matrix simply zeros out the lower triangle
# of the R8GE matrix.
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# The R8PP storage format is appropriate for a symmetric positive
# definite matrix. Only the upper triangle of the matrix is stored,
# by successive partial columns, in an array of length (N*(N+1))/2,
# which contains (A11,A12,A22,A13,A23,A33,A14,...,ANN)
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 12 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix.
#
# real A(N,N), the R8GE matrix.
#
# Output:
#
# real B(N*(N+1)/2), the R8PP matrix.
#
import numpy as np
b = np.zeros ( (n*(n+1)//2) )
k = 0
for j in range ( 0, n ):
for i in range ( 0, j + 1 ):
b[k] = a[i,j]
k = k + 1
return b
def r8ge_to_r8pp_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8pp_test() tests r8ge_to_r8pp().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 12 August 2022
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
n = 5
seed = 123456789
print ( '' )
print ( 'r8ge_to_r8pp_test():' )
print ( ' r8ge_to_r8pp() converts an R8GE matrix to R8PP format.' )
print ( '' )
print ( ' Matrix order N = ', n )
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
a[i,j] = 2.0 * min ( i, j ) - 1.0
print ( '' )
print ( ' The positive definite symmetric R8GE matrix:' )
print ( a )
b = r8ge_to_r8pp ( n, a )
print ( '' )
print ( ' The RPP matrix:' )
print ( b )
return
def r8ge_to_r8ut ( m, n, a_ge ):
#*****************************************************************************80
#
## r8ge_to_r8ut() copies an R8GE matrix to an R8UT matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# The R8UT storage format is used for an M by N upper triangular
# matrix. The format stores all M*N entries, even those which are zero.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 18 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# real A_GE(M,N), the R8GE matrix.
#
# Output:
#
# real A_UT(M,N), the R8UT matrix.
#
import numpy as np
a_ut = np.zeros ( [ m, n ] )
for j in range ( 0, n ):
for i in range ( 0, min ( j + 1, m ) ):
a_ut[i,j] = a_ge[i,j]
return a_ut
def r8ge_to_r8ut_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8ut_test() tests r8ge_to_r8ut().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 18 August 2015
#
# Author:
#
# John Burkardt
#
import platform
from r8ut import r8ut_print
m = 5
n = 4
print ( '' )
print ( 'r8ge_to_r8ut_test():' )
print ( ' r8ge_to_r8ut() converts an R8GE matrix to R8UT format.' )
a_ge = r8ge_random ( m, n )
r8ge_print ( m, n, a_ge, ' The random R8GE matrix:' )
a_ut = r8ge_to_r8ut ( m, n, a_ge )
r8ut_print ( m, n, a_ut, ' The R8UT matrix:' )
return
def r8ge_to_r8vec ( m, n, a ):
#*****************************************************************************80
#
## r8ge_to_r8vec() copies an R8GE matrix to an R8VEC.
#
# Discussion:
#
# In C++ and FORTRAN, this routine is not really needed. In MATLAB,
# a data item carries its dimensionality implicitly, and so cannot be
# regarded sometimes as a vector and sometimes as an array.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns in the array.
#
# real A(M,N), the array to be copied.
#
# Output:
#
# real X(M*N), the vector.
#
import numpy as np
x = np.zeros ( m * n )
k = 0
for j in range ( 0, n ):
for i in range ( 0, m ):
x[k] = a[i,j]
k = k + 1
return x
def r8ge_to_r8vec_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8vec_test() tests r8ge_to_r8vec().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 February 2016
#
# Author:
#
# John Burkardt
#
import platform
m = 4
n = 3
print ( '' )
print ( 'r8ge_to_r8vec_test():' )
print ( ' r8ge_to_r8vec() converts an R8GE matrix to an R8VEC vector.' )
a_r8ge = r8ge_indicator ( m, n )
r8ge_print ( m, n, a_r8ge, ' R8GE matrix:' )
a_r8vec = r8ge_to_r8vec ( m, n, a_r8ge )
r8vec_print ( m * n, a_r8vec, ' Corresponding R8VEC vector:' )
return
def r8ge_to_r8vm ( m, n, a_ge ):
#*****************************************************************************80
#
## r8ge_to_r8vm() copies an R8GE matrix to an R8VM matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# The R8VM storage format is used for an M by N Vandermonde matrix.
# An M by N Vandermonde matrix is defined by the values in its second
# row, which will be written here as X(1:N). The matrix has a first
# row of 1's, a second row equal to X(1:N), a third row whose entries
# are the squares of the X values, up to the M-th row whose entries
# are the (M-1)th powers of the X values. The matrix can be stored
# compactly by listing just the values X(1:N).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 24 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# real A_GE(M,N), the R8GE matrix.
#
# Output:
#
# real A_VM(N), the R8VM matrix.
#
import numpy as np
a_vm = np.zeros ( n )
for j in range ( 0, n ):
a_vm[j] = a_ge[1,j]
return a_vm
def r8ge_to_r8vm_test ( ):
#*****************************************************************************80
#
## r8ge_to_r8vm_test() tests r8ge_to_r8vm().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 02 July 2016
#
# Author:
#
# John Burkardt
#
from r8vm import r8vm_print
m = 5
n = 4
print ( '' )
print ( 'r8ge_to_r8vm_test():' )
print ( ' r8ge_to_r8vm() converts an R8GE matrix to R8VM format.' )
a_ge = r8ge_random ( m, n )
r8ge_print ( m, n, a_ge, ' The random R8GE matrix:' )
a_vm = r8ge_to_r8vm ( m, n, a_ge )
r8vm_print ( m, n, a_vm, ' The R8VM matrix:' )
return
def r8ge_transpose ( m, n, a ):
#*****************************************************************************80
#
## r8ge_transpose() makes a transposed copy of an R8GE matrix.
#
# Discussion:
#
# The R8GE storage format is used for a general M by N matrix. A storage
# space is made for each entry. The two dimensional logical
# 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.
#
# R8GE storage is used by LINPACK and LAPACK.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
#
# real A(M,N), the matrix to be copied.
#
# Output:
#
# real B(N,M), a copy of the transposed matrix.
#
import numpy as np
b = r8ge_zeros ( n, m )
for j in range ( 0, n ):
for i in range ( 0, m ):
b[j,i] = a[i,j]
return b
def r8ge_transpose_test ( ):
#*****************************************************************************80
#
## r8ge_transpose_test() tests r8ge_transpose().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 05 February 2016
#
# Author:
#
# John Burkardt
#
m = 5
n = 4
print ( '' )
print ( 'r8ge_transpose_test():' )
print ( ' r8ge_transpose() makes a transposed copy of an R8GE matrix.' )
a = r8ge_indicator ( m, n )
r8ge_print ( m, n, a, ' Indicator matrix A:' )
b = r8ge_transpose ( m, n, a )
r8ge_print ( n, m, b, ' B = A\':' )
return
def r8ge_transpose_print ( m, n, a, title ):
#*****************************************************************************80
#
## r8ge_transpose_print() prints an R8GE matrix, transposed.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 August 2014
#
# 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_transpose_print_some ( m, n, a, 0, 0, m - 1, n - 1, title )
return
def r8ge_transpose_print_test ( ):
#*****************************************************************************80
#
## r8ge_transpose_print_test() tests r8ge_transpose_print().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 October 2014
#
# Author:
#
# John Burkardt
#
import numpy as np
print ( '' )
print ( 'r8ge_transpose_print_test():' )
print ( ' r8ge_transpose_print prints the transpose of an R8GE matrix.' )
m = 4
n = 3
v = np.array ( [ \
[ 11.0, 12.0, 13.0 ],
[ 21.0, 22.0, 23.0 ],
[ 31.0, 32.0, 33.0 ],
[ 41.0, 42.0, 43.0 ] ], dtype = np.float64 )
r8ge_transpose_print ( m, n, v, ' Here is an R8GE matrix, transposed:' )
return
def r8ge_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## r8ge_transpose_print_some() prints a portion of an R8GE matrix, transposed.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 13 November 2014
#
# 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 i2lo in range ( max ( ilo, 0 ), min ( ihi, m - 1 ), incx ):
i2hi = i2lo + incx - 1
i2hi = min ( i2hi, m - 1 )
i2hi = min ( i2hi, ihi )
print ( '' )
print ( ' Row: ' ),
for i in range ( i2lo, i2hi + 1 ):
print ( '%7d ' % ( i ) ),
print ( '' )
print ( ' Col' )
j2lo = max ( jlo, 0 )
j2hi = min ( jhi, n - 1 )
for j in range ( j2lo, j2hi + 1 ):
print ( '%7d :' % ( j ) ),
for i in range ( i2lo, i2hi + 1 ):
print ( '%12g ' % ( a[i,j] ) ),
print ( '' )
return
def r8ge_transpose_print_some_test ( ):
#*****************************************************************************80
#
## r8ge_transpose_print_some_test() tests r8ge_transpose_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 October 2014
#
# Author:
#
# John Burkardt
#
import numpy as np
print ( '' )
print ( 'r8ge_transpose_print_some_test():' )
print ( ' r8ge_transpose_print_some() prints some of an R8GE matrix, transposed.' )
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_transpose_print_some ( m, n, v, 0, 3, 2, 5, \
' R8GE matrix, rows 0:2, cols 3:5:' )
return
def r8ge_trf ( m, n, a ):
#*****************************************************************************80
#
## r8ge_trf() performs a LAPACK-style PLU factorization of 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.
#
# r8ge_trf is a standalone version of the LAPACK routine R8GETRF.
#
# The factorization uses partial pivoting with row interchanges,
# and has the form
# A = P * L * U
# where P is a permutation matrix, L is lower triangular with unit
# diagonal elements (lower trapezoidal if N < M), and U is upper
# triangular (upper trapezoidal if M < N).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt.
#
# Reference:
#
# Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum,
# Hammarling, McKenney, Ostrouchov, Sorensen,
# LAPACK User's Guide,
# Second Edition,
# SIAM, 1995.
#
# Input:
#
# integer M, the number of rows of the matrix A. 0 <= M.
#
# integer N, the number of columns of the matrix A. 0 <= N.
#
# real A(M,N), the M by N matrix to be factored.
#
# Output:
#
# real A_LU(M,N), the factors L and U from the factorization
# A = P*L*U the unit diagonal elements of L are not stored.
#
# integer PIVOT(min(M,N)), the pivot indices
# for 1 <= I <= min(M,N), row i of the matrix was interchanged with
# row PIVOT(I).
#
# integer INFO.
# = 0: successful computation.
# < 0: if INFO = -K, the K-th argument had an illegal value
# > 0: if INFO = K, U(K,K) is exactly zero. The factorization
# has been completed, but the factor U is exactly
# singular, and division by zero will occur if it is used
# to solve a system of equations.
#
import numpy as np
a_lu = a.copy ( )
pivot = np.zeros ( n, dtype = np.int32 )
#
# Test the input parameters.
#
info = 0
if ( m < 0 ):
info = -1
return a_lu, pivot, info
if ( n < 0 ):
info = -2
return a_lu, pivot, info
if ( m == 0 or n == 0 ):
return a_lu, pivot, info
for j in range ( 0, min ( m, n ) ):
#
# Find the pivot.
#
t = abs ( a_lu[j,j] )
jp = j
for i in range ( j + 1, m ):
if ( t < abs ( a_lu[i,j] ) ):
t = abs ( a_lu[i,j] )
jp = i
pivot[j] = jp
#
# Apply the interchange to columns 1:N.
# Compute elements J+1:M of the J-th column.
#
if ( a_lu[jp,j] != 0.0 ):
if ( jp != j ):
for jj in range ( 0, n ):
t = a_lu[j,jj]
a_lu[j,jj] = a_lu[jp,jj]
a_lu[jp,jj] = t
if ( j < m - 1 ):
for i in range ( j + 1, m ):
a_lu[i,j] = a_lu[i,j] / a_lu[j,j]
elif ( info == 0 ):
info = j
#
# Update the trailing submatrix.
#
if ( j < min ( m, n ) - 1 ):
for ii in range ( j + 1, m ):
for k in range ( j + 1, n ):
a_lu[ii,k] = a_lu[ii,k] - a_lu[ii,j] * a_lu[j,k]
return a_lu, pivot, info
def r8ge_trf_test ( ):
#*****************************************************************************80
#
## r8ge_trf_test() tests r8ge_trf().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 5
m = n
nrhs = 1
print ( '' )
print ( 'r8ge_trf_test():' )
print ( ' r8ge_trf() computes the LU factors of an R8GE matrix,' )
print ( ' so that r8ge_trs() can solve the factored system.' )
print ( '' )
print ( ' Number of matrix rows M = %d' % ( m ) )
print ( ' Number of matrix columns N = %d' % ( n ) )
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
if ( i == j ):
a[i,j] = 2.0
elif ( i == j - 1 ):
a[i,j] = - 1.0
elif ( i == j + 1 ):
a[i,j] = - 1.0
else:
a[i,j] = 0.0
a_lu, pivot, info = r8ge_trf ( m, n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trf_test(): Fatal error!' )
print ( ' r8ge_trf declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
return
b = np.zeros ( [ n, nrhs ] )
b[n-1,0] = float ( n + 1 )
x, info = r8ge_trs ( n, nrhs, 'N', a_lu, pivot, b )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trf_test(): Fatal error!' )
print ( ' r8ge_trs() returned an error condition!' )
print ( ' The value of INFO is %d' % ( info ) )
return
r8vec_print ( n, x, ' Solution:' )
b = np.zeros ( [ n, nrhs ] )
b[n-1,0] = float ( n + 1 )
x, info = r8ge_trs ( n, nrhs, 'T', a_lu, pivot, b )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trf_test(): Fatal error!' )
print ( ' r8ge_trs returned an error condition!' )
print ( ' The value of INFO is %d' % ( info ) )
return
r8vec_print ( n, x, ' Solution to transposed system:' )
return
def r8ge_trs ( n, nrhs, trans, a_lu, pivot, b ):
#*****************************************************************************80
#
## r8ge_trs() solves a system of linear equations factored by r8ge_trf.
#
# Discussion:
#
# Note that in MATLAB we will have peculiar and maddening problems
# if our input data B is actually a vector in fact, if B is a vector,
# we must do something like call r8vec_to_r8ge in order to make it
# look like a 2D array to MATLAB.
#
# 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.
#
# r8ge_trs is a standalone version of the LAPACK routine R8GETRS.
#
# r8ge_trs solves a system of linear equations
# A * x = b or A' * X = B
# with a general N by N matrix A using the PLU factorization computed
# by r8ge_trf.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt.
#
# Reference:
#
# Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum,
# Hammarling, McKenney, Ostrouchov, Sorensen,
# LAPACK User's Guide,
# Second Edition,
# SIAM, 1995.
#
# Input:
#
# integer N, the order of the matrix A. 0 <= N.
#
# integer NRHS, the number of right hand sides. 0 <= NRHS.
#
# character TRANS, specifies the form of the system of equations:
# 'N': A * x = b (No transpose)
# 'T': A'* X = B (Transpose)
# 'C': A'* X = B (Conjugate transpose = Transpose)
#
# real A_LU(N,N), the LU factors from r8ge_trf.
#
# integer PIVOT(N), the pivot indices from r8ge_trf.
#
# real B(N,NRHS), the right hand side matrix B.
#
# Output:
#
# real X(N,NRHS), the solution matrix X.
#
# integer INFO
# = 0: successful computation.
# < 0: if INFO = -I, the I-th argument had an illegal value.
#
x = b.copy ( )
info = 0
if ( trans != 'n' and trans != 'N' and trans != 't' and trans != 'T' and \
trans != 'c' and trans != 'C' ):
info = -1
return x, info
elif ( n < 0 ):
info = -2
return x, info
elif ( nrhs < 0 ):
info = -3
return x, info
if ( n == 0 or nrhs == 0 ):
return x, info
if ( trans == 'n' or trans == 'N' ):
#
# Apply row interchanges to the right hand sides.
#
for i in range ( 0, n ):
if ( pivot[i] != i ):
for k in range ( 0, nrhs ):
t = x[i,k]
x[i,k] = x[pivot[i],k]
x[pivot[i],k] = t
#
# Solve L * x = b, overwriting b with x.
#
for k in range ( 0, nrhs ):
for j in range ( 0, n - 1 ):
for i in range ( j + 1, n ):
x[i,k] = x[i,k] - a_lu[i,j] * x[j,k]
#
# Solve U * x = b, overwriting b with x.
#
for k in range ( 0, nrhs ):
for j in range ( n - 1, -1, -1 ):
x[j,k] = x[j,k] / a_lu[j,j]
else:
#
# Solve U' * x = b, overwriting b with x.
#
for k in range ( 0, nrhs ):
for j in range ( 0, n ):
x[j,k] = x[j,k] / a_lu[j,j]
for i in range ( j + 1, n ):
x[i,k] = x[i,k] - a_lu[j,i] * x[j,k]
#
# Solve L' * x = b, overwriting b with x.
#
for k in range ( 0, nrhs ):
for j in range ( n - 1, 0, -1 ):
for i in range ( 0, j ):
x[i,k] = x[i,k] - a_lu[j,i] * x[j,k]
#
# Apply row interchanges to the solution vectors.
#
for i in range ( n - 1, -1, -1 ):
if ( pivot[i] != i ):
for k in range ( 0, nrhs ):
t = x[i,k]
x[i,k] = x[pivot[i],k]
x[pivot[i],k] = t
return x, info
def r8ge_trs_test ( ):
#*****************************************************************************80
#
## r8ge_trs_test() tests r8ge_trs().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 5
m = n
nrhs = 1
print ( '' )
print ( 'r8ge_trs_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_trs solves a linear system' )
print ( ' that has been factored by r8ge_trf.' )
print ( '' )
print ( ' Number of matrix rows M = %d' % ( m ) )
print ( ' Number of matrix columns N = %d' % ( n ) )
a = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
for j in range ( 0, n ):
if ( i == j ):
a[i,j] = 2.0
elif ( i == j - 1 ):
a[i,j] = - 1.0
elif ( i == j + 1 ):
a[i,j] = - 1.0
else:
a[i,j] = 0.0
a_lu, pivot, info = r8ge_trf ( m, n, a )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trs_test(): Fatal error!' )
print ( ' r8ge_trf declares the matrix is singular!' )
print ( ' The value of INFO is %d' % ( info ) )
return
b = np.zeros ( [ n, nrhs ] )
b[n-1,0] = float ( n + 1 )
x, info = r8ge_trs ( n, nrhs, 'N', a_lu, pivot, b )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trs_test(): Fatal error!' )
print ( ' r8ge_trs returned an error condition!' )
print ( ' The value of INFO is %d' % ( info ) )
return
r8vec_print ( n, x, ' Solution:' )
b = np.zeros ( [ n, nrhs ] )
b[n-1,0] = float ( n + 1 )
x, info = r8ge_trs ( n, nrhs, 'T', a_lu, pivot, b )
if ( info != 0 ):
print ( '' )
print ( 'r8ge_trs_test(): Fatal error!' )
print ( ' r8ge_trs returned an error condition!' )
print ( ' The value of INFO is %d' % ( info ) )
return
r8vec_print ( n, x, ' Solution to transposed system:' )
return
def r8ge_zeros ( m, n ):
#*****************************************************************************80
#
## r8ge_zeros() zeroes an 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:
#
# 01 August 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the order of the matrix.
# N must be positive.
#
# Output:
#
# real A(M,N), the zeroed out matrix.
#
import numpy as np
a = np.zeros ( [ m, n ] )
return a
def r8ge_zeros_test ( ):
#*****************************************************************************80
#
## r8ge_zeros_test() tests r8ge_zeros().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 01 August 2015
#
# Author:
#
# John Burkardt
#
import platform
m = 5
n = 4
print ( '' )
print ( 'r8ge_zeros_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8ge_zeros zeros out space for a general matrix.' )
print ( '' )
print ( ' Matrix order M, N = %d, %d' % ( m, n ) )
a = r8ge_zeros ( m, n )
r8ge_print ( m, n, a, ' Matrix A:' )
return
def r8vec_house_column ( n, a_vec, k ):
#*****************************************************************************80
#
## r8vec_house_column() defines a Householder premultiplier that "packs" a column.
#
# Discussion:
#
# The routine returns a vector V that defines a Householder
# premultiplier matrix H(V) that zeros out the subdiagonal entries of
# column K of the matrix A.
#
# H(V) = I - 2 * v * v'
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 19 February 2015
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the order of the matrix A.
#
# real A_VEC(N), a row or column of the matrix A.
#
# integer K, the "special" entry in A_VEC.
# The Householder matrix will zero out the entries after this.
#
# Output:
#
# real V(N), a vector of unit L2 norm which defines an
# orthogonal Householder premultiplier matrix H with the property
# that the K-th column of H*A is zero below the diagonal.
#
import numpy as np
v = np.zeros ( n )
if ( k < 0 or n - 1 <= k ):
return v
s = 0.0
for i in range ( k, n ):
s = s + a_vec[i] ** 2
s = np.sqrt ( s )
if ( s == 0.0 ):
return v
if ( a_vec[k] < 0.0 ):
v[k] = a_vec[k] - abs ( s )
else:
v[k] = a_vec[k] + abs ( s )
for i in range ( k + 1, n ):
v[i] = a_vec[i]
s = 0.0
for i in range ( k, n ):
s = s + v[i] ** 2
s = np.sqrt ( s )
for i in range ( k, n ):
v[i] = v[i] / s
return v
def r8vec_house_column_test ( ):
#*****************************************************************************80
#
## r8vec_house_column_test() tests r8vec_house_column().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 14 February 2015
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
print ( '' )
print ( 'r8vec_house_column_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8vec_house_column returns the compact form of' )
print ( ' a Householder matrix that "packs" a column' )
print ( ' of a matrix.' )
#
# Get a random matrix.
#
n = 4
r8_lo = 0.0
r8_hi = 5.0
a = r8ge_random_ab ( n, n, r8_lo, r8_hi )
r8ge_print ( n, n, a, ' Matrix A:' )
a_col = np.zeros ( n )
for k in range ( 0, n - 1 ):
print ( '' )
print ( ' Working on column K = %d' % ( k ) )
for i in range ( 0, n ):
a_col[i] = a[i,k]
v = r8vec_house_column ( n, a_col, k )
h = r8ge_house_form ( n, v )
r8ge_print ( n, n, h, ' Householder matrix H:' )
ha = r8ge_mm ( n, n, n, h, a )
r8ge_print ( n, n, ha, ' Product H*A:' )
#
# If we set A := HA, then we can successively convert A to upper
# triangular form.
#
a = ha
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] ) )
def r8vec_print_some ( n, a, max_print, title ):
#*****************************************************************************80
#
## r8vec_print_some() prints "some" of an R8VEC.
#
# Discussion:
#
# The user specifies MAX_print, the maximum number of lines to print.
#
# If N, the size of the vector, is no more than MAX_print, then
# the entire vector is printed, one entry 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:
#
# 03 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the number of entries of the vector.
#
# real A(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' % ( i, a[i] ) )
elif ( 3 <= max_print ):
for i in range ( 0, max_print - 2 ):
print ( ' %6d %14g' % ( i, a[i] ) )
print ( ' ...... ..............' )
i = n - 1
print ( ' %6d %14g' % ( i, a[i] ) )
else:
for i in range ( 0, max_print - 1 ):
print ( ' %6d %14g' % ( i, a[i] ) )
i = max_print - 1
print ( ' %6d %14g ...more entries...' % ( i, a[i] ) )
return
def r8vec_print_some_test ( ):
#*****************************************************************************80
#
## r8vec_print_some_test() tests r8vec_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 03 February 2016
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8vec_print_some_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8vec_print_some prints some of an R8VEC.' )
n = 100
a = r8vec_indicator1 ( n )
max_print = 10
r8vec_print_some ( n, a, max_print, ' No more than 10 lines of this vector:' )
return
def r8vec_to_r8ge ( m, n, x ):
#*****************************************************************************80
#
## r8vec_to_r8ge() copies an R8VEC into a R8GE matrix.
#
# Discussion:
#
# In C++ and FORTRAN, this routine is not really needed. In MATLAB,
# a data item carries its dimensionality implicitly, and so cannot be
# regarded sometimes as a vector and sometimes as an array.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, N, the number of rows and columns in the array.
#
# real X(M*N), the vector to be copied into the array.
#
# Output:
#
# real A(M,N), the array.
#
import numpy as np
a = np.zeros ( [ m, n ] )
k = 0
for j in range ( 0, n ):
for i in range ( 0, m ):
a[i,j] = x[k]
k = k + 1
return a
def r8vec_to_r8ge_test ( ):
#*****************************************************************************80
#
## r8vec_to_r8ge_test() tests r8vec_to_r8ge().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
m = 4
n = 3
print ( '' )
print ( 'r8vec_to_r8ge_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' r8vec_to_r8ge converts an R8VEC vector to an R8GE matrix.' )
a_r8vec = r8vec_indicator1 ( m * n )
r8vec_print ( m * n, a_r8vec, ' The R8VEC vector:' )
a_r8ge = r8vec_to_r8ge ( m, n, a_r8vec )
r8ge_print ( m, n, a_r8ge, ' Corresponding R8GE matrix:' )
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 r8vec2_print_some_test ( ):
#*****************************************************************************80
#
## r8vec2_print_some_test() tests r8vec2_print_some().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 January 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
n = 100
a = np.zeros ( n )
b = np.zeros ( n )
for i in range ( 0, n ):
x = float ( i + 1 )
a[i] = x * x
b[i] = np.sqrt ( x )
print ( '' )
print ( 'r8vec2_print_some_test():' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' r8vec2_print_some prints some of a pair of R8VEC\'s.' )
r8vec2_print_some ( n, a, b, 10, ' Square and square root:' )
return
def timestamp ( ):
#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 April 2013
#
# Author:
#
# John Burkardt
#
import time
t = time.time ( )
print ( time.ctime ( t ) )
return None
def r8ge_test ( ):
#*****************************************************************************80
#
## r8ge_test() tests r8ge().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 February 2016
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8ge_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' Test r8ge().' )
i4_log_10_test ( )
r8_sign_test ( )
r8ge_cg_test ( )
r8ge_co_test ( )
r8ge_det_test ( )
r8ge_dif2_test ( )
r8ge_dilu_test ( )
r8ge_fa_test01 ( )
r8ge_fa_test02 ( )
r8ge_fs_test ( )
r8ge_fss_test ( )
r8ge_hilbert_test ( )
r8ge_hilbert_inverse_test ( )
r8ge_house_axh_test ( )
r8ge_house_form_test ( )
r8ge_identity_test ( )
r8ge_ilu_test ( )
r8ge_indicator_test ( )
r8ge_inverse_test ( )
r8ge_ml_test ( )
r8ge_mm_test ( )
r8ge_mtm_test ( )
r8ge_mtv_test ( )
r8ge_mu_test ( )
r8ge_mv_test ( )
r8ge_orth_random_test ( )
r8ge_spd_random_test ( )
r8ge_plu_test ( )
r8ge_poly_test ( )
r8ge_print_test ( )
r8ge_print_some_test ( )
r8ge_random_test ( )
r8ge_random_ab_test ( )
r8ge_res_test ( )
r8ge_sl_test ( )
r8ge_sl_it_test ( )
# r8ge_to_r8lt_test ( )
r8ge_to_r8po_test ( )
r8ge_to_r8pp_test ( )
# r8ge_to_r8ut_test ( )
r8ge_to_r8vec_test ( )
# r8ge_to_r8vm_test ( )
r8ge_transpose_test ( )
r8ge_transpose_print_test ( )
r8ge_transpose_print_some_test ( )
r8ge_trf_test ( )
r8ge_trs_test ( )
r8ge_zeros_test ( )
r8vec_house_column_test ( )
r8vec_indicator1_test ( )
r8vec_to_r8ge_test ( )
r8vec2_print_some_test ( )
#
# Terminate.
#
print ( '' )
print ( 'r8ge_test():' )
print ( ' Normal end of execution.' )
return
if ( __name__ == '__main__' ):
timestamp ( )
r8ge_test ( )
timestamp ( )