#! /usr/bin/env python3 # def r8bb_test ( ): #*****************************************************************************80 # ## r8bb_test() tests r8bb(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 August 2022 # # Author: # # John Burkardt # print ( '' ) print ( 'r8bb_test():' ) print ( ' Python version:' ) print ( ' Test r8bb().' ) r8bb_add_test ( ) r8bb_dif2_test ( ) r8bb_fa_test ( ) r8bb_get_test ( ) r8bb_indicator_test ( ) r8bb_mtv_test ( ) r8bb_mv_test ( ) r8bb_print_test ( ) r8bb_print_some_test ( ) r8bb_random_test ( ) r8bb_set_test ( ) r8bb_sl_test ( ) r8bb_to_r8ge_test ( ) r8bb_zeros_test ( ) # # Terminate. # print ( '' ) print ( 'r8bb_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 # # Parameters: # # 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 = np.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 r8bb_add ( n1, n2, ml, mu, a, i, j, value ): #*****************************************************************************80 # ## r8bb_add() adds a value to an entry in a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Input, integer I, J, the row and column of the entry to be incremented. # Some combinations of I and J are illegal. # # Input, real VALUE, the value to be added to the (I,J)-th entry. # # Output, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the updated R8BB matrix. # if ( value == 0.0 ): return a if ( i < 0 or n1 + n2 <= i ): print ( '' ) print ( 'R8BB_ADD - Fatal error!' ) print ( ' Illegal input value of row index I = ', i ) raise Exception ( 'R8BB_ADD - Fatal error!' ) if ( j < 0 or n1 + n2 <= j ): print ( '' ) print ( 'R8BB_ADD - Fatal error!' ) print ( ' Illegal input value of column index J = ', j ) raise Exception ( 'R8BB_ADD - Fatal error!' ) # # The A1 block of the matrix. # # Check for out of band problems. # # Normally, we would check the condition MU < (J-I), but the storage # format requires extra entries be set aside in case of pivoting, which # means that the condition becomes MU+ML < (J-I). # if ( i < n1 and j < n1 ): if ( mu + ml < j - i or ml < i - j ): print ( '' ) print ( 'R8BB_ADD - Fatal error!' ) print ( ' Unable to add to entry (%d,%d).' % ( i, j ) ) exit ( ) else: ij = ( i - j + ml + mu ) + j * ( 2 * ml + mu + 1 ) # # The A2 block of the matrix. # elif ( i < n1 and n1 <= j ): ij = ( 2 * ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i # # The A3 and A4 blocks of the matrix. # elif ( n1 <= i ): ij = ( 2 * ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 ); a[ij] = a[ij] + value return a def r8bb_add_test ( ): #*****************************************************************************80 # ## R8BB_ADD_TEST tests R8BB_ADD. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 0 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_ADD_TEST' ) print ( ' R8BB_ADD adds a value to elements of an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Initialize matrix to indicator matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) # # Print initial matrix. # r8bb_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 = r8bb_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 = r8bb_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 = r8bb_add ( n1, n2, ml, mu, a, i, j, value ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix after additions:' ) return def r8bb_dif2 ( n1, n2, ml, mu ): #*****************************************************************************80 # ## R8BB_DIF2 sets up an R8BB second difference matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # import numpy as np a = np.zeros ( ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 ) if ( ml < 1 or mu < 1 ): print ( '' ) print ( 'R8BB_DIF2 - Fatal error!' ) print ( ' 1 <= ML and 1 <= MU required.' ) raise Exception ( 'R8BB_DIF2 - Fatal error!' ) for i in range ( 1, n1 + n2 ): j = i - 1 value = - 1.0 a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) for i in range ( 0, n1 + n2 ): j = i value = 2.0 a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) for i in range ( 0, n1 + n2 - 1 ): j = i + 1 value = - 1.0 a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) return a def r8bb_dif2_test ( ): #*****************************************************************************80 # ## R8BB_DIF2_TEST tests R8BB_DIF2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_DIF2_TEST' ) print ( ' R8BB_DIF2 sets up an R8BB second difference matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_dif2 ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB second difference matrix:' ) return def r8bb_fa ( n1, n2, ml, mu, a ): #*****************************************************************************80 # ## R8BB_FA factors a R8BB 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 R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1). # (same formula used for A3). # # Once the matrix has been factored by R8BB_FA, R8BB_SL may be called # to solve linear systems involving the matrix. # # R8BB_FA uses LINPACK routines to carry out the factorization. # # # 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 + inv(A1) A2 X2 = inv(A1) B1 # # A3 X1 + A4 X2 = B2 # # and then rewrites the second equation as # # ( A4 - A3 inv(A1) A2 ) X2 = B2 - A3 inv(A1) B1 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 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( (2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2 ), the border-banded # matrix to be factored. # # Output, real A_LU( (2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2 ). # Information describing a partial factorization # of the original coefficient matrix. This information is required # by R8BB_SL in order to solve linear systems associated with that # matrix. # # Output, integer PIVOT(N1+N2), contains pivoting information. # # Output, integer INFO, singularity flag. # 0, no singularity detected. # nonzero, the factorization failed on the INFO-th step. # import numpy as np nband = ( 2 * ml + mu + 1 ) * n1 a_lu = a.copy ( ) pivot = np.zeros ( n1 + n2, dtype = np.int32 ) # # Copy the (2*ML+MU+1 x N1 ) A1 band matrix out of A. # # Factor the A1 band matrix, overwriting A1 by its factors. # if ( 0 < n1 ): a1 = r8vec_to_r8gb ( n1, n1, ml, mu, a[0:nband] ) a1_lu, pivot[0:n1], info = r8gb_fa ( n1, ml, mu, a1 ) if ( info != 0 ): print ( '' ) print ( 'R8BB_FA - Fatal error!' ) print ( ' R8GB_FA returned INFO = %d' % ( info ) ) print ( ' Factoring failed for column INFO.' ) print ( ' The band matrix A1 is singular.' ) print ( ' This algorithm cannot continue!' ) raise Exception ( 'R8BB_FA(): Fatal error!' ) a_lu[0:nband] = r8gb_to_r8vec ( n1, n1, ml, mu, a1_lu ) if ( 0 < n1 and 0 < n2 ): # # Solve A1 * x = -A2 for x, and overwrite A2 by the results. # job = 0 for j in range ( 0, n2 ): b2 = - a[nband+j*n1:nband+j*n1+n1] x2 = r8gb_sl ( n1, ml, mu, a1_lu, pivot, b2, job ) a_lu[nband+j*n1:nband+j*n1+n1] = x2[0:n1] # # 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 temp = a_lu[ik] 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, pivot2, info = r8ge_fa ( n2, a4 ) pivot[n1:n1+n2] = pivot2[0:n2] if ( info != 0 ): print ( '' ) print ( 'R8BB_FA - Fatal error!' ) print ( ' R8GE_FA returned INFO = %d' % (info ) ) print ( ' This indicates singularity in column INFO.' ) print ( ' of the A4 submatrix, which is column %d' % ( n1 + info ) ) print ( ' of the full matrix.' ) print ( '' ) print ( ' It is possible that the full matrix is' ) print ( ' nonsingular, but the algorithm R8BB_FA may' ) print ( ' not be used for this matrix.' ) exit ( ) a_lu[nband+2*n1*n2:nband+2*n1*n2+n2*n2] = r8ge_to_r8vec ( n2, n2, a4_lu ) return a_lu, pivot, info def r8bb_fa_test ( ): #*****************************************************************************80 # ## R8BB_FA_TEST tests R8BB_FA. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 seed = 123456789 print ( '' ) print ( 'R8BB_FA_TEST' ) print ( ' R8BB_FA factors an R8BB matrix' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a, seed = r8bb_random ( n1, n2, ml, mu, seed ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8bb_mv ( n1, n2, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8bb_fa ( n1, n2, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8BB_FA_TEST - Fatal error!' ) print ( ' R8BB_FA claims the matrix is singular.' ) print ( ' The value of INFO is %d' % ( info ) ) exit ( ) # # Solve the system. # x = r8bb_sl ( n1, n2, ml, mu, a_lu, pivot, b ) r8vec_print ( n, x, ' Solution:' ) return def r8bb_get ( n1, n2, ml, mu, a, i, j ): #*****************************************************************************80 # ## R8BB_GET returns an entry of a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Input, integer I, J, the row and column of the entry to be retrieved. # # Output, real VALUE, the value of the (I,J) entry. # value = 0.0 if ( i < 0 or n1 + n2 <= i ): print ( '' ) print ( 'R8BB_GET - Fatal error!' ) print ( ' Illegal input value of row index I = %d' % ( i ) ) raise Exception ( 'R8BB_GET - Fatal error!' ) if ( j < 0 or n1 + n2 <= j ): print ( '' ) print ( 'R8BB_GET - Fatal error!' ) print ( ' Illegal input value of colum index J = %d' % ( j ) ) raise Exception ( 'R8BB_GET - Fatal error!' ) # # The A1 block of the matrix. # # Check for out of band problems. # # Normally, we would check the condition MU < (J-I), but the storage # format requires extra entries be set aside in case of pivoting, which # means that the condition becomes MU+ML < (J-I). # if ( i < n1 and j < n1 ): if ( mu + ml < j - i or ml < i - j ): value = 0.0 return value else: ij = ( i - j + ml + mu ) + j * ( 2 * ml + mu + 1 ) # # The A2 block of the matrix. # elif ( i < n1 and n1 <= j ): ij = ( 2 * ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i # # The A3 and A4 blocks of the matrix. # elif ( n1 <= i ): ij = ( 2 * ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 ) value = a[ij] return value def r8bb_get_test ( ): #*****************************************************************************80 # ## R8BB_GET_TEST tests R8BB_GET. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # from numpy.random import default_rng import numpy as np rng = default_rng ( ) n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 0 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_GET_TEST' ) print ( ' R8BB_GET gets a value of an element of an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set matrix to indicator matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) # # Print matrix. # r8bb_print ( n1, n2, ml, mu, a, ' The matrix to be queried:' ) # # Request random entries. # print ( '' ) for k in range ( 0, 10 ): i = rng.integers ( low = 0, high = n1 + n2, size = 1, endpoint = False ) j = rng.integers ( low = 0, high = n1 + n2, size = 1, endpoint = False ) value = r8bb_get ( n1, n2, ml, mu, a, i, j ) print ( ' A(%d,%d) = %g' % ( i, j, value ) ) return def r8bb_indicator ( n1, n2, ml, mu ): #*****************************************************************************80 # ## R8BB_INDICATOR sets up a R8BB indicator matrix. # # Discussion: # # The R8BB storage format is for a 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. # # Example: # # With N1 = 4, N2 = 1, ML = 1, MU = 2, the matrix entries would be: # # 00 # 00 00 # 00 00 00 --- --- # A11 A12 A13 00 --- A16 A17 # A21 A22 A23 A24 00 A26 A27 # --- A32 A33 A34 A35 A36 A37 # --- --- A43 A44 A45 A46 A47 # --- --- --- A54 A55 A56 A57 # 00 # # A61 A62 A63 A64 A65 A66 A67 # A71 A72 A73 A74 A75 A76 A77 # # The matrix is actually stored as a vector, and we will simply suggest # the structure and values of the indicator matrix as: # # 00 00 00 00 00 # 00 00 13 24 35 16 17 61 62 63 64 65 66 67 # 00 12 23 34 45 + 26 27 + 71 72 73 74 75 + 76 77 # 11 22 33 44 55 36 37 # 21 32 43 54 00 46 47 # 56 57 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # import numpy as np a = np.zeros ( ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 ) fac = 10 ** ( i4_log_10 ( n1 + n2 ) + 1 ) # # Set the banded matrix A1. # for i in range ( 0, n1 ): jlo = max ( i - ml, 0 ) jhi = min ( i + mu, n1 ) for j in range ( jlo, jhi + 1 ): value = float ( fac * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Set the N1 by N2 rectangular strip A2. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): value = float ( fac * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Set the N2 by N1 rectangular strip A3. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): value = float ( fac * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Set the N2 by N2 square A4. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): value = float ( fac * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) return a def r8bb_indicator_test ( ): #*****************************************************************************80 # ## R8BB_INDICATOR_TEST tests R8BB_INDICATOR. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_INDICATOR_TEST' ) print ( ' R8BB_INDICATOR sets up an indicator matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The indicator matrix:' ) return def r8bb_mtv ( n1, n2, ml, mu, a, x ): #*****************************************************************************80 # ## R8BB_MTV multiplies a vector by a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8BB matrix. # # Input, real X(N1+N2), the vector to multiply A. # # Output, real B(N1+N2), the product X times A. # import numpy as np # # Initialize B. # b = np.zeros ( n1 + n2 ) # # Multiply by A1. # for i in range ( 0, n1 ): jlo = max ( i - ml, 0 ) jhi = min ( i + mu, n1 ) for j in range ( jlo, jhi + 1 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[j] = b[j] + x[i] * aij # # Multiply by A2. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[j] = b[j] + x[i] * aij # # Multiply by A3. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[j] = b[j] + x[i] * aij # # Multiply by A4. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[j] = b[j] + x[i] * aij return b def r8bb_mtv_test ( ): #*****************************************************************************80 # ## R8BB_MTV_TEST tests R8BB_MTV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 6 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_MTV_TEST' ) print ( ' R8BB_MTV computes b=A''*x, where A is an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix A:' ) x = r8vec_indicator1 ( n ) r8vec_print ( n, x, ' The vector x:' ) b = r8bb_mtv ( n1, n2, ml, mu, a, x ) r8vec_print ( n, b, ' The product b=A''*x:' ) return def r8bb_mv ( n1, n2, ml, mu, a, x ): #*****************************************************************************80 # ## R8BB_MV multiplies a R8BB matrix times a vector. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB 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 # # Initialize B. # b = np.zeros ( n1 + n2 ) # # Multiply by A1. # for i in range ( 0, n1 ): jlo = max ( i - ml, 0 ) jhi = min ( i + mu, n1 ) for j in range ( jlo, jhi + 1 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[i] = b[i] + aij * x[j] # # Multiply by A2. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[i] = b[i] + aij * x[j] # # Multiply by A3. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[i] = b[i] + aij * x[j] # # Multiply by A4. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) b[i] = b[i] + aij * x[j] return b def r8bb_mv_test ( ): #*****************************************************************************80 # ## R8BB_MV_TEST tests R8BB_MV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 6 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_MV_TEST' ) print ( ' R8BB_MV computes b=A*x, where A is an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix A:' ) x = r8vec_indicator1 ( n ) r8vec_print ( n, x, ' The vector x:' ) b = r8bb_mv ( n1, n2, ml, mu, a, x ) r8vec_print ( n, b, ' The product b=A*x:' ) return def r8bb_print ( n1, n2, ml, mu, a, title ): #*****************************************************************************80 # ## R8BB_PRINT prints a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Input, string TITLE, a title to be printed. # r8bb_print_some ( n1, n2, ml, mu, a, 0, 0, n1 + n2 - 1, n1 + n2 - 1, title ) return def r8bb_print_test ( ): #*****************************************************************************80 # ## R8BB_PRINT_TEST tests R8BB_PRINT. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 8 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 seed = 123456789 print ( '' ) print ( 'R8BB_PRINT_TEST' ) print ( ' R8BB_PRINT prints an R8BB matrix' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a, seed = r8bb_random ( n1, n2, ml, mu, seed ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix:' ) return def r8bb_print_some ( n1, n2, ml, mu, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## R8BB_PRINT_SOME prints some of a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB 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 ) inc = j2hi + 1 - j2lo 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 = '' ) # # Print out (up to) 5 entries in row I, that lie in the current strip. # for j in range ( j2lo, j2hi + 1 ): aij = r8bb_get ( n1, n2, ml, mu, a, i, j ) print ( '%12g ' % ( aij ), end = '' ) print ( '' ) return def r8bb_print_some_test ( ): #*****************************************************************************80 # ## R8BB_PRINT_SOME_TEST tests R8BB_PRINT_SOME. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 6 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_PRINT_SOME_TEST' ) print ( ' R8BB_PRINT_SOME prints some of an R8BB matrix' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) r8bb_print_some ( n1, n2, ml, mu, a, 6, 6, 7, 7, ' The Lower Right Block:' ) return def r8bb_random ( n1, n2, ml, mu, seed ): #*****************************************************************************80 # ## R8BB_RANDOM randomizes a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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, integer SEED, a seed for the random number generator. # # Output, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Output, integer SEED, an updated seed for the random number generator. # from numpy.random import default_rng import numpy as np rng = default_rng ( ) a = np.zeros ( ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 ) # # A1. # for i in range ( 0, n1 ): jlo = max ( i - ml, 0 ) jhi = min ( i + mu, n1 ) for j in range ( jlo, jhi + 1 ): aij = rng.random ( ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, aij ) # # A2. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): aij = rng.random ( ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, aij ) # # A3. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): aij = rng.random ( ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, aij ) # # A4. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): aij = rng.random ( ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, aij ) return a, seed def r8bb_random_test ( ): #*****************************************************************************80 # ## R8BB_RANDOM_TEST tests R8BB_RANDOM. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 seed = 123456789 print ( '' ) print ( 'R8BB_RANDOM_TEST' ) print ( ' R8BB_RANDOM returns a random R8BB matrix' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a, seed = r8bb_random ( n1, n2, ml, mu, seed ) r8bb_print ( n1, n2, ml, mu, a, ' The border-banded matrix:' ) return def r8bb_set ( n1, n2, ml, mu, a, i, j, value ): #*****************************************************************************80 # ## R8BB_SET sets an entry of a R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Input, integer I, J, the row and column of the entry to be set. # # Input, real VALUE, the value to be assigned to the (I,J) entry. # # Output, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the updated R8BB matrix. # if ( i < 0 or n1 + n2 <= i ): print ( '' ) print ( 'R8BB_SET - Fatal error!' ) print ( 'R8BB_SET - Illegal value of row index I = %d' % ( i ) ) raise Exception ( 'R8BB_SET(): Fatal error!' ) if ( j < 0 or n1 + n2 <= j ): print ( '' ) print ( 'R8BB_SET - Fatal error!' ) print ( 'R8BB_SET - Illegal value of column index J = %d' % ( j ) ) raise Exception ( 'R8BB_SET(): Fatal error!' ) # # The A1 block of the matrix. # # Check for out of band problems. # # Normally, we would check the condition MU < (J-I), but the storage # format requires extra entries be set aside in case of pivoting, which # means that the condition becomes MU+ML < (J-I). # if ( i < n1 and j < n1 ): if ( mu + ml < j - i or ml < i - j ): print ( '' ) print ( 'R8BB_SET - Warning!' ) print ( 'R8BB_SET - Unable to set entry A(%d,%d).' % ( i, j ) ) exit ( ) else: ij = ( i - j + ml + mu ) + j * ( 2 * ml + mu + 1 ) # # The A2 block of the matrix. # elif ( i < n1 and n1 <= j ): ij = ( 2 * ml + mu + 1 ) * n1 + ( j - n1 ) * n1 + i # # The A3 and A4 blocks of the matrix. # elif ( n1 <= i ): ij = ( 2 * ml + mu + 1 ) * n1 + n2 * n1 + j * n2 + ( i - n1 ) a[ij] = value return a def r8bb_set_test ( ): #*****************************************************************************80 # ## R8BB_SET_TEST tests R8BB_SET. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 4 n2 = 1 n = n1 + n2 ml = 2 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_SET_TEST' ) print ( ' R8BB_SET sets elements of an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Initialize matrix to zero. # a = r8bb_zeros ( n1, n2, ml, mu ) # # Fill in band matrix. # for i in range ( 0, n1 ): for j in range ( 0, n1 ): if ( i - ml <= j and j <= i + mu ): value = float ( 10 * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Fill in right border vector. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): value = float ( 10 * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Fill in lower border vector. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): value = float ( 10 * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) # # Fill in lower right dense matrix. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): value = float ( 10 * ( i + 1 ) + ( j + 1 ) ) a = r8bb_set ( n1, n2, ml, mu, a, i, j, value ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix:' ) return def r8bb_sl ( n1, n2, ml, mu, a_lu, pivot, b ): #*****************************************************************************80 # ## R8BB_SL solves a R8BB system factored by R8BB_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 R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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) # # 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 exit from R8BB_FA. # # If more than one right hand side is to be solved, with the same matrix, # R8BB_SL should be called repeatedly. However, R8BB_FA only needs to be # called once to create the factorization. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 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( (2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the factor information # computed by R8BB_FA. # # Input, integer PIVOT(N1+N2), the pivoting information from R8BB_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 = ( 2 * ml + mu + 1 ) * n1 # # Set B1 := inverse(A1) * B1. # Copy the banded matrix out of A_LU and into A1_LU. # if ( 0 < n1 ): a1_lu = r8vec_to_r8gb ( n1, n1, ml, mu, a_lu[0:nband] ) job = 0 x[0:n1] = r8gb_sl ( n1, ml, mu, a1_lu, pivot[0:n1], b[0:n1], job ) # # Modify the right hand side of the second linear subsystem. # Set B2 := B2 - A3*B1. # for i in range ( 0, n2 ): for j in range ( 0, n1 ): ij = nband + n1 * n2 + j * n2 + i b[n1+i] = b[n1+i] - a_lu[ij] * x[j] # # Set B2 := inverse(A4) * B2. # Copy the dense matrix out of A_LU and into A4_LU. # 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_sl ( n2, a4_lu, pivot[n1:n1+n2], 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 r8bb_sl_test ( ): #*****************************************************************************80 # ## R8BB_SL_TEST tests R8BB_SL. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 8 n2 = 2 n = n1 + n2 ml = 0 mu = 0 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 seed = 123456789 print ( '' ) print ( 'R8BB_SL_TEST' ) print ( ' R8BB_SL solves a linear system factored by R8BB_FA.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a, seed = r8bb_random ( n1, n2, ml, mu, seed ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8bb_mv ( n1, n2, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8bb_fa ( n1, n2, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8BB_SL_TEST - Fatal error!' ) print ( ' R8BB_FA claims the matrix is singular.' ) print ( ' The value of INFO is %d' % ( info ) ) exit ( ) # # Solve the system. # x = r8bb_sl ( n1, n2, ml, mu, a_lu, pivot, b ) r8vec_print ( n, x, ' Solution:' ) return def r8bb_to_r8ge ( n1, n2, ml, mu, a ): #*****************************************************************************80 # ## R8BB_TO_R8GE copies a R8BB matrix to a R8GE matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Output, real B(N1+N2,N1+N2), the R8GE matrix. # import numpy as np b = np.zeros ( [ n1 + n2, n1 + n2 ] ) # # Fill in band matrix. # for i in range ( 0, n1 ): for j in range ( 0, n1 ): if ( i - ml <= j and j <= i + mu ): b[i,j] = r8bb_get ( n1, n2, ml, mu, a, i, j ) # # Fill in right border vector. # for i in range ( 0, n1 ): for j in range ( n1, n1 + n2 ): b[i,j] = r8bb_get ( n1, n2, ml, mu, a, i, j ) # # Fill in lower border vector. # for i in range ( n1, n1 + n2 ): for j in range ( 0, n1 ): b[i,j] = r8bb_get ( n1, n2, ml, mu, a, i, j ) # # Fill in lower right dense matrix. # for i in range ( n1, n1 + n2 ): for j in range ( n1, n1 + n2 ): b[i,j] = r8bb_get ( n1, n2, ml, mu, a, i, j ) return b def r8bb_to_r8ge_test ( ): #*****************************************************************************80 # ## R8BB_TO_R8GE_TEST tests R8BB_TO_R8GE. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_TO_R8GE_TEST' ) print ( ' R8BB_TO_R8GE converts an R8BB matrix to R8GE format.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_indicator ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The R8BB matrix:' ) a_r8ge = r8bb_to_r8ge ( n1, n2, ml, mu, a ) print ( '' ) print ( ' The R8GE matrix:' ) print ( a_r8ge ) return def r8bb_zeros ( n1, n2, ml, mu ): #*****************************************************************************80 # ## R8BB_ZEROS zeros an R8BB matrix. # # Discussion: # # The R8BB storage format is for a 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. # # 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 (2*ML+MU+1)*N1 entries of A, using standard LINPACK # general band format. The reason for the factor of 2 in front of # ML is to allocate space that may be required if pivoting occurs. # # 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+ML+MU+1)+(J-1)*(2*ML+MU+1). # # Entries of A2: # # 1 <= I <= N1, N1+1 <= J <= N1+N2. # # Store the I, J entry into location # (2*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 # (2*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 # (2*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: # # 12 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, integer SEED, a seed for the random number generator. # # Output, real A((2*ML+MU+1)*N1+2*N1*N2+N2*N2), the R8BB matrix. # # Output, integer SEED, an updated seed for the random number generator. # import numpy as np nn = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 a = np.zeros ( nn ) return a def r8bb_zeros_test ( ): #*****************************************************************************80 # ## R8BB_ZEROS_TEST tests R8BB_ZEROS. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 12 July 2016 # # Author: # # John Burkardt # n1 = 3 n2 = 2 n = n1 + n2 ml = 1 mu = 1 na = ( 2 * ml + mu + 1 ) * n1 + 2 * n1 * n2 + n2 * n2 print ( '' ) print ( 'R8BB_ZEROS_TEST' ) print ( ' R8BB_ZEROS zeros an R8BB matrix.' ) print ( '' ) print ( ' Matrix order N = %d' % ( n ) ) print ( ' Matrix suborder N1 = %d' % ( n1 ) ) print ( ' Matrix suborder N2 = %d' % ( n2 ) ) print ( ' Lower bandwidth ML = %d' % ( ml ) ) print ( ' Upper bandwidth MU = %d' % ( mu ) ) # # Set the matrix. # a = r8bb_zeros ( n1, n2, ml, mu ) r8bb_print ( n1, n2, ml, mu, a, ' The zero R8BB matrix:' ) return def r8gb_fa ( n, ml, mu, a ): #*****************************************************************************80 # ## r8gb_fa() performs a LINPACK-style PLU factorization of a R8GB matrix. # # Discussion: # # An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU # is assumed to be entirely zero, except for the main diagonal, and # entries in the ML nearest subdiagonals, and MU nearest superdiagonals. # # LINPACK and LAPACK "R8GB" storage for such a matrix generally includes # room for ML extra superdiagonals, which may be required to store # nonzero entries generated during Gaussian elimination. # # 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 following program segment will set up the input. # # m = ml + 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 # # This uses rows ML+1 through 2*ML+MU+1 of the array A. # In addition, the first ML rows in the array are used for # elements generated during the triangularization. # The total number of rows needed in A is 2*ML+MU+1. # The ML+MU by ML+MU upper left triangle and the # ML by ML lower right triangle are not referenced. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 June 2016 # # Author: # # John Burkardt. # # Reference: # # Dongarra, Bunch, Moler, Stewart, # LINPACK User's Guide, # SIAM, 1979. # # 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(2*ML+MU+1,N), the matrix in band storage. The # columns of the matrix are stored in the columns of the array, # and the diagonals of the matrix are stored in rows ML+1 through # 2*ML+MU+1. # # Output, real ALU(2*ML+MU+1,N), the LU factors in band storage. # The L and U matrices are stored in a single array. # # Output, integer PIVOT(N), the pivot vector. # # Output, integer INFO, singularity flag. # 0, no singularity detected. # nonzero, the factorization failed on the INFO-th step. # import numpy as np alu = a.copy ( ) pivot = np.zeros ( n, dtype = np.int32 ) m = ml + mu + 1 info = 0 # # Zero out the initial fill-in columns. # j0 = mu + 2 j1 = min ( n, m ) - 1 for jz in range ( j0, j1 + 1 ): i0 = m + 1 - jz for i in range ( i0, ml + 1 ): alu[i-1,jz-1] = 0.0 jz = j1 ju = 0 for k in range ( 1, n ): # # Zero out the next fill-in column. # jz = jz + 1 if ( jz <= n ): for i in range ( 1, ml + 1 ): alu[i-1,jz-1] = 0.0 # # Find L = pivot index. # lm = min ( ml, n - k ) l = m for j in range ( m + 1, m + lm + 1 ): if ( abs ( alu[l-1,k-1] ) < abs ( alu[j-1,k-1] ) ): l = j pivot[k-1] = l + k - m # # Zero pivot implies this column already triangularized. # if ( alu[l-1,k-1] == 0.0 ): info = k print ( '' ) print ( 'R8GB_FA - Fatal error!' ) print ( ' Zero pivot on step %d' % ( info ) ) return alu, pivot, info # # Interchange if necessary. # t = alu[l-1,k-1] alu[l-1,k-1] = alu[m-1,k-1] alu[m-1,k-1] = t # # Compute multipliers. # for i in range ( m + 1, m + lm + 1 ): alu[i-1,k-1] = - alu[i-1,k-1] / alu[m-1,k-1] # # Row elimination with column indexing. # ju = max ( ju, mu + pivot[k-1] ) ju = min ( ju, n ) mm = m for j in range ( k + 1, ju + 1 ): l = l - 1 mm = mm - 1 if ( l != mm ): t = alu[l-1,j-1] alu[l-1,j-1] = alu[mm-1,j-1] alu[mm-1,j-1] = t for i in range ( 0, lm ): alu[mm+i,j-1] = alu[mm+i,j-1] + alu[mm-1,j-1] * alu[m+i,k-1] pivot[n-1] = n if ( alu[m-1,n-1] == 0.0 ): info = n print ( '' ) print ( 'R8GB_FA - Fatal error!' ) print ( ' Zero pivot on step %d' % ( info ) ) return alu, pivot, info def r8gb_sl ( n, ml, mu, a_lu, pivot, b, job ): #*****************************************************************************80 # ## r8gb_sl() solves a system factored by R8GB_FA. # # Discussion: # # An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU # is assumed to be entirely zero, except for the main diagonal, and # entries in the ML nearest subdiagonals, and MU nearest superdiagonals. # # LINPACK and LAPACK "R8GB" storage for such a matrix generally includes # room for ML extra superdiagonals, which may be required to store # nonzero entries generated during Gaussian elimination. # # 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. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt. # # Reference: # # Dongarra, Bunch, Moler, Stewart, # LINPACK User's Guide, # SIAM, 1979. # # 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(2*ML+MU+1,N), the LU factors from R8GB_FA. # # Input, integer PIVOT(N), the pivot vector from R8GB_FA. # # Input, real B(N), the right hand side vector. # # Input, integer JOB. # 0, solve A * x = b. # nonzero, solve A' * x = b. # # Output, real X(N), the solution. # x = b.copy ( ) m = mu + ml + 1 # # Solve A * x = b. # if ( job == 0 ): # # Solve L * Y = B. # if ( 1 <= ml ): for k in range ( 1, n ): lm = min ( ml, n - k ) l = pivot[k-1] if ( l != k ): t = x[l-1] x[l-1] = x[k-1] x[k-1] = t 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 L' * 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] l = pivot[k-1] if ( l != k ): t = x[l-1] x[l-1] = x[k-1] x[k-1] = t return x def r8gb_to_r8vec ( m, n, ml, mu, a ): #*****************************************************************************80 # ## r8gb_to_r8vec() copies a R8GB 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: # # 20 June 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(2*ML+MU+1,N), the array to be copied. # # Output, real X((2*ML+MU+1)*N), the vector. # import numpy as np x = np.zeros ( ( 2 * ml + mu + 1 ) * n ) for j in range ( 1, n + 1 ): ihi = min ( ml + mu, ml + mu + 1 - j ) for i in range ( 1, ihi + 1 ): x[i+(j-1)*(2*ml+mu+1)-1] = 0.0 ilo = max ( ihi + 1, 1 ) ihi = min ( 2 * ml + mu + 1, ml + mu + m + 1 - j ) for i in range ( ilo, ihi + 1 ): x[i+(j-1)*(2*ml+mu+1)-1] = a[i-1,j-1] ilo = ihi + 1 ihi = 2 * ml + mu + 1 for i in range ( ilo, ihi + 1 ): x[i+(j-1)*(2*ml+mu+1)-1] = 0.0 return x 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_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_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_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 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_r8gb ( m, n, ml, mu, x ): #*****************************************************************************80 # ## r8vec_to_r8gb() copies an R8VEC into a R8GB 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: # # 16 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((2*ML+MU+1)*N), the vector to be copied into the array. # # Output, real A(2*ML+MU+1,N), the array. # import numpy as np a = np.zeros ( [ 2 * ml + mu + 1, n ] ) for j in range ( 0, n ): for i in range ( 0, 2 * ml + mu + 1 ): if ( 0 <= i + j - ml - mu and i + j - ml - mu <= m - 1 ): a[i,j] = x[i+j*(2*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: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) r8bb_test ( ) timestamp ( )