#! /usr/bin/env python3
#
def r8cbb_test ( ):
#*****************************************************************************80
#
## r8cbb_test() tests r8cbb().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 August 2022
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'r8cbb_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' Test r8cbb().' )
r8cbb_add_test ( )
r8cbb_dif2_test ( )
r8cbb_fa_test ( )
r8cbb_get_test ( )
r8cbb_indicator_test ( )
r8cbb_mtv_test ( )
r8cbb_mv_test ( )
r8cbb_print_test ( )
r8cbb_print_some_test ( )
r8cbb_random_test ( )
r8cbb_set_test ( )
r8cbb_sl_test ( )
r8cbb_to_r8ge_test ( )
r8cbb_zeros_test ( )
#
# Terminate.
#
print ( '' )
print ( 'r8cbb_test():' )
print ( ' Normal end of execution.' )
return
def i4_log_10 ( i ):
#*****************************************************************************80
#
## i4_log_10() returns the integer part of the logarithm base 10 of ABS(X).
#
# Example:
#
# I VALUE
# ----- --------
# 0 0
# 1 0
# 2 0
# 9 0
# 10 1
# 11 1
# 99 1
# 100 2
# 101 2
# 999 2
# 1000 3
# 1001 3
# 9999 3
# 10000 4
#
# Discussion:
#
# i4_log_10 ( I ) + 1 is the number of decimal digits in I.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 08 May 2013
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer I, the number whose logarithm base 10 is desired.
#
# Output:
#
# integer VALUE, the integer part of the logarithm base 10 of
# the absolute value of X.
#
import numpy as np
i = np.floor ( i )
if ( i == 0 ):
value = 0
else:
value = 0
ten_pow = 10
i_abs = abs ( i )
while ( ten_pow <= i_abs ):
value = value + 1
ten_pow = ten_pow * 10
return value
def r8cb_np_fa ( n, ml, mu, a ):
#*****************************************************************************80
#
## r8cb_np_fa() factors a R8CB matrix by Gaussian elimination.
#
# Discussion:
#
# The R8CB storage format is appropriate for a compact banded matrix.
# It is assumed that the matrix has lower and upper bandwidths ML and MU,
# respectively. The matrix is stored in a way similar to that used
# by LINPACK and LAPACK for a general banded matrix, except that in
# this mode, no extra rows are set aside for possible fillin during pivoting.
# Thus, this storage mode is suitable if you do not intend to factor
# the matrix, or if you can guarantee that the matrix can be factored
# without pivoting.
#
# R8CB_NP_FA is a version of the LINPACK routine R8GBFA, modifed to use
# no pivoting, and to be applied to the R8CB compressed band matrix storage
# format. It will fail if the matrix is singular, or if any zero
# pivot is encountered.
#
# If R8CB_NP_FA successfully factors the matrix, R8CB_NP_SL may be called
# to solve linear systems involving the matrix.
#
# The matrix is stored in a compact version of LINPACK general
# band storage, which does not include the fill-in entires.
# The following program segment will store the entries of a banded
# matrix in the compact format used by this routine:
#
# m = mu+1
# do j = 1, n
# i1 = max ( 1, j-mu )
# i2 = min ( n, j+ml )
# do i = i1, i2
# k = i-j+m
# a(k,j) = afull(i,j)
# end do
# end do
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 26 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N-1.
#
# Input, real A(ML+MU+1,N), the compact band matrix.
#
# Output, real A_LU(ML+MU+1,N), the LU factors of the matrix.
#
# Output, integer INFO, singularity flag.
# 0, no singularity detected.
# nonzero, the factorization failed on the INFO-th step.
#
a_lu = a.copy ( )
#
# The value of M is MU + 1 rather than ML + MU + 1.
#
m = mu + 1
info = 0
ju = 0
for k in range ( 1, n ):
#
# If our pivot entry A(MU+1,K) is zero, then we must give up.
#
if ( a_lu[m,k-1] == 0.0 ):
info = k
print ( '' )
print ( 'R8CB_NP_FA(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'R8CB_NP_FA(): Fatal error!' )
#
# LM counts the number of nonzero elements that lie below the current
# diagonal entry, A(K,K).
#
# Multiply the LM entries below the diagonal by -1/A(K,K), turning
# them into the appropriate "multiplier" terms in the L matrix.
#
lm = min ( ml, n - k )
for i in range ( m + 1, m + lm + 1 ):
a_lu[i-1,k-1] = - a_lu[i-1,k-1] / a_lu[m-1,k-1]
#
# MM points to the row in which the next entry of the K-th row is, A(K,J).
# We then add L(I,K)*A(K,J) to A(I,J) for rows I = K+1 to K+LM.
#
ju = max ( ju, mu + k )
ju = min ( ju, n )
mm = m
for j in range ( k + 1, ju + 1 ):
mm = mm - 1
for i in range ( 0, lm ):
a_lu[mm+i,j-1] = a_lu[mm+i,j-1] + a_lu[mm-1,j-1] * a_lu[m+i,k-1]
if ( a_lu[m-1,n-1] == 0.0 ):
info = n
print ( '' )
print ( 'R8CB_NP_FA(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'R8CB_NP_FA(): Fatal error!' )
return a_lu, info
def r8cb_np_sl ( n, ml, mu, a_lu, b, job ):
#*****************************************************************************80
#
## R8CB_NP_SL solves a R8CB system factored by R8CB_NP_FA.
#
# Discussion:
#
# The R8CB storage format is appropriate for a compact banded matrix.
# It is assumed that the matrix has lower and upper bandwidths ML and MU,
# respectively. The matrix is stored in a way similar to that used
# by LINPACK and LAPACK for a general banded matrix, except that in
# this mode, no extra rows are set aside for possible fillin during pivoting.
# Thus, this storage mode is suitable if you do not intend to factor
# the matrix, or if you can guarantee that the matrix can be factored
# without pivoting.
#
# R8CB_NP_SL can also solve the related system A' * x = b.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 26 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N-1.
#
# Input, real A_LU(ML+MU+1,N), the compact band matrix, factored by R8CB_NP_FA.
#
# Input, real B(N), the right hand side of the linear system.
#
# Input, integer JOB.
# If JOB is zero, the routine will solve A * x = b.
# If JOB is nonzero, the routine will solve A' * x = b.
#
# Output, real X(N), the solution of the linear system.
#
x = b.copy ( )
#
# The value of M is ML + 1, rather than MU + ML + 1.
#
m = mu + 1
#
# Solve A * x = b.
#
if ( job == 0 ):
#
# Solve PL * Y = B.
#
if ( 0 < ml ):
for k in range ( 1, n ):
lm = min ( ml, n - k )
for i in range ( 1, lm + 1 ):
x[k+i-1] = x[k+i-1] + x[k-1] * a_lu[m+i-1,k-1]
#
# Solve U * X = Y.
#
for k in range ( n, 0, -1 ):
x[k-1] = x[k-1] / a_lu[m-1,k-1]
lm = min ( k, m ) - 1
la = m - lm
lb = k - lm
for i in range ( 0, lm ):
x[lb+i-1] = x[lb+i-1] - x[k-1] * a_lu[la+i-1,k-1]
#
# Solve A' * X = B.
#
else:
#
# Solve U' * Y = B.
#
for k in range ( 1, n + 1 ):
lm = min ( k, m ) - 1
la = m - lm
lb = k - lm
for i in range ( 0, lm ):
x[k-1] = x[k-1] - a_lu[la+i-1,k-1] * x[lb+i-1]
x[k-1] = x[k-1] / a_lu[m-1,k-1]
#
# Solve ( PL )' * X = Y.
#
if ( 1 <= ml ):
for k in range ( n - 1, 0, -1 ):
lm = min ( ml, n - k )
for i in range ( 1, lm + 1 ):
x[k-1] = x[k-1] + a_lu[m+i-1,k-1] * x[k+i-1]
return x
def r8cb_to_r8vec ( m, n, ml, mu, a ):
#*****************************************************************************80
#
## r8cb_to_r8vec() copies a R8CB 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:
#
# 26 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the number of rows and columns in the array.
#
# Input, integer ML, MU, the lower and upper bandwidths.
#
# Input, real A(ML+MU+1,N), the array to be copied.
#
# Output, real X((ML+MU+1)*N), the vector.
#
import numpy as np
x = np.zeros ( ( ml + mu + 1 ) * n )
for j in range ( 0, n ):
ilo = max ( mu - j, 0 )
ihi = mu + min ( ml, m - j - 1 )
for i in range ( ilo, ihi + 1 ):
x[i+j*(ml+mu+1)] = a[i,j]
return x
def r8cbb_add ( n1, n2, ml, mu, a, i, j, value ):
#*****************************************************************************80
#
## r8cbb_add() adds a value to an entry of a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= n1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
# Input, integer I, J, the indices of the entry to be incremented.
#
# Input, real VALUE, the value to be added to the (I,J) entry.
#
# Output, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the modified R8CBB matrix.
#
import numpy as np
# a = np.zeros ( ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 )
if ( value == 0.0 ):
return a
#
# Check for I or J out of bounds.
#
if ( i < 0 or n1 + n2 <= i ):
print ( '' )
print ( 'R8CBB_ADD(): Fatal error!' )
print ( ' Illegal input value of row index I = ', i )
raise Exception ( 'R8CBB_ADD(): Fatal error!' )
if ( j < 0 or n1 + n2 <= j ):
print ( '' )
print ( 'R8CBB_ADD(): Fatal error!' )
print ( ' Illegal input value of column index J = ', j )
raise Exception ( 'R8CBB_ADD(): Fatal error!' )
#
# The A1 block of the matrix.
#
# Check for out of band problems.
#
if ( i < n1 and j < n1 ):
if ( mu < j - i or ml < i - j ):
print ( '' )
print ( 'R8CBB_ADD(): Fatal error!' )
print ( ' Unable to add to entry (', i, ',', j, ').' )
raise Exception ( 'R8CBB_ADD(): Fatal error!' )
ij = ( i - j + mu ) + j * ( ml + mu + 1 )
#
# The A2 block of the matrix:
#
elif ( i < n1 and n1 <= j ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i
#
# The A3 and A4 blocks of the matrix.
#
elif ( n1 <= i ):
ij = ( ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 )
a[ij] = a[ij] + value
return a
def r8cbb_add_test ( ):
#*****************************************************************************80
#
## R8CBB_ADD_TEST tests R8CBB_ADD.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 3
n2 = 2
n = n1 + n2
ml = 1
mu = 0
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_ADD_TEST' )
print ( ' R8CBB_ADD adds a value to elements of an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Initialize matrix to indicator matrix.
#
a = r8cbb_indicator ( n1, n2, ml, mu )
#
# Print initial matrix.
#
r8cbb_print ( n1, n2, ml, mu, a, ' Matrix before additions:' )
#
# Add 100 to band diagonal.
#
for i in range ( 0, n1 ):
j = i
value = 100.0
a = r8cbb_add ( n1, n2, ml, mu, a, i, j, value )
#
# Add 200 to right border.
#
for i in range ( 0, n1 ):
for j in range ( n1, n1 + n2 ):
value = 200.0
a = r8cbb_add ( n1, n2, ml, mu, a, i, j, value )
#
# Add 400 to offdiagonals in lower right dense matrix.
#
for i in range ( n1, n1 + n2 ):
for j in range ( n1, n1 + n2 ):
if ( i != j ):
value = 400.0
a = r8cbb_add ( n1, n2, ml, mu, a, i, j, value )
r8cbb_print ( n1, n2, ml, mu, a, ' The matrix after additions:' )
return
def r8cbb_dif2 ( n1, n2, ml, mu ):
#*****************************************************************************80
#
## R8CBB_DIF2 sets up an R8CBB second difference matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense
# blocks. N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# 1 <= ML, 1 <= MU.
#
# Output, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
#
import numpy as np
a = np.zeros ( ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 )
if ( ml < 1 or mu < 1 ):
print ( '' )
print ( 'R8CBB_DIF2(): Fatal error!' )
print ( ' 1 <= ML and 1 <= MU required.' )
raise Exception ( 'R8CBB_DIF2(): Fatal error!' )
for i in range ( 1, n1 + n2 ):
j = i - 1
value = - 1.0
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
for i in range ( 0, n1 + n2 ):
j = i
value = 2.0
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
for i in range ( 0, n1 + n2 - 1 ):
j = i + 1
value = - 1.0
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
return a
def r8cbb_dif2_test ( ):
#*****************************************************************************80
#
## R8CBB_DIF2_TEST tests R8CBB_DIF2.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_DIF2_TEST' )
print ( ' R8CBB_DIF2 sets up an R8CBB second difference matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
a = r8cbb_dif2 ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB second difference matrix:' )
return
def r8cbb_fa ( n1, n2, ml, mu, a ):
#*****************************************************************************80
#
## R8CBB_FA factors a R8CBB matrix.
#
# Discussion:
#
# Note that in C++ and FORTRAN, we can look at A as an abstract
# vector, but then look at parts of A as storing a two dimensional
# array. MATLAB assigns an inherent dimensionality to a data object,
# and gets very unhappy when you try to manipulate the data yourself.
# This means that the MATLAB implementation of this routine requires
# the use of temporary 2D arrays.
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
#
# Once the matrix has been factored by R8CBB_FA, R8CBB_SL may be called
# to solve linear systems involving the matrix.
#
# R8CBB_FA uses special non-pivoting versions of LINPACK routines to
# carry out the factorization. The special version of the banded
# LINPACK solver also results in a space saving, since no entries
# need be set aside for fill in due to pivoting.
#
# The linear system must be border banded, of the form:
#
# ( A1 A2 ) (X1) = (B1)
# ( A3 A4 ) (X2) (B2)
#
# where A1 is a (usually big) banded square matrix, A2 and A3 are
# column and row strips which may be nonzero, and A4 is a dense
# square matrix.
#
# The algorithm rewrites the system as:
#
# X1 + inverse(A1) A2 X2 = inverse(A1) B1
#
# A3 X1 + A4 X2 = B2
#
# and then rewrites the second equation as
#
# ( A4 - A3 inverse(A1) A2 ) X2 = B2 - A3 inverse(A1) B1
#
# The algorithm will certainly fail if the matrix A1 is singular,
# or requires pivoting. The algorithm will also fail if the A4 matrix,
# as modified during the process, is singular, or requires pivoting.
# All these possibilities are in addition to the failure that will
# if the total matrix A is singular.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 30 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A( (ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the compact
# border-banded coefficient matrix.
#
# Output, real A_LU( (ML+MU+1)*N1 + 2*N1*N2 + N2*N2).
# information describing a partial factorization
# of the original coefficient matrix.
#
# Output, integer INFO, singularity flag.
# 0, no singularity detected.
# nonzero, the factorization failed on the INFO-th step.
#
import numpy as np
nband = ( ml + mu + 1 ) * n1
#
# A_LU should be as big as A.
#
a_lu = a.copy ( )
#
# A1 is the compact band matrix portion of A.
# A1_LU is the LU factorization of A1.
# A_LU can store A1_LU.
#
if ( 0 < n1 ):
a1 = r8vec_to_r8cb ( n1, n1, ml, mu, a )
a1_lu, info = r8cb_np_fa ( n1, ml, mu, a1 )
if ( info != 0 ):
print ( '' )
print ( 'R8CBB_FA(): Fatal error!' )
print ( ' R8CB_NP_FA returned INFO = ', info )
print ( ' Factoring failed for column INFO.' )
print ( ' The band matrix A1 is singular.' )
print ( ' This algorithm cannot continue!' )
raise Exception ( 'R8CBB_FA(): Fatal error!' )
a_lu[0:nband] = r8cb_to_r8vec ( n1, n1, ml, mu, a1_lu )
if ( 0 < n1 and 0 < n2 ):
#
# A2 is the N1 x N2 strip to the right of A1.
# Let B2 be column J of A2.
# Let X2 solve A1 * X2 = - B2.
# Store each X2 into A_LU.
#
b2 = np.zeros ( n1 )
job = 0
for j in range ( 0, n2 ):
b2[0:n1] = - a[nband+j*n1:nband+j*n1+n1]
x2 = r8cb_np_sl ( n1, ml, mu, a1_lu, b2, job )
a_lu[nband+j*n1:nband+j*n1+n1] = x2[0:n1]
#
# Set A4 := A4 + A3*A2
#
for i in range ( 0, n2 ):
for j in range ( 0, n1 ):
ij = nband + n1 * n2 + j * n2 + i
for k in range ( 0, n2 ):
ik = nband + 2 * n1 * n2 + k * n2 + i
jk = nband + k * n1 + j
a_lu[ik] = a_lu[ik] + a_lu[ij] * a_lu[jk]
#
# Factor A4.
#
if ( 0 < n2 ):
a4 = r8vec_to_r8ge ( n2, n2, a_lu[nband+2*n1*n2:nband+2*n1*n2+n2*n2] )
a4_lu, info = r8ge_np_fa ( n2, a4 )
if ( info != 0 ):
print ( '' )
print ( 'R8CBB_FA(): Fatal error!' )
print ( ' R8GE_NP_FA returned INFO = ',info )
print ( ' This indicates singularity in column INFO' )
info = n1 + info
print ( ' of the A4 submatrix, which is column ', info )
print ( ' of the full matrix.' )
print ( '' )
print ( ' It is possible that the full matrix is ' )
print ( ' nonsingular, but the algorithm R8CBB_FA may' )
print ( ' not be used for this matrix.' )
raise Exception ( 'R8CBB_FA(): Fatal error!' )
a_lu[nband+2*n1*n2:nband+2*n1*n2+n2*n2] = r8ge_to_r8vec ( n2, n2, a4_lu )
return a_lu, info
def r8cbb_fa_test ( ):
#*****************************************************************************80
#
## R8CBB_FA_TEST tests R8CBB_FA.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_FA_TEST' )
print ( ' R8CBB_FA factors an R8CBB matrix, with no pivoting.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_random ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix:' )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r8cbb_mv ( n1, n2, ml, mu, a, x )
#
# Factor the matrix
#
a_lu, info = r8cbb_fa ( n1, n2, ml, mu, a )
r8cbb_print ( n1, n2, ml, mu, a_lu, ' The factored matrix:' )
if ( info != 0 ):
print ( '' )
print ( 'R8CBB_FA_TEST(): Fatal error!' )
print ( ' R8CBB_FA claims the matrix is singular.' )
print ( ' The value of INFO is ', info )
raise Exception ( 'R8CBB_FA_TEST(): Fatal error!' )
#
# Solve the system.
#
r8vec_print ( n, b, ' The right hand side vector:' )
x = r8cbb_sl ( n1, n2, ml, mu, a_lu, b )
r8vec_print ( n, x, ' Solution:' )
return
def r8cbb_get ( n1, n2, ml, mu, a, i, j ):
#*****************************************************************************80
#
## R8CBB_GET returns the value of an entry of a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
# Input, integer I, J, the row and column of the entry to retrieve.
#
# Output, real VALUE, the value of the (I,J) entry.
#
#
# Check for I or J out of bounds.
#
if ( i < 0 or n1 + n2 <= i ):
print ( '' )
print ( 'R8CBB_GET(): Fatal error!' )
print ( ' Illegal input value of row index I = ', i )
raise Exception ( 'R8CBB_GET(): Fatal error!' )
if ( j < 0 or n1 + n2 <= j ):
print ( '' )
print ( 'R8CBB_GET(): Fatal error!' )
print ( ' Illegal input value of column index J = ', j )
raise Exception ( 'R8CBB_GET(): Fatal error!' )
#
# The A1 block of the matrix.
#
# Check for out of band problems.
#
if ( i <= n1 and j <= n1 ):
if ( mu < ( j - i ) or ml < ( i - j ) ):
value = 0.0
return value
ij = ( i - j + mu ) + j * ( ml + mu + 1 )
#
# The A2 block of the matrix:
#
elif ( i <= n1 and n1 < j ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i
#
# The A3 and A4 blocks of the matrix.
#
elif ( n1 < i ):
ij = ( ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 )
value = a[ij]
return value
def r8cbb_get_test ( ):
#*****************************************************************************80
#
## R8CBB_GET_TEST tests R8CBB_GET.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
from numpy.random import default_rng
rng = default_rng ( )
n1 = 3
n2 = 2
n = n1 + n2
ml = 1
mu = 0
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_GET_TEST' )
print ( ' R8CBB_GET gets a value of an element of an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set matrix to indicator matrix.
#
a = r8cbb_indicator ( n1, n2, ml, mu )
#
# Print matrix.
#
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix to be queried:' )
#
# Request random entries.
#
print ( '' )
for k in range ( 0, 10 ):
i = rng.integers ( 0, n1 + n2, size = 1, endpoint = False )
j = rng.integers ( 0, n1 + n2, size = 1, endpoint = False )
value = r8cbb_get ( n1, n2, ml, mu, a, i, j )
print ( ' A(', i, ',', j, ') = ', value )
return
def r8cbb_indicator ( n1, n2, ml, mu ):
#*****************************************************************************80
#
## R8CBB_INDICATOR sets up a R8CBB indicator matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the R8BB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative and no greater than N1-1.
#
# Output, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
#
import numpy as np
a = np.zeros ( ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 )
fac = 10 ** ( i4_log_10 ( n1 + n2 ) + 1 )
#
# Set the banded matrix A1.
#
for j in range ( 0, n1 ):
for row in range ( 0, ml + mu + 1 ):
i = row + j - mu
if ( 1 <= i and i <= n1 ):
a[row+j*(ml+mu+1)] = float ( fac * ( i + 1 ) + ( j + 1 ) )
#
# Set the N1 by N2 rectangular strip A2.
#
base = ( ml + mu + 1 ) * n1
for i in range ( 0, n1 ):
for j in range ( n1, n1 + n2 ):
a[base + i + (j-n1)*n1 ] = float ( fac * ( i + 1 ) + ( j + 1 ) )
#
# Set the N2 by N1 rectangular strip A3.
#
base = ( ml + mu + 1 ) * n1 + n1 * n2
for i in range ( n1, n1 + n2 ):
for j in range ( 0, n1 ):
a[base + i - n1 + j*n2 ] = float ( fac * ( i + 1 ) + ( j + 1 ) )
#
# Set the N2 by N2 square A4.
#
base = ( ml + mu + 1 ) * n1 + n1 * n2 + n2 * n1
for i in range ( n1, n1 + n2 ):
for j in range ( n1, n1 + n2 ):
a[base + i-n1 + (j-n1)*n2 ] = float ( fac * ( i + 1 ) + ( j + 1 ) )
return a
def r8cbb_indicator_test ( ):
#*****************************************************************************80
#
## R8CBB_INDICATOR_TEST tests R8CBB_INDICATOR.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_INDICATOR_TEST' )
print ( ' R8CBB_INDICATOR sets up an R8CBB indicator matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
a = r8cbb_indicator ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The compact border-banded matrix:' )
return
def r8cbb_mtv ( n1, n2, ml, mu, a, x ):
#*****************************************************************************80
#
## R8CBB_MTV multiplies a vector by a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= n1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
# Input, real X(N1+N2), the vector to multiply the matrix.
#
# Output, real B(N1+N2), the product X * A.
#
import numpy as np
#
# Set B to zero.
#
b = np.zeros ( n1 + n2 )
#
# Multiply by A1.
#
for j in range ( 0, n1 ):
ilo = max ( 0, j - mu )
ihi = min ( n1 - 1, j + ml )
ij = j * ( ml + mu + 1 ) - j + mu
for k in range ( ilo, ihi + 1 ):
b[j] = b[j] + x[k] * a[k+ij]
#
# Multiply by A2.
#
for j in range ( n1, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1
for k in range ( 0, n1 ):
b[j] = b[j] + x[k] * a[k+ij]
#
# Multiply by A3 and A4.
#
for j in range ( 0, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + n1 * n2 + j * n2 - n1
for k in range ( n1, n1 + n2 ):
b[j] = b[j] + x[k] * a[k+ij]
return b
def r8cbb_mtv_test ( ):
#*****************************************************************************80
#
## R8CBB_MTV_TEST tests R8CBB_MTV.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 6
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_MTV_TEST' )
print ( ' R8CBB_MTV computes b=A\'*x, where A is an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_indicator ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix A:' )
x = r8vec_indicator1 ( n )
r8vec_print ( n, x, ' The vector x:' )
b = r8cbb_mtv ( n1, n2, ml, mu, a, x )
r8vec_print ( n, b, ' The product b=A\'*x:' )
return
def r8cbb_mv ( n1, n2, ml, mu, a, x ):
#*****************************************************************************80
#
## R8CBB_MV multiplies a R8CBB matrix times a vector.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
# Input, real X(N1+N2), the vector to be multiplied by A.
#
# Output, real B(N1+N2), the result of multiplying A by X.
#
import numpy as np
#
# Set B to zero.
#
b = np.zeros ( n1 + n2 )
#
# Multiply by A1.
#
for j in range ( 0, n1 ):
ilo = max ( 0, j - mu )
ihi = min ( n1 - 1, j + ml )
ij = j * ( ml + mu + 1 ) - j + mu
for i in range ( ilo, ihi + 1 ):
b[i] = b[i] + a[i+ij] * x[j]
#
# Multiply by A2.
#
for j in range ( n1, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1
for i in range ( 0, n1 ):
b[i] = b[i] + a[i+ij] * x[j]
#
# Multiply by A3 and A4.
#
for j in range ( 0, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + n1 * n2 + j * n2 - n1
for i in range ( n1, n1 + n2 ):
b[i] = b[i] + a[i+ij] * x[j]
return b
def r8cbb_mv_test ( ):
#*****************************************************************************80
#
## R8CBB_MV_TEST tests R8CBB_MV.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 6
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_MV_TEST' )
print ( ' R8CBB_MV computes b=A*x, where A is an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_indicator ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix A:' )
x = r8vec_indicator1 ( n )
r8vec_print ( n, x, ' The vector x:' )
b = r8cbb_mv ( n1, n2, ml, mu, a, x )
r8vec_print ( n, b, ' The product b=A*x:' )
return
def r8cbb_print ( n1, n2, ml, mu, a, title ):
#*****************************************************************************80
#
## R8CBB_PRINT prints a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
#
# Input, string TITLE, a title to be printed.
#
r8cbb_print_some ( n1, n2, ml, mu, a, 0, 0, n1 + n2 - 1, n1 + n2 - 1, title )
return
def r8cbb_print_test ( ):
#*****************************************************************************80
#
## R8CBB_PRINT_TEST tests R8CBB_PRINT.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_PRINT_TEST' )
print ( ' R8CBB_PRINT prints an R8CBB matrix' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_random ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix:' )
return
def r8cbb_print_some ( n1, n2, ml, mu, a, ilo, jlo, ihi, jhi, title ):
#*****************************************************************************80
#
## R8CBB_PRINT_SOME prints some of a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
#
# Input, integer ILO, JLO, IHI, JHI, the first row and
# column, and the last row and column to be printed.
#
# Input, string TITLE, a title.
#
print ( '' )
print ( title )
incx = 5
#
# Print the columns of the matrix, in strips of 5.
#
for j2lo in range ( jlo, jhi + 1, incx ):
j2hi = j2lo + incx - 1
j2hi = min ( j2hi, n1 + n2 - 1 )
j2hi = min ( j2hi, jhi )
print ( '' )
print ( ' Col: ', end = '' )
for j in range ( j2lo, j2hi + 1 ):
print ( '%7d ' % ( j ), end = '' )
print ( '' )
print ( ' Row' )
#
# Determine the range of the rows in this strip.
#
i2lo = max ( ilo, 0 )
i2hi = min ( ihi, n1 + n2 - 1 )
for i in range ( i2lo, i2hi + 1 ):
print ( '%7d :' % ( i ), end = '' )
for j in range ( j2lo, j2hi + 1 ):
#
# Print out (up to) 5 entries in row I, that lie in the current strip.
#
aij = r8cbb_get ( n1, n2, ml, mu, a, i, j )
print ( '%12g ' % ( aij ), end = '' )
print ( '' )
return
def r8cbb_print_some_test ( ):
#*****************************************************************************80
#
## R8CBB_PRINT_SOME_TEST tests R8CBB_PRINT_SOME.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_PRINT_SOME_TEST' )
print ( ' R8CBB_PRINT_SOME prints some of an R8CBB matrix' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_random ( n1, n2, ml, mu )
r8cbb_print_some ( n1, n2, ml, mu, a, 1, 9, 10, 10, ' Rows 1-10, Cols 9-10' )
return
def r8cbb_random ( n1, n2, ml, mu ):
#*****************************************************************************80
#
## R8CBB_RANDOM randomizes a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative and no greater than N1-1.
#
# Output, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
a = np.zeros ( ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 )
for j in range ( 0, n1 ):
for row in range ( 0, ml + mu + 1 ):
i = row + j - mu
if ( 0 <= i and i < n1 ):
a[row+j*(ml+mu+1)] = rng.random ( size = 1 )
#
# Randomize the rectangular strips A2+A3+A4.
#
klo = ( ml + mu + 1 ) * n1
khi = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
for k in range ( klo, khi ):
a[k] = rng.random ( size = 1 )
return a
def r8cbb_random_test ( ):
#*****************************************************************************80
#
## R8CBB_RANDOM_TEST tests R8CBB_RANDOM.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_RANDOM_TEST' )
print ( ' R8CBB_RANDOM generates a random R8CBB matrix' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_random ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix:' )
return
def r8cbb_set ( n1, n2, ml, mu, a, i, j, value ):
#*****************************************************************************80
#
## R8CBB_SET sets the value of an entry in a R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 28 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
# Input, integer I, J, the row and column of the entry to set.
#
# Input, real VALUE, the value to be assigned to the (I,J) entry.
#
# Output, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the updated R8CBB matrix.
#
#
# Check for I or J out of bounds.
#
if ( i < 0 or n1 + n2 <= i ):
print ( '' )
print ( 'R8CBB_SET(): Fatal error!' )
print ( ' Illegal input value of row index I = ', i )
raise Exception ( 'R8CBB_SET(): Fatal error!' )
if ( j < 0 or n1 + n2 <= j ):
print ( '' )
print ( 'R8CBB_SET(): Fatal error!' )
print ( ' Illegal input value of column index J = ', j )
raise Exception ( 'R8CBB_SET(): Fatal error!' )
#
# The A1 block of the matrix.
#
# Check for out of band problems.
#
if ( i < n1 and j < n1 ):
if ( mu < j - i or ml < i - j ):
print ( '' )
print ( 'R8CBB_SET(): Fatal error!' )
print ( ' Unable to set entry (', i, ',', j, ').' )
raise Exception ( 'R8CBB_SET(): Fatal error!' )
ij = ( i - j + mu ) + j * ( ml + mu + 1 )
#
# The A2 block of the matrix:
#
elif ( i < n1 and n1 <= j ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i
#
# The A3 and A4 blocks of the matrix.
#
elif ( n1 <= i ):
ij = ( ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 )
else:
print ( '' )
print ( 'R8CBB_SET(): Fatal error!' )
print ( ' Unable to set entry (', i, ',', j, ').' )
raise Exception ( 'R8CBB_SET(): Fatal error!' )
a[ij] = value
return a
def r8cbb_set_test ( ):
#*****************************************************************************80
#
## R8CBB_SET_TEST tests R8CBB_SET.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 4
n2 = 1
n = n1 + n2
ml = 2
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_SET_TEST' )
print ( ' R8CBB_SET sets the value of an element of an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Initialize matrix to zero.
#
a = r8cbb_zeros ( n1, n2, ml, mu )
#
# Set matrix entries to the "indicator" values.
#
for i in range ( 0, n1 ):
jlo = max ( 0, i - ml )
jhi = min ( n1, i + mu )
for j in range ( jlo, jhi + 1 ):
value = float ( 10 * ( i + 1 ) + ( j + 1 ) )
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
for i in range ( 0, n1 ):
for j in range ( n1, n1 + n2 ):
value = float ( 10 * ( i + 1 ) + ( j + 1 ) )
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
for i in range ( n1, n1 + n2 ):
for j in range ( 0, n1 ):
value = float ( 10 * ( i + 1 ) + ( j + 1 ) )
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
for i in range ( n1, n1 + n2 ):
for j in range ( n1, n1 + n2 ):
value = float ( 10 * ( i + 1 ) + ( j + 1 ) )
a = r8cbb_set ( n1, n2, ml, mu, a, i, j, value )
r8cbb_print ( n1, n2, ml, mu, a, ' The matrix after additions:' )
return
def r8cbb_sl ( n1, n2, ml, mu, a_lu, b ):
#*****************************************************************************80
#
## R8CBB_SL solves a R8CBB system factored by R8CBB_FA.
#
# Discussion:
#
# Note that in C++ and FORTRAN, we can look at A as an abstract
# vector, but then look at parts of A as storing a two dimensional
# array. MATLAB assigns an inherent dimensionality to a data object,
# and gets very unhappy when you try to manipulate the data yourself.
# This means that the MATLAB implementation of this routine requires
# the use of temporary 2D arrays.
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# The linear system A * x = b is decomposable into the block system:
#
# ( A1 A2 ) * (X1) = (B1)
# ( A3 A4 ) (X2) (B2)
#
# where A1 is a (usually big) banded square matrix, A2 and A3 are
# column and row strips which may be nonzero, and A4 is a dense
# square matrix.
#
# All the arguments except B are input quantities only, which are
# not changed by the routine. They should have exactly the same values
# they had on return from R8CBB_FA.
#
# If more than one right hand side is to be solved, with the same
# matrix, R8CBB_SL should be called repeatedly. However, R8CBB_FA only
# needs to be called once to create the factorization.
#
# See the documentation of R8CBB_FA for details on the matrix storage.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 30 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A_LU( (ML+MU+1)*N1 + 2*N1*N2 + N2*N2).
# the compact border banded matrix, as factored by R8CBB_FA.
#
# Input, real B(N1+N2), the right hand side of the linear system.
#
# Output, real X(N1+N2), the solution.
#
import numpy as np
x = np.zeros ( n1 + n2 )
nband = ( ml + mu + 1 ) * n1
#
# Set B1 := inverse(A1) * B1.
#
if ( 0 < n1 ):
a1_lu = r8vec_to_r8cb ( n1, n1, ml, mu, a_lu[0:nband] )
job = 0
x[0:n1] = r8cb_np_sl ( n1, ml, mu, a1_lu, b[0:n1], job )
#
# Modify the right hand side of the second linear subsystem.
# Replace B2 by B2-A3*B1.
#
for j in range ( 0, n1 ):
for i in range ( 0, n2 ):
ij = nband + n1 * n2 + j * n2 + i
b[n1+i] = b[n1+i] - a_lu[ij] * x[j]
#
# Solve A4*B2 = B2.
#
if ( 0 < n2 ):
a4_lu = r8vec_to_r8ge ( n2, n2, a_lu[nband+2*n1*n2:nband+2*n1*n2+n2*n2] )
job = 0
x[n1:n1+n2] = r8ge_np_sl ( n2, a4_lu, b[n1:n1+n2], job )
#
# Modify the first subsolution.
# Set B1 = B1+A2*B2.
#
for i in range ( 0, n1 ):
for j in range ( 0, n2 ):
ij = nband + j * n1 + i
x[i] = x[i] + a_lu[ij] * x[n1+j]
return x
def r8cbb_sl_test ( ):
#*****************************************************************************80
#
## R8CBB_SL_TEST tests R8CBB_SL.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_SL_TEST' )
print ( ' R8CBB_SL solves a linear system factored by R8CBB_FA' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
#
# Set the matrix.
#
a = r8cbb_random ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix:' )
#
# Set the desired solution.
#
x = r8vec_indicator1 ( n )
#
# Compute the corresponding right hand side.
#
b = r8cbb_mv ( n1, n2, ml, mu, a, x )
#
# Factor the matrix
#
a_lu, info = r8cbb_fa ( n1, n2, ml, mu, a )
r8cbb_print ( n1, n2, ml, mu, a_lu, ' The factored matrix:' )
if ( info != 0 ):
print ( '' )
print ( 'R8CBB_SL_TEST(): Fatal error!' )
print ( ' R8CBB_FA claims the matrix is singular.' )
print ( ' The value of INFO is ', info )
raise Exception ( 'R8CBB_SL_TEST(): Fatal error!' )
#
# Solve the system.
#
r8vec_print ( n, b, ' The right hand side vector:' )
x = r8cbb_sl ( n1, n2, ml, mu, a_lu, b )
r8vec_print ( n, x, ' Solution:' )
return
def r8cbb_to_r8ge ( n1, n2, ml, mu, a ):
#*****************************************************************************80
#
## R8CBB_TO_R8GE copies a R8CBB matrix to a R8GE matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= n1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative, and no greater than N1-1.
#
# Input, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
#
# Output, real B(N1+N2,N1+N2), the R8GE matrix.
#
import numpy as np
b = np.zeros ( [ n1 + n2, n1 + n2 ] )
for i in range ( 0, n1 ):
for j in range ( 0, n1 ):
if ( i - ml <= j and j <= i + mu ):
ij = ( i - j + mu ) + j * ( ml + mu + 1 )
b[i,j] = a[ij]
for i in range ( 0, n1 ):
for j in range ( n1, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i
b[i,j] = a[ij]
for i in range ( n1, n1 + n2 ):
for j in range ( 0, n1 + n2 ):
ij = ( ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 )
b[i,j] = a[ij]
return b
def r8cbb_to_r8ge_test ( ):
#*****************************************************************************80
#
## R8CBB_TO_R8GE_TEST tests R8CBB_TO_R8GE.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_TO_R8GE_TEST' )
print ( ' R8CBB_TO_R8GE converts an R8CBB matrix to R8GE format.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
a = r8cbb_indicator ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB matrix:' )
a_r8ge = r8cbb_to_r8ge ( n1, n2, ml, mu, a )
print ( '' )
print ( ' The R8GE matrix:' )
print ( a_r8ge )
return
def r8cbb_zeros ( n1, n2, ml, mu ):
#*****************************************************************************80
#
## R8CBB_ZEROS zeros an R8CBB matrix.
#
# Discussion:
#
# The R8CBB storage format is for a compressed border banded matrix.
# Such a matrix has the logical form:
#
# A1 | A2
# ---+---
# A3 | A4
#
# with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
# are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
# respectively.
#
# The R8CBB format is the same as the DBB format, except that the banded
# matrix A1 is stored in compressed band form rather than standard
# banded form. In other words, we do not include the extra room
# set aside for fill in during pivoting.
#
# A should be defined as a vector. The user must then store
# the entries of the four blocks of the matrix into the vector A.
# Each block is stored by columns.
#
# A1, the banded portion of the matrix, is stored in
# the first (ML+MU+1)*N1 entries of A, using the obvious variant
# of the LINPACK general band format.
#
# The following formulas should be used to determine how to store
# the entry corresponding to row I and column J in the original matrix:
#
# Entries of A1:
#
# 1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
#
# Store the I, J entry into location
# (I-J+MU+1)+(J-1)*(ML+MU+1).
#
# Entries of A2:
#
# 1 <= I <= N1, N1+1 <= J <= N1+N2.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+(J-N1-1)*N1+I.
#
# Entries of A3:
#
# N1+1 <= I <= N1+N2, 1 <= J <= N1.
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
#
# Entries of A4:
#
# N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
#
# Store the I, J entry into location
# (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
# (same formula used for A3).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N1, N2, the order of the banded and dense blocks.
# N1 and N2 must be nonnegative, and at least one must be positive.
#
# Input, integer ML, MU, the lower and upper bandwidths.
# ML and MU must be nonnegative and no greater than N1-1.
#
# Output, real A((ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8CBB matrix.
#
import numpy as np
nn = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
a = np.zeros ( nn )
return a
def r8cbb_zeros_test ( ):
#*****************************************************************************80
#
## R8CBB_ZEROS_TEST tests R8CBB_ZEROS.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 25 July 2016
#
# Author:
#
# John Burkardt
#
n1 = 8
n2 = 2
n = n1 + n2
ml = 1
mu = 1
na = ( ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2
print ( '' )
print ( 'R8CBB_ZEROS_TEST' )
print ( ' R8CBB_ZEROS zeros an R8CBB matrix.' )
print ( '' )
print ( ' Matrix order N = ', n )
print ( ' Matrix suborder N1 = ', n1 )
print ( ' Matrix suborder N2 = ', n2 )
print ( ' Lower bandwidth ML = ', ml )
print ( ' Upper bandwidth MU = ', mu )
a = r8cbb_zeros ( n1, n2, ml, mu )
r8cbb_print ( n1, n2, ml, mu, a, ' The R8CBB zero matrix:' )
return
def r8ge_np_fa ( n, a ):
#*****************************************************************************80
#
## r8ge_np_fa() factors a R8GE matrix by nonpivoting Gaussian elimination.
#
# 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_np_fa is a version of the LINPACK routine R8GEFA, but uses no
# pivoting. It will fail if the matrix is singular, or if any zero
# pivot is encountered.
#
# If r8ge_np_fa successfully factors the matrix, R8GE_NP_SL may be called
# to solve linear systems involving the matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 June 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be positive.
#
# Input, real A(N,N), the matrix to be factored.
#
# Output, real A_LU(N,N), information about the factorization,
# which must be passed unchanged to R8GE_NP_SL for solutions.
#
# Output, integer INFO, singularity flag.
# 0, no singularity detected.
# nonzero, the factorization failed on the INFO-th step.
#
a_lu = a.copy ( )
info = 0
for k in range ( 1, n ):
if ( a_lu[k-1,k-1] == 0.0 ):
info = k
print ( '' )
print ( 'r8ge_np_fa(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'r8ge_np_fa(): Fatal error!' )
for i in range ( k + 1, n + 1 ):
a_lu[i-1,k-1] = - a_lu[i-1,k-1] / a_lu[k-1,k-1]
for j in range ( k + 1, n + 1 ):
for i in range ( k + 1, n + 1 ):
a_lu[i-1,j-1] = a_lu[i-1,j-1] + a_lu[i-1,k-1] * a_lu[k-1,j-1]
if ( a_lu[n-1,n-1] == 0.0 ):
info = n
print ( '' )
print ( 'r8ge_np_fa(): Fatal error!' )
print ( ' Zero pivot on step ', info )
raise Exception ( 'r8ge_np_fa(): Fatal error!' )
return a_lu, info
def r8ge_np_sl ( n, a_lu, b, job ):
#*****************************************************************************80
#
## r8ge_np_sl() solves a system factored by r8ge_np_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.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 06 June 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
# N must be positive.
#
# Input, real A_LU(N,N), the matrix as factored by r8ge_np_fa.
#
# Input, real B(N), the right hand side vector B.
#
# Input, integer JOB.
# If JOB is zero, the routine will solve A * x = b.
# If JOB is nonzero, the routine will solve A' * x = b.
#
# Output, real X(N), the solution.
#
x = b.copy ( )
#
# Solve A * x = b.
#
if ( job == 0 ):
for k in range ( 1, n ):
for i in range ( k + 1, n + 1 ):
x[i-1] = x[i-1] + a_lu[i-1,k-1] * x[k-1]
for k in range ( n, 0, -1 ):
x[k-1] = x[k-1] / a_lu[k-1,k-1]
for i in range ( 1, k ):
x[i-1] = x[i-1] - a_lu[i-1,k-1] * x[k-1]
#
# Solve A' * X = B.
#
else:
for k in range ( 1, n + 1 ):
t = 0.0
for i in range ( 1, k ):
t = t + x[i-1] * a_lu[i-1,k-1]
x[k-1] = ( x[k-1] - t ) / a_lu[k-1,k-1]
for k in range ( n - 1, 0, -1 ):
for i in range ( k + 1, n + 1 ):
x[k-1] = x[k-1] + x[i-1] * a_lu[i-1,k-1]
return x
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 r8vec_indicator1 ( n ):
#*****************************************************************************80
#
## r8vec_indicator1() sets an R8VEC to the indicator vector (1,2,3,...).
#
# Discussion:
#
# An R8VEC is a vector of R8's.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 27 September 2014
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the number of elements of the vector.
#
# Output:
#
# real A(N), the indicator array.
#
import numpy as np
a = np.zeros ( n );
for i in range ( 0, n ):
a[i] = i + 1
return a
def r8vec_print ( n, a, title ):
#*****************************************************************************80
#
## r8vec_print() prints an R8VEC.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 31 August 2014
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer N, the dimension of the vector.
#
# real A(N), the vector to be printed.
#
# string TITLE, a title.
#
print ( '' )
print ( title )
print ( '' )
for i in range ( 0, n ):
print ( '%6d: %12g' % ( i, a[i] ) )
return
def r8vec_to_r8cb ( m, n, ml, mu, x ):
#*****************************************************************************80
#
## r8vec_to_r8cb() copies an R8VEC into a R8CB 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:
#
# 26 July 2016
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer M, N, the number of rows and columns in the array.
#
# Input, integer ML, MU, the lower and upper bandwidths.
#
# Input, real X((ML+MU+1)*N), the vector to be copied into the array.
#
# Output, real A(ML+MU+1,N), the array.
#
import numpy as np
a = np.zeros ( [ ml + mu + 1, n ] )
for j in range ( 0, n ):
for k in range ( 0, ml + mu + 1 ):
i = k + j - mu
if ( 0 <= i and i < m ):
a[k,j] = x[k+j*(ml+mu+1)]
return a
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 timestamp ( ):
#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 August 2019
#
# Author:
#
# John Burkardt
#
import time
t = time.time ( )
print ( time.ctime ( t ) )
return
if ( __name__ == '__main__' ):
timestamp ( )
r8cbb_test ( )
timestamp ( )