#! /usr/bin/env python3 # def r8gb_test ( ): #*****************************************************************************80 # ## r8gb_test() tests r8gb(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 August 2022 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'r8gb_test():' ) print ( ' Python version: ' + platform.python_version ( ) ) print ( ' Test r8gb().' ) r8gb_det_test ( ) r8gb_dif2_test ( ) r8gb_fa_test ( ) r8gb_indicator_test ( ) r8gb_ml_test ( ) r8gb_mtv_test ( ) r8gb_mu_test ( ) r8gb_mv_test ( ) r8gb_nz_num_test ( ) r8gb_print_test ( ) r8gb_print_some_test ( ) r8gb_random_test ( ) r8gb_sl_test ( ) r8gb_to_r8ge_test ( ) r8gb_to_r8st_test ( ) # r8gb_to_r8sp_test ( ) r8gb_to_r8vec_test ( ) r8gb_trf_test ( ) r8gb_trs_test ( ) r8gb_zeros_test ( ) # r8ge_to_r8gb_test ( ) r8vec_to_r8gb_test ( ) # # Terminate. # print ( '' ) print ( 'r8gb_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 r8gb_det ( n, ml, mu, a_lu, pivot ): #*****************************************************************************80 # ## r8gb_det() computes the determinant of a matrix factored by r8gb_fa or R8GB_TRF. # # 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: # # Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum, # Hammarling, McKenney, Ostrouchov, Sorensen, # LAPACK User's Guide, # Second Edition, # SIAM, 1995. # # Dongarra, Bunch, Moler, Stewart, # LINPACK User's Guide, # SIAM, 1979. # # 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 or R8GB_TRF. # # Input, integer PIVOT(N), the pivot vector, as computed by r8gb_fa # or R8GB_TRF. # # Output, real DET, the determinant of the matrix. # det = 1.0 for j in range ( 1, n + 1 ): det = det * a_lu[ml+mu,j-1] for i in range ( 1, n + 1 ): if ( pivot[i-1] != i ): det = - det return det def r8gb_det_test ( ): #*****************************************************************************80 # ## r8gb_det_test() tests r8gb_det(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 August 2022 # # Author: # # John Burkardt # import numpy as np m = 10 n = m ml = 3 mu = 2 print ( '' ) print ( 'r8gb_det_test():' ) print ( ' r8gb_det() computes the determinant of an R8GB matrix' ) print ( ' which has been factored by r8gb_fa.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' A random R8GB matrix:' ) # # Copy the matrix into a general array. # a2 = r8gb_to_r8ge ( m, n, ml, mu, a ) # # Factor the matrix. # a_lu, pivot, info = r8gb_fa ( n, ml, mu, a ) # # Compute the determinant. # det = r8gb_det ( n, ml, mu, a_lu, pivot ) print ( '' ) print ( ' r8gb_det() computes the determinant = ', det ) # # Recompute the determinant, using the R8GE matrix. # determ = np.linalg.det ( a2 ) print ( ' np.linalg.det() computes the determinant = ', determ ) return def r8gb_dif2 ( m, n, ml, mu ): #*****************************************************************************80 # ## r8gb_dif2() sets up an R8GB second difference matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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 # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Output, real A(2*ML+MU+1,N), the R8GB matrix. # import numpy as np a = np.zeros ( [ 2 * ml + mu + 1, n ] ) for j in range ( 1, n + 1 ): for diag in range ( 1, 2 * ml + mu + 2 ): i = diag + j - ml - mu - 1 if ( i == j ): a[diag-1,j-1] = 2.0 elif ( i == j - 1 or i == j + 1 ): a[diag-1,j-1] = -1.0 return a def r8gb_dif2_test ( ): #*****************************************************************************80 # ## r8gb_dif2_test() tests r8gb_dif2(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 5 ml = 1 mu = 1 print ( '' ) print ( 'r8gb_dif2_test():' ) print ( ' r8gb_dif2() returns an R8GB second difference matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_dif2 ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB second difference 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 ', info ) raise Exception ( 'r8gb_fa(): Fatal error!' ) # # 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 ', info ) raise Exception ( 'r8gb_fa(): Fatal error!' ) return alu, pivot, info def r8gb_fa_test ( ): #*****************************************************************************80 # ## r8gb_fa_test() tests r8gb_fa(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 June 2016 # # Author: # # John Burkardt # m = 5 n = m ml = 1 mu = 2 print ( '' ) print ( 'r8gb_fa_test():' ) print ( ' r8gb_fa() computes the PLU factors of an R8GB matrix.' ) print ( '' ) print ( ' Number of matrix rows M = ', m ) print ( ' Number of matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The banded matrix:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8gb_mv ( m, n, ml, mu, a, x ) # # Factor the matrix. # alu, pivot, info = r8gb_fa ( n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'r8gb_fa_test - Fatal error!' ) print ( ' r8gb_fa declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gb_fa_test(): Fatal error!' ) # # Solve the linear system. # job = 0 x = r8gb_sl ( n, ml, mu, alu, pivot, b, job ) r8vec_print ( n, x, ' Solution to A*x=b:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # job = 1 b = r8gb_ml ( n, ml, mu, alu, pivot, x, job ) r8vec_print ( n, b, ' Right hand side of transposed system:' ) # # Solve the linear system. # job = 1 x = r8gb_sl ( n, ml, mu, alu, pivot, b, job ) r8vec_print ( n, x, ' Solution to A''x=b:' ) return def r8gb_indicator ( m, n, ml, mu ): #*****************************************************************************80 # ## r8gb_indicator() sets up a R8GB indicator matrix. # # Discussion: # # Note that the R8GB storage format includes extra room for # fillin entries that occur during Gauss elimination. This routine # will supply zero values for those entries. # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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 # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Output, real A(2*ML+MU+1,N), the R8GB matrix. # import numpy as np a = np.zeros ( [ 2 * ml + mu + 1, n ] ) fac = 10 ** ( i4_log_10 ( n ) + 1 ) k = 0 for j in range ( 1, n + 1 ): for diag in range ( 1, 2 * ml + mu + 2 ): i = diag + j - ml - mu - 1 if ( 1 <= i and i <= m and i - ml <= j and j <= i + mu ): a[diag-1,j-1] = float ( fac * i + j ) elif ( 1 <= i and i <= m and i - ml <= j and j <= i + mu + ml ): value = 0.0 else: k = k + 1 a[diag-1,j-1] = - float ( k ) return a def r8gb_indicator_test ( ): #*****************************************************************************80 # ## r8gb_indicator_test() tests r8gb_indicator(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 10 n = 8 ml = 3 mu = 2 print ( '' ) print ( 'r8gb_indicator_test():' ) print ( ' r8gb_indicator() returns an R8GB indicator matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB indicator matrix:' ) return def r8gb_ml ( n, ml, mu, a_lu, pivot, x, job ): #*****************************************************************************80 # ## R8GB_ML computes A * x or A' * X, using r8gb_fa factors. # # 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. # # It is assumed that r8gb_fa has overwritten the original matrix # information by LU factors. R8GB_ML is able to reconstruct the # original matrix from the LU factor data. # # R8GB_ML allows the user to check that the solution of a linear # system is correct, without having to save an unfactored copy # of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 June 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(2*ML+MU+1,N), the LU factors from r8gb_fa. # # Input, integer PIVOT(N), the pivot vector computed by r8gb_fa. # # Input, real X(N), the vector to be multiplied. # # Input, integer JOB, specifies the operation to be done: # JOB = 0, compute A * x. # JOB nonzero, compute A' * X. # # Output, real B(N), the result of the multiplication. # b = x.copy ( ) if ( job == 0 ): # # Y = U * X. # for j in range ( 1, n + 1 ): ilo = max ( 1, j - ml - mu ) for i in range ( ilo, j ): b[i-1] = b[i-1] + a_lu[i-j+ml+mu,j-1] * b[j-1] b[j-1] = a_lu[j-j+ml+mu,j-1] * b[j-1] # # B = PL * Y = PL * U * X = A * x. # for j in range ( n - 1, 0, -1 ): ihi = min ( n, j + ml ) for i in range ( j + 1, ihi + 1 ): b[i-1] = b[i-1] - a_lu[i-j+ml+mu,j-1] * b[j-1] k = pivot[j-1] if ( k != j ): temp = b[k-1] b[k-1] = b[j-1] b[j-1] = temp else: # # Y = ( PL )' * X. # for j in range ( 1, n ): k = pivot[j-1] if ( k != j ): temp = b[k-1] b[k-1] = b[j-1] b[j-1] = temp jhi = min ( n, j + ml ) for i in range ( j + 1, jhi + 1 ): b[j-1] = b[j-1] - b[i-1] * a_lu[i-j+ml+mu,j-1] # # B = U' * Y = ( PL * U )' * X = A' * X. # for i in range ( n, 0, -1 ): jhi = min ( n, i + ml + mu ) for j in range ( i + 1, jhi + 1 ): b[j-1] = b[j-1] + b[i-1] * a_lu[i-j+ml+mu,j-1] b[i-1] = b[i-1] * a_lu[i-i+ml+mu,i-1] return b def r8gb_ml_test ( ): #*****************************************************************************80 # ## R8GB_ML_test tests R8GB_ML. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 10 n = m ml = 1 mu = 2 print ( '' ) print ( 'R8GB_ML_test' ) print ( ' R8GB_ML computes A*x or A''*X' ) print ( ' where A has been factored by r8gb_fa().' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) for job in range ( 0, 2 ): # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # if ( job == 0 ): b = r8gb_mv ( m, n, ml, mu, a, x ) else: b = r8gb_mtv ( m, n, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8gb_fa ( n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8GB_ML_test - Fatal error!' ) print ( ' r8gb_fa declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gm_ml(): Fatal error!' ) # # Now multiply factored matrix times solution to get right hand side again. # b2 = r8gb_ml ( n, ml, mu, a_lu, pivot, x, job ) if ( job == 0 ): r8vec2_print ( b, b2, ' A*x and PLU*x' ) else: r8vec2_print ( b, b2, ' A''*x and (PLU)''*x' ) return def r8gb_mtv ( m, n, ml, mu, a, x ): #*****************************************************************************80 # ## R8GB_MTV multiplies a vector by 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. # # LINPACK and LAPACK storage of general band matrices requires # an extra ML upper diagonals for possible fill in entries during # Gauss elimination. This routine does not access any entries # in the fill in diagonals, because it assumes that the matrix # has NOT had Gauss elimination applied to it. If the matrix # has been Gauss eliminated, then the routine R8GB_MU must be # used instead. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Input, real X(M), the vector to be multiplied by A. # # Output, real B(N), the product X*A. # import numpy as np b = np.zeros ( n ) for j in range ( 1, n + 1 ): ilo = max ( 1, j - mu ) ihi = min ( m, j + ml ) for i in range ( ilo, ihi + 1 ): b[j-1] = b[j-1] + x[i-1] * a[i-j+ml+mu,j-1] return b def r8gb_mtv_test ( ): #*****************************************************************************80 # ## R8GB_MTV_test tests R8GB_MTV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 5 ml = 1 mu = 2 print ( '' ) print ( 'R8GB_MTV_test' ) print ( ' R8GB_MTV computes b=A''*x, where A is an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The random R8GB matrix:' ) x = r8vec_indicator1 ( m ) r8vec_print ( m, x, ' x:' ) b = r8gb_mtv ( m, n, ml, mu, a, x ) r8vec_print ( n, b, ' b=A''*x:' ) return def r8gb_mu ( n, ml, mu, a_lu, pivot, x, job ): #*****************************************************************************80 # ## R8GB_MU computes A * x or A' * X, using R8GB_TRF factors. # # Warning: # # This routine must be updated to allow for rectangular matrices. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # It is assumed that R8GB_TRF has overwritten the original matrix # information by LU factors. R8GB_MU is able to reconstruct the # original matrix from the LU factor data. # # R8GB_MU allows the user to check that the solution of a linear # system is correct, without having to save an unfactored copy # of the matrix. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Reference: # # Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum, # Hammarling, McKenney, Ostrouchov, Sorensen, # LAPACK User's Guide, # Second Edition, # SIAM, 1995. # # 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_TRF. # # Input, integer PIVOT(N), the pivot vector computed by R8GB_TRF. # # Input, real X(N), the vector to be multiplied. # # Input, integer JOB, specifies the operation to be done: # JOB = 0, compute A * x. # JOB nonzero, compute A' * X. # # Output, real B(N), the result of the multiplication. # b = x.copy ( ) if ( job == 0 ): # # Y = U * X. # for j in range ( 1, n + 1 ): ilo = max ( 1, j - ml - mu ) for i in range ( ilo, j ): b[i-1] = b[i-1] + a_lu[i-j+ml+mu,j-1] * b[j-1] b[j-1] = a_lu[j-j+ml+mu,j-1] * b[j-1] # # B = PL * Y = PL * U * X = A * x. # for j in range ( n - 1, 0, -1 ): ihi = min ( n, j + ml ) for i in range ( j + 1, ihi + 1 ): b[i-1] = b[i-1] + a_lu[i-j+ml+mu,j-1] * b[j-1] k = pivot[j-1] if ( k != j ): t = b[k-1] b[k-1] = b[j-1] b[j-1] = t else: # # Y = ( PL )' * X. # for j in range ( 1, n ): k = pivot[j-1] if ( k != j ): t = b[k-1] b[k-1] = b[j-1] b[j-1] = t jhi = min ( n, j + ml ) for i in range ( j + 1, jhi + 1 ): b[j-1] = b[j-1] + b[i-1] * a_lu[i-j+ml+mu,j-1] # # B = U' * Y = ( PL * U )' * X = A' * X. # for i in range ( n, 0, -1 ): jhi = min ( n, i + ml + mu ) for j in range ( i + 1, jhi + 1 ): b[j-1] = b[j-1] + b[i-1] * a_lu[i-j+ml+mu,j-1] b[i-1] = b[i-1] * a_lu[i-i+ml+mu,i-1] return b def r8gb_mu_test ( ): #*****************************************************************************80 # ## R8GB_MU_test tests R8GB_MU. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 10 n = m ml = 1 mu = 2 print ( '' ) print ( 'R8GB_MU_test' ) print ( ' R8GB_MU computes A*x or A''*X' ) print ( ' where A has been factored by R8GB_TRF.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) for job in range ( 0, 2 ): # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # if ( job == 0 ): b = r8gb_mv ( m, n, ml, mu, a, x ) else: b = r8gb_mtv ( m, n, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8gb_trf ( m, n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8GB_MU_test - Fatal error!' ) print ( ' R8GB_TRF declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gb_mu_test(): Fatal error!' ) # # Now multiply factored matrix times solution to get right hand side again. # b2 = r8gb_mu ( n, ml, mu, a_lu, pivot, x, job ) if ( job == 0 ): r8vec2_print ( b, b2, ' A*x and PLU*x' ) else: r8vec2_print ( b, b2, ' A''*x and (PLU)''*x' ) return def r8gb_mv ( m, n, ml, mu, a, x ): #*****************************************************************************80 # ## R8GB_MV multiplies a R8GB matrix times a vector. # # 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. # # LINPACK and LAPACK storage of general band matrices requires # an extra ML upper diagonals for possible fill in entries during # Gauss elimination. This routine does not access any entries # in the fill in diagonals, because it assumes that the matrix # has NOT had Gauss elimination applied to it. If the matrix # has been Gauss eliminated, then the routine R8GB_MU must be # used instead. # # 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 M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Input, real X(N), the vector to be multiplied by A. # # Output, real B(M), the product A * x. # import numpy as np b = np.zeros ( m ) for i in range ( 1, m + 1 ): jlo = max ( 1, i - ml ) jhi = min ( n, i + mu ) for j in range ( jlo, jhi + 1 ): b[i-1] = b[i-1] + a[i-j+ml+mu,j-1] * x[j-1] return b def r8gb_mv_test ( ): #*****************************************************************************80 # ## R8GB_MV_test tests R8GB_MV. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 5 ml = 1 mu = 2 print ( '' ) print ( 'R8GB_MV_test' ) print ( ' R8GB_MV computes b=A*x, where A is an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The random R8GB matrix:' ) x = r8vec_indicator1 ( m ) r8vec_print ( m, x, ' x:' ) b = r8gb_mv ( m, n, ml, mu, a, x ) r8vec_print ( n, b, ' b=A*x:' ) return def r8gb_nz_num ( m, n, ml, mu, a ): #*****************************************************************************80 # ## R8GB_NZ_NUM counts the nonzeros in a R8GB matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine will examine # values it finds in these extra bands, so that both unfactored # and factored matrices can be handled. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Output, integer NZ_NUM, the number of nonzero entries in A. # nz_num = 0 for i in range ( 1, m + 1 ): jlo = max ( 1, i - ml ) jhi = min ( n, i + mu + ml ) for j in range ( jlo, jhi + 1 ): if ( a[i-j+ml+mu,j-1] != 0.0 ): nz_num = nz_num + 1 return nz_num def r8gb_nz_num_test ( ): #*****************************************************************************80 # ## R8GB_NZ_NUM_test tests R8GB_NZ_NUM. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 10 n = m ml = 1 mu = 2 print ( '' ) print ( 'R8GB_NZ_NUM_test' ) print ( ' R8GB_NZ_NUM counts the nonzero entries in an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) # # Make some zero entries. # for j in range ( 1, n + 1 ): for diag in range ( 1, 2 * ml + mu + 2 ): if ( a[diag-1,j-1] < 0.3 ): a[diag-1,j-1] = 0.0 r8gb_print ( m, n, ml, mu, a, ' The R8GB matrix:' ) nz_num = r8gb_nz_num ( m, n, ml, mu, a ) print ( '' ) print ( ' Nonzero entries = ', nz_num ) return def r8gb_print ( m, n, ml, mu, a, title ): #*****************************************************************************80 # ## R8GB_PRINT prints a banded 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. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1.. # # Input, real A(2*ML+MU+1,N), the M by N band matrix, stored in LINPACK # or LAPACK general band storage mode. # # Input, string TITLE, a title to be printed. # r8gb_print_some ( m, n, ml, mu, a, 0, 0, m - 1, n - 1, title ) return def r8gb_print_test ( ): #*****************************************************************************80 # ## R8GB_PRINT_test tests R8GB_PRINT. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 March 2009 # # Author: # # John Burkardt # m = 8 n = 10 ml = 1 mu = 3 print ( '' ) print ( 'R8GB_PRINT_test' ) print ( ' R8GB_PRINT prints an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB matrix:' ) return def r8gb_print_some ( m, n, ml, mu, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## R8GB_PRINT_SOME prints some of a banded 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. # # Only entries in rows ILO to IHI, columns JLO to JHI are considered. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1.. # # Input, real A(2*ML+MU+1,N), the M by N band matrix, stored in LINPACK # or LAPACK general band storage mode. # # Input, integer ILO, JLO, IHI, JHI, designate the first row and # column, and the last row and column to be printed. # # Input, string TITLE, a title. # print ( '' ) print ( title ) incx = 5 for j2lo in range ( max ( jlo, 0 ), min ( jhi + 1, n ), incx ): j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n - 1 ) j2hi = min ( j2hi, jhi ) print ( '' ) print ( ' Col: ', end = '' ) for j in range ( j2lo, j2hi + 1 ): print ( '%7d ' % ( j ), end = '' ) print ( '' ) print ( ' Row' ) i2lo = max ( ilo, 0 ) i2lo = max ( i2lo, j2lo - mu ) i2hi = min ( ihi, m - 1 ) i2hi = min ( i2hi, j2hi + ml ) for i in range ( i2lo, i2hi + 1 ): print ( '%7d :' % ( i ), end = '' ) for j in range ( j2lo, j2hi + 1 ): if ( mu < j - i or ml < i - j ): print ( ' ', end = '' ) else: print ( '%14g' % ( a[i-j+ml+mu,j] ), end = '' ) print ( '' ) return def r8gb_print_some_test ( ): #*****************************************************************************80 # ## R8GB_PRINT_SOME_test tests R8GB_PRINT_SOME. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 8 n = 10 ml = 1 mu = 3 print ( '' ) print ( 'R8GB_PRINT_SOME_test' ) print ( ' R8GB_PRINT_SOME prints some of an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print_some ( m, n, ml, mu, a, 4, 3, 6, 9, ' Rows(4-6), Cols (3-9)' ) return def r8gb_random ( m, n, ml, mu ): #*****************************************************************************80 # ## r8gb_random() randomizes 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine assumes it is setting # up an unfactored matrix, so it only uses the first MU upper bands, # and does not place nonzero values in the fillin bands. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Output, real A(2*ML+MU+1,N), the R8GB matrix. # from numpy.random import default_rng import numpy as np rng = default_rng ( ) a = np.zeros ( [ 2 * ml + mu + 1, n ] ) for j in range ( 1, n + 1 ): for irow in range ( 1, 2 * ml + mu + 2 ): i = irow + j - ml - mu - 1 if ( ml < irow and 1 <= i and i <= m ): r = rng.random ( size = 1 ) a[irow-1,j-1] = r return a def r8gb_random_test ( ): #*****************************************************************************80 # ## R8GB_RANDOM_test tests R8GB_RANDOM. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 5 ml = 1 mu = 2 print ( '' ) print ( 'R8GB_RANDOM_test' ) print ( ' R8GB_RANDOM returns a random R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The random R8GB matrix:' ) return 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_sl_test ( ): #*****************************************************************************80 # ## R8GB_SL_test tests R8GB_SL. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = m ml = 1 mu = 2 print ( '' ) print ( 'R8GB_SL_test' ) print ( ' R8GB_SL solves a linear system factored by r8gb_fa.' ) print ( '' ) print ( ' Number of matrix rows M = ', m ) print ( ' Number of matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The banded matrix:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8gb_mv ( m, n, ml, mu, a, x ) # # Factor the matrix. # alu, pivot, info = r8gb_fa ( n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8GB_SL_test - Fatal error!' ) print ( ' r8gb_fa declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gb_sl_test(): Fatal error!' ) # # Solve the linear system. # job = 0 x = r8gb_sl ( n, ml, mu, alu, pivot, b, job ) r8vec_print ( n, x, ' Solution to A*x=b:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # job = 1 b = r8gb_ml ( n, ml, mu, alu, pivot, x, job ) r8vec_print ( n, b, ' Right hand side of transposed system:' ) # # Solve the linear system. # job = 1 x = r8gb_sl ( n, ml, mu, alu, pivot, b, job ) r8vec_print ( n, x, ' Solution to A''x=b:' ) return def r8gb_to_r8ge ( m, n, ml, mu, a ): #*****************************************************************************80 # ## R8GB_TO_R8GE copies a R8GB matrix to a R8GE matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine will copy nonzero # values it finds in these extra bands, so that both unfactored # and factored matrices can be handled. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrices. # M must be positive. # # Input, integer N, the number of columns of the matrices. # N must be positive. # # Input, integer ML, MU, the lower and upper bandwidths of A1. # ML and MU must be nonnegative, and no greater than min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Output, real B(M,N), the R8GE matrix. # import numpy as np b = np.zeros ( [ m, n ] ) for i in range ( 1, m + 1 ): for j in range ( 1, n + 1 ): if ( i - ml <= j and j <= i + mu ): b[i-1,j-1] = a[ml+mu+i-j,j-1] return b def r8gb_to_r8ge_test ( ): #*****************************************************************************80 # ## R8GB_TO_R8GE_test tests R8GB_TO_R8GE. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 8 ml = 2 mu = 1 print ( '' ) print ( 'R8GB_TO_R8GE_test' ) print ( ' R8GB_TO_R8GE copies an R8GB matrix to an R8GE matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB matrix:' ) a_r8ge = r8gb_to_r8ge ( m, n, ml, mu, a ) print ( '' ) print ( ' The R8GE matrix:' ) print ( a_r8ge ) return def r8gb_to_r8st ( m, n, ml, mu, a, nz_num ): #*****************************************************************************80 # ## r8gb_to_r8st() copies a R8GB matrix to a R8ST matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine will copy nonzero # values it finds in these extra bands, so that both unfactored # and factored matrices can be handled. # # The r8st storage format corresponds to the SLAP Triad format. # # The r8st storage format stores the row, column and value of each nonzero # entry of a sparse matrix. The entries may be given in any order. No # check is made for the erroneous case in which a given matrix entry is # specified more than once. # # There is a symmetry option for square matrices. If the symmetric storage # option is used, the format specifies that only nonzeroes on the diagonal # and lower triangle are stored. However, this routine makes no attempt # to enforce this. The only thing it does is to "reflect" any nonzero # offdiagonal value. Moreover, no check is made for the erroneous case # in which both A(I,J) and A(J,I) are specified, but with different values. # # This routine reorders the entries of A so that the first N entries # are exactly the diagonal entries of the matrix, in order. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrices. # M must be positive. # # Input, integer N, the number of columns of the matrices. # N must be positive. # # Input, integer ML, MU, the lower and upper bandwidths of A1. # ML and MU must be nonnegative, and no greater than min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Input, integer NZ_NUM, the number of nonzero entries in A. # # Output, integer SYM, is 0 if the matrix is not symmetric, and 1 # if the matrix is symmetric. If the matrix is symmetric, then # only the nonzeroes on the diagonal and in the lower triangle are stored. # For this routine, SYM is always output 0. # # Output, integer ROW(NZ_NUM), the row indices. # # Output, integer COL(NZ_NUM), the column indices. # # Output, real B(NZ_NUM), the r8st matrix. # import numpy as np row = np.zeros ( nz_num, dtype = int ) col = np.zeros ( nz_num, dtype = int ) b = np.zeros ( nz_num, dtype = float ) sym = 0 nz = 0 for i in range ( 1, m + 1 ): for j in range ( 1, n + 1 ): if ( i - ml <= j and j <= i + mu + ml ): if ( a[ml+mu+i-j,j-1] != 0.0 ): if ( nz_num <= nz ): print ( '' ) print ( 'R8GB_TO_r8st - Fatal error!' ) print ( ' NZ_NUM = ', nz_num ) print ( ' But the matrix has more nonzeros than that!' ) raise Exception ( 'R8GB_TO_r8st - Fatal error!' ) row[nz] = i col[nz] = j b[nz] = a[ml+mu+i-j,j-1] nz = nz + 1 if ( nz < nz_num ): print ( '' ) print ( 'R8GB_TO_r8st - Warning!' ) print ( ' NZ_NUM = ', nz_num ) print ( ' But the number of nonzeros is ', nz ) return sym, row, col, b def r8gb_to_r8st_test ( ): #*****************************************************************************80 # ## r8gb_to_r8st_test() tests r8gb_to_r8st(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 22 August 2022 # # Author: # # John Burkardt # m = 5 n = 8 ml = 2 mu = 1 print ( '' ) print ( 'R8GB_TO_r8st_test' ) print ( ' R8GB_TO_r8st copies a R8GB matrix to a r8st matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB matrix:' ) nz_num = r8gb_nz_num ( m, n, ml, mu, a ) print ( ' Nonzeros NZ_NUM = ', nz_num ) sym, row, col, b = r8gb_to_r8st ( m, n, ml, mu, a, nz_num ) r8st_print ( m, n, nz_num, sym, row, col, b, ' The r8st matrix:' ) return def r8gb_to_r8sp ( m, n, ml, mu, a, nz_num ): #*****************************************************************************80 # ## r8gb_to_r8sp() copies a R8GB matrix to a R8SP matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine will copy nonzero # values it finds in these extra bands, so that both unfactored # and factored matrices can be handled. # # The R8SP storage format stores the row, column and value of each nonzero # entry of a sparse matrix. # # It is possible that a pair of indices (I,J) may occur more than # once. Presumably, in this case, the intent is that the actual value # of A(I,J) is the sum of all such entries. This is not a good thing # to do, but I seem to have come across this in MATLAB. # # The R8SP format is used by CSPARSE ("sparse triplet"), SLAP # ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrices. # M must be positive. # # Input, integer N, the number of columns of the matrices. # N must be positive. # # Input, integer ML, MU, the lower and upper bandwidths of A1. # ML and MU must be nonnegative, and no greater than min(M,N)-1. # # Input, real A(2*ML+MU+1,N), the R8GB matrix. # # Input, integer NZ_NUM, the number of nonzero entries in A. # This number can be obtained by calling R8GB_NZ_NUM. # # Output, integer ROW(NZ_NUM), the row indices. # # Output, integer COL(NZ_NUM), the column indices. # # Output, real B(NZ_NUM), the R8SP matrix. # import numpy as np b = np.zeros ( nz_num ) col = np.zeros ( nz_num, dtype = np.int32 ) row = np.zeros ( nz_num, dtype = np.int32 ) nz = 0 for i in range ( 1, m + 1 ): jlo = max ( 1, i - ml ) jhi = min ( n, i + mu + ml ) for j in range ( jlo, jhi + 1 ): if ( a[ml+mu+i-j,j-1] == 0.0 ): continue if ( nz_num <= nz ): print ( '' ) print ( 'R8GB_TO_R8SP - Fatal error!' ) print ( ' NZ_NUM = ', nz_num ) print ( ' But the matrix has more nonzeros than that!' ) raise Exception ( 'R8GB_TO_DS3 - Fatal error!' ) row[nz] = i col[nz] = j b[nz] = a[ml+mu+i-j,j-1] nz = nz + 1 if ( nz < nz_num ): print ( '' ) print ( 'R8GB_TO_R8SP - Warning!' ) print ( ' NZ_NUM = ', nz_num ) print ( ' But the number of nonzeros is ', nz ) return row, col, b def r8gb_to_r8sp_test ( ): #*****************************************************************************80 # ## R8GB_TO_R8SP_test tests R8GB_TO_R8SP. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 8 ml = 2 mu = 1 print ( '' ) print ( 'R8GB_TO_R8SP_test' ) print ( ' R8GB_TO_R8SP copies an R8GB matrix to an R8SP matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB matrix:' ) nz_num = r8gb_nz_num ( m, n, ml, mu, a ) print ( ' Nonzeros NZ_NUM = ', nz_num ) row, col, a_r8sp = r8gb_to_r8sp ( m, n, ml, mu, a, nz_num ) r8sp_print ( m, n, nz_num, row, col, a_r8sp, ' The R8SP matrix:' ) return 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 r8gb_to_r8vec_test ( ): #*****************************************************************************80 # ## R8GB_TO_R8VEC_test tests R8GB_TO_R8VEC. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 14 March 2009 # # Author: # # John Burkardt # m = 5 n = 8 ml = 2 mu = 1 print ( '' ) print ( 'R8GB_TO_R8VEC_test' ) print ( ' R8GB_TO_R8VEC converts an R8GB matrix to an R8VEC.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB indicator matrix:' ) x = r8gb_to_r8vec ( m, n, ml, mu, a ) k = 0 for j in range ( 1, n + 1 ): for i in range ( 1, 2 * ml + mu + 2 ): k = k + 1 print ( '%4d %4d %4d %14f' % ( i, j, k, x[k-1] ) ) return def r8gb_trf ( m, n, ml, mu, a ): #*****************************************************************************80 # ## R8GB_TRF performs a LAPACK-style PLU factorization of a R8GB matrix. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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. # # This is a simplified, standalone version of the LAPACK # routine R8GBTRF. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt. # # Reference: # # Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum, # Hammarling, McKenney, Ostrouchov, Sorensen, # LAPACK User's Guide, # Second Edition, # SIAM, 1995. # # Parameters: # # Input, integer M, the number of rows of the matrix A. 0 <= M. # # Input, integer N, the number of columns of the matrix A. 0 <= N. # # Input, integer ML, the number of subdiagonals within the band of A. # 0 <= ML. # # Input, integer MU, the number of superdiagonals within the band of A. # 0 <= MU. # # Input, real A(2*ML+MU+1,N), the matrix A in band storage. # # Output, real A_LU(2*ML+MU+1,N), information about the PLU factorization. # # Output, integer PIVOT(min(M,N)), the pivot indices # for 1 <= i <= min(M,N), row i of the matrix was interchanged with # row IPIV(i). # # Output, integer INFO, error flag. # = 0: successful exit # < 0: an input argument was illegal # > 0: if INFO = +i, U(i,i) is exactly zero. The factorization # has been completed, but the factor U is exactly # singular, and division by zero will occur if it is used # to solve a system of equations. # import numpy as np info = 0 a_lu = a.copy ( ) pivot = np.zeros ( n, dtype = np.int32 ) # # KV is the number of superdiagonals in the factor U, allowing for fill-in. # kv = mu + ml # # Set fill-in elements in columns MU+2 to KV to zero. # for j in range ( mu + 2, min ( kv, n ) + 1 ): for i in range ( kv - j + 2, ml + 1 ): a_lu[i-1,j-1] = 0.0 # # JU is the index of the last column affected by the current stage # of the factorization. # ju = 1 for j in range ( 1, min ( m, n ) + 1 ): # # Set the fill-in elements in column J+KV to zero. # if ( j + kv <= n ): for i in range ( 1, ml + 1 ): a_lu[i-1,j+kv-1] = 0.0 # # Find the pivot and test for singularity. # KM is the number of subdiagonal elements in the current column. # km = min ( ml, m - j ) piv = abs ( a_lu[kv,j-1] ) jp = kv + 1 for i in range ( kv + 2, kv + km + 2 ): if ( piv < abs ( a_lu[i-1,j-1] ) ): piv = abs ( a_lu[i-1,j-1] ) jp = i jp = jp - kv pivot[j-1] = jp + j - 1 if ( a_lu[kv+jp-1,j-1] != 0.0 ): ju = max ( ju, min ( j + mu + jp - 1, n ) ) # # Apply interchange to columns J to JU. # if ( jp != 1 ): for i in range ( 0, ju - j + 1 ): t = a_lu[kv+jp-i-1,j+i-1] a_lu[kv+jp-i-1,j+i-1] = a_lu[kv-i,j+i-1] a_lu[kv-i,j+i-1] = t # # Compute the multipliers. # if ( 0 < km ): for i in range ( kv + 2, kv + km + 2 ): a_lu[i-1,j-1] = a_lu[i-1,j-1] / a_lu[kv,j-1] # # Update the trailing submatrix within the band. # if ( j < ju ): for k in range ( 1, ju - j + 1 ): if ( a_lu[kv-k,j+k-1] != 0.0 ): for i in range ( 1, km + 1 ): a_lu[kv+i-k,j+k-1] = a_lu[kv+i-k,j+k-1] \ - a_lu[kv+i,j-1] * a_lu[kv-k,j+k-1] else: # # If pivot is zero, set INFO to the index of the pivot # unless a zero pivot has already been found. # if ( info == 0 ): info = j return a_lu, pivot, info def r8gb_trf_test ( ): #*****************************************************************************80 # ## R8GB_TRF_test tests R8GB_TRF. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # n = 10 m = n ml = 1 mu = 2 nrhs = 1 print ( '' ) print ( 'R8GB_TRF_test' ) print ( ' R8GB_TRF computes the PLU factors of an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8gb_mv ( m, n, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8gb_trf ( m, n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8GB_TRF_test - Fatal error!' ) print ( ' R8GB_TRF declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gb_trf_test(): Fatal error!' ) # # Solve the linear system. # Note that, because of quirks in MATLAB, we need to copy our vector-based # data to and from 2D arrays. # b_mat = r8vec_to_r8ge ( n, nrhs, b ) x_mat, info = r8gb_trs ( n, ml, mu, nrhs, 'N', a_lu, pivot, b_mat ) x = r8ge_to_r8vec ( n, nrhs, x_mat ) r8vec_print ( n, x, ' Solution:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # job = 1 b = r8gb_mu ( n, ml, mu, a_lu, pivot, x, job ) # # Solve the linear system. # b_mat = r8vec_to_r8ge ( n, nrhs, b ) x_mat, info = r8gb_trs ( n, ml, mu, nrhs, 'T', a_lu, pivot, b_mat ) x = r8ge_to_r8vec ( n, nrhs, x_mat ) r8vec_print ( n, x, ' Solution to transposed system:' ) return def r8gb_trs ( n, ml, mu, nrhs, trans, a_lu, pivot, b ): #*****************************************************************************80 # ## R8GB_TRS solves a linear system factored by R8GB_TRF. # # Discussion: # # The R8GB storage format is for an M by N banded matrix, with lower # bandwidth ML and upper bandwidth MU. Storage includes room for ML # extra superdiagonals, which may be required to store nonzero entries # generated during Gaussian elimination. # # 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: # # Anderson, Bai, Bischof, Demmel, Dongarra, Du Croz, Greenbaum, # Hammarling, McKenney, Ostrouchov, Sorensen, # LAPACK User's Guide, # Second Edition, # SIAM, 1995. # # Parameters: # # Input, integer N, the order of the matrix A. # N must be positive. # # Input, integer ML, the number of subdiagonals within the band of A. # ML must be at least 0, and no greater than N - 1. # # Input, integer MU, the number of superdiagonals within the band of A. # MU must be at least 0, and no greater than N - 1. # # Input, integer NRHS, the number of right hand sides and the number of # columns of the matrix B. NRHS must be positive. # # Input, character TRANS, specifies the form of the system. # 'N': A * x = b (No transpose) # 'T': A'* X = B (Transpose) # 'C': A'* X = B (Conjugate transpose = Transpose) # # Input, real A_LU(2*ML+MU+1,N), the LU factors from R8GB_TRF. # # Input, integer PIVOT(N), the pivot indices for 1 <= I <= N, row I # of the matrix was interchanged with row PIVOT(I). # # Input, real B(N,NRHS), the right hand side vectors. # # Output, real X(N,NRHS), the solution vectors, X. # # Output, integer INFO, error flag. # = 0: successful exit # < 0: if INFO = -K, the K-th argument had an illegal value # # # Test the input parameters. # info = 0 if ( trans != 'N' and trans != 'n' and \ trans != 'T' and trans != 't' and \ trans != 'C' and trans != 'c' ): info = -1 return elif ( n <= 0 ): info = -2 return elif ( ml < 0 ): info = -3 return elif ( mu < 0 ): info = -4 return elif ( nrhs <= 0 ): info = -5 return x = b.copy ( ) kd = mu + ml + 1 # # Solve A * x = b. # # Solve L * x = b, overwriting b with x. # # L is represented as a product of permutations and unit lower # triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1), # where each transformation L(i) is a rank-one modification of # the identity matrix. # if ( trans == 'N' or trans == 'n' ): if ( 0 < ml ): for j in range ( 1, n ): lm = min ( ml, n - j ) l = pivot[j-1] for k in range ( 1, nrhs + 1 ): t = x[l-1,k-1] x[l-1,k-1] = x[j-1,k-1] x[j-1,k-1] = t for k in range ( 1, nrhs + 1 ): if ( x[j-1,k-1] != 0.0 ): for i in range ( 0, lm ): x[j+i,k-1] = x[j+i,k-1] - a_lu[kd+i,j-1] * x[j-1,k-1] # # Solve U * x = b, overwriting b with x. # for k in range ( 0, nrhs ): for j in range ( n, 0, -1 ): if ( x[j-1,k] != 0.0 ): l = ml + mu + 1 - j x[j-1,k] = x[j-1,k] / a_lu[ml+mu,j-1] for i in range ( j - 1, max ( 1, j - ml - mu ) - 1, -1 ): x[i-1,k] = x[i-1,k] - a_lu[l+i-1,j-1] * x[j-1,k] else: # # Solve A' * x = b. # # Solve U' * x = b, overwriting b with x. # for i in range ( 1, nrhs + 1 ): for j in range ( 1, n + 1 ): temp = x[j-1,i-1] l = ml + mu + 1 - j for k in range ( max ( 1, j - ml - mu ), j ): temp = temp - a_lu[l+k-1,j-1] * x[k-1,i-1] x[j-1,i-1] = temp / a_lu[ml+mu,j-1] # # Solve L' * x = b, overwriting b with x. # if ( 0 < ml ): for j in range ( n - 1, 0, -1 ): lm = min ( ml, n - j ) for k in range ( 1, nrhs + 1 ): t = 0.0 for i in range ( 0, lm ): t = t + x[j+i,k-1] * a_lu[kd+i,j-1] x[j-1,k-1] = x[j-1,k-1] - t l = pivot[j-1] for i in range ( 1, nrhs + 1 ): t = x[l-1,i-1] x[l-1,i-1] = x[j-1,i-1] x[j-1,i-1] = t return x, info def r8gb_trs_test ( ): #*****************************************************************************80 # ## R8GB_TRS_test tests R8GB_TRS. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # n = 10 m = n ml = 1 mu = 2 nrhs = 1 print ( '' ) print ( 'R8GB_TRS_test' ) print ( ' R8GB_TRS solves a linear system factored by R8GB_TRS.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) # # Set the matrix. # a = r8gb_random ( m, n, ml, mu ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # b = r8gb_mv ( m, n, ml, mu, a, x ) # # Factor the matrix. # a_lu, pivot, info = r8gb_trf ( m, n, ml, mu, a ) if ( info != 0 ): print ( '' ) print ( 'R8GB_TRS_test - Fatal error!' ) print ( ' R8GB_TRF declares the matrix is singular!' ) print ( ' The value of INFO is ', info ) raise Exception ( 'r8gb_trs_test(): Fatal error!' ) # # Solve the linear system. # Note that, because of quirks in MATLAB, we need to copy our vector-based # data to and from 2D arrays. # b_mat = r8vec_to_r8ge ( n, nrhs, b ) x_mat, info = r8gb_trs ( n, ml, mu, nrhs, 'N', a_lu, pivot, b_mat ) x = r8ge_to_r8vec ( n, nrhs, x_mat ) r8vec_print ( n, x, ' Solution:' ) # # Set the desired solution. # x = r8vec_indicator1 ( n ) # # Compute the corresponding right hand side. # job = 1 b = r8gb_mu ( n, ml, mu, a_lu, pivot, x, job ) # # Solve the linear system. # b_mat = r8vec_to_r8ge ( n, nrhs, b ) x_mat, info = r8gb_trs ( n, ml, mu, nrhs, 'T', a_lu, pivot, b_mat ) x = r8ge_to_r8vec ( n, nrhs, x_mat ) r8vec_print ( n, x, ' Solution to transposed system:' ) return def r8gb_zeros ( m, n, ml, mu ): #*****************************************************************************80 # ## R8GB_ZEROS zeros an 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. # # LINPACK and LAPACK band storage requires that an extra ML # superdiagonals be supplied to allow for fillin during Gauss # elimination. Even though a band matrix is described as # having an upper bandwidth of MU, it effectively has an # upper bandwidth of MU+ML. This routine assumes it is setting # up an unfactored matrix, so it only uses the first MU upper bands, # and does not place nonzero values in the fillin bands. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows of the matrix. # M must be positive. # # Input, integer N, the number of columns 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 min(M,N)-1. # # Output, real A(2*ML+MU+1,N), the R8GB matrix. # import numpy as np a = np.zeros ( [ 2 * ml + mu + 1, n ] ) return a def r8gb_zeros_test ( ): #*****************************************************************************80 # ## R8GB_ZEROS_test tests R8GB_ZEROS. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2016 # # Author: # # John Burkardt # m = 5 n = 5 ml = 1 mu = 2 print ( '' ) print ( 'R8GB_ZEROS_test' ) print ( ' R8GB_ZEROS returns an R8GB zero matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_zeros ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB zero matrix:' ) return def r8ge_to_r8vec ( m, n, a ): #*****************************************************************************80 # ## r8ge_to_r8vec() copies an R8GE matrix to an R8VEC. # # Discussion: # # In C++ and FORTRAN, this routine is not really needed. In MATLAB, # a data item carries its dimensionality implicitly, and so cannot be # regarded sometimes as a vector and sometimes as an array. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 07 February 2016 # # Author: # # John Burkardt # # Input: # # integer M, N, the number of rows and columns in the array. # # real A(M,N), the array to be copied. # # Output: # # real X(M*N), the vector. # import numpy as np x = np.zeros ( m * n ) k = 0 for j in range ( 0, n ): for i in range ( 0, m ): x[k] = a[i,j] k = k + 1 return x def r8st_print ( m, n, nz_num, sym, row, col, a, title ): #*****************************************************************************80 # ## r8st_print() prints a r8st matrix. # # Discussion: # # The r8st storage format corresponds to the SLAP Triad format. # # The r8st storage format stores the row, column and value of each nonzero # entry of a sparse matrix. The entries may be given in any order. No # check is made for the erroneous case in which a given matrix entry is # specified more than once. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2015 # # Author: # # John Burkardt # # Input: # # integer M, N, the order of the matrix. # # integer NZ_NUM, the number of nonzero elements in the matrix. # # integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices # of the nonzero elements. # # integer SYM: the symmetry option. (ignored for now) # # real A(NZ_NUM), the nonzero elements of the matrix. # # string TITLE, a title. # r8st_print_some ( m, n, nz_num, row, col, a, 0, 0, m - 1, n - 1, title ) return def r8st_print_some ( m, n, nz_num, row, col, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## r8st_print_some() prints some of an r8st matrix. # # Discussion: # # The r8st storage format corresponds to the SLAP Triad format. # # The r8st storage format stores the row, column and value of each nonzero # entry of a sparse matrix. The entries may be given in any order. No # check is made for the erroneous case in which a given matrix entry is # specified more than once. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 September 2015 # # Author: # # John Burkardt # # Input: # # integer M, N, the order of the matrix. # # integer NZ_NUM, the number of nonzero elements in the matrix. # # integer ROW[NZ_NUM], COL[NZ_NUM], the row and column indices # of the nonzero elements. # # real A[NZ_NUM], the nonzero elements # of the matrix. # # integer ILO, JLO, IHI, JHI, the first row and # column, and the last row and column to be printed. # # string TITLE, a title. # import numpy as np incx = 5 index = np.zeros ( incx, dtype = np.int32 ) print ( '' ) print ( title ) # # Print the columns of the matrix, in strips of 5. # for j2lo in range ( jlo, jhi + 1, incx ): j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n - 1 ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo print ( '' ) print ( ' Col: ', end = '' ) for j in range ( j2lo, j2hi + 1 ): print ( '%7d ' % ( j ), end = '' ) print ( '' ) print ( ' Row' ) print ( ' ---' ) # # Determine the range of the rows in this strip. # i2lo = max ( ilo, 0 ) i2hi = min ( ihi, m - 1 ) for i in range ( i2lo, i2hi + 1 ): # # Print out (up to) 5 entries in row I, that lie in the current strip. # nonzero = False for j2 in range ( 0, inc ): index[j2] = -1 for k in range ( 0, nz_num ): if ( i == row[k] and j2lo <= col[k] and col[k] <= j2hi ): j2 = col[k] - j2lo + 1 if ( a[k] != 0.0 ): index[j2-1] = k nonzero = True if ( nonzero ): print ( '%5d ' % ( i ), end = '' ) for j2 in range ( 0, inc ): if ( 0 <= index[j2] ): aij = a[index[j2]] else: aij = 0.0 print ( '%14g' % ( aij ), end = '' ) print ( '' ) return def 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_r8gb_test ( ): #*****************************************************************************80 # ## r8vec_to_r8gb_test() tests r8vec_to_r8gb(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2016 # # Author: # # John Burkardt # m = 5 n = 8 ml = 2 mu = 1 print ( '' ) print ( 'r8vec_to_r8gb_test()' ) print ( ' r8vec_to_r8gb() converts an R8VEC to an R8GB matrix.' ) print ( '' ) print ( ' Matrix rows M = ', m ) print ( ' Matrix columns N = ', n ) print ( ' Lower bandwidth ML = ', ml ) print ( ' Upper bandwidth MU = ', mu ) a = r8gb_indicator ( m, n, ml, mu ) r8gb_print ( m, n, ml, mu, a, ' The R8GB indicator matrix:' ) x = r8gb_to_r8vec ( m, n, ml, mu, a ) k = 0 for j in range ( 0, n ): for i in range ( 0, 2 * ml + mu + 1 ): print ( '%4d %4d %4d %14g' % ( i, j, k, x[k] ) ) k = k + 1 a = r8vec_to_r8gb ( m, n, ml, mu, x ) r8gb_print ( m, n, ml, mu, a, ' The recovered R8GB indicator matrix:' ) return def r8vec_to_r8ge ( m, n, x ): #*****************************************************************************80 # ## r8vec_to_r8ge() copies an R8VEC into a R8GE matrix. # # Discussion: # # In C++ and FORTRAN, this routine is not really needed. In MATLAB, # a data item carries its dimensionality implicitly, and so cannot be # regarded sometimes as a vector and sometimes as an array. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 08 February 2016 # # Author: # # John Burkardt # # Input: # # integer M, N, the number of rows and columns in the array. # # real X(M*N), the vector to be copied into the array. # # Output: # # real A(M,N), the array. # import numpy as np a = np.zeros ( [ m, n ] ) k = 0 for j in range ( 0, n ): for i in range ( 0, m ): a[i,j] = x[k] k = k + 1 return a def r8vec2_print ( a1, a2, title ): #*****************************************************************************80 # ## r8vec2_print() prints an R8VEC2. # # Discussion: # # An R8VEC2 is a dataset consisting of N pairs of real values, stored # as two separate vectors A1 and A2. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 20 June 2020 # # Author: # # John Burkardt # # Input: # # integer N, the number of components of the vector. # # real A1(N), A2(N), the vectors to be printed. # # string TITLE, a title. # n = len ( a1 ) print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( ' %6d: %12g %12g' % ( i, a1[i], a2[i] ) ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 August 2019 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) r8gb_test ( ) timestamp ( )