#! /usr/bin/env python
#
def r8mat_solve ( n, nrhs, a ):
#*****************************************************************************80
#
## R8MAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 17 October 2015
#
# Author:
#
# John Burkardt
#
# Parameters:
#
# Input, integer N, the order of the matrix.
#
# Input, integer NRHS, the number of right hand sides. NRHS
# must be at least 0.
#
# Input, real A(N,N+NRHS), contains in rows and
# columns 1 to N the coefficient matrix, and in columns N+1 through
# N+NRHS, the right hand sides.
#
# Output, real A(N,N+NRHS), the coefficient matrix
# area has been destroyed, while the right hand sides have
# been overwritten with the corresponding solutions.
#
# Output, integer INFO, singularity flag.
# 0, the matrix was not singular, the solutions were computed;
# J, factorization failed on step J, and the solutions could not
# be computed.
#
info = 0
for j in range ( 0, n ):
#
# Choose a pivot row IPIVOT.
#
ipivot = j
apivot = a[j,j]
for i in range ( j + 1, n ):
if ( abs ( apivot ) < abs ( a[i,j] ) ):
apivot = a[i,j]
ipivot = i
if ( apivot == 0.0 ):
info = j
return a, info
#
# Interchange.
#
for k in range ( 0, n + nrhs ):
temp = a[ipivot,k]
a[ipivot,k] = a[j,k]
a[j,k] = temp
#
# A(J,J) becomes 1.
#
a[j,j] = 1.0
for k in range ( j + 1, n + nrhs ):
a[j,k] = a[j,k] / apivot;
#
# A(I,J) becomes 0.
#
for i in range ( 0, n ):
if ( i != j ):
factor = a[i,j]
a[i,j] = 0.0
for k in range ( j + 1, n + nrhs ):
a[i,k] = a[i,k] - factor * a[j,k]
return a, info
def r8mat_solve_test ( ):
#*****************************************************************************80
#
## R8MAT_SOLVE_TEST tests R8MAT_SOLVE.
#
# Licensing:
#
# This code is distributed under the GNU LGPL license.
#
# Modified:
#
# 29 February 2016
#
# Author:
#
# John Burkardt
#
import numpy as np
import platform
from r8mat_print import r8mat_print
n = 3
rhs_num = 2
#
# Each row of this definition is a COLUMN of the matrix.
#
a = np.array ( [ \
[ 1.0, 2.0, 3.0, 14.0, 7.0 ], \
[ 4.0, 5.0, 6.0, 32.0, 16.0 ], \
[ 7.0, 8.0, 0.0, 23.0, 7.0 ] ] )
print ( '' )
print ( 'R8MAT_SOLVE_TEST' )
print ( ' Python version: %s' % ( platform.python_version ( ) ) )
print ( ' R8MAT_SOLVE solves linear systems.' )
#
# Print out the matrix to be inverted.
#
r8mat_print ( n, n + rhs_num, a, ' The linear system:' )
#
# Solve the systems.
#
a, info = r8mat_solve ( n, rhs_num, a )
if ( info != 0 ):
print ( '' )
print ( ' The input matrix was singular.' )
print ( ' The solutions could not be computed.' )
return
r8mat_print ( n, n + rhs_num, a, ' Factored matrix and solutions:' )
#
# Terminate.
#
print ( '' )
print ( 'R8MAT_SOLVE_TEST:' )
print ( ' Normal end of execution.' )
return
if ( __name__ == '__main__' ):
from timestamp import timestamp
timestamp ( )
r8mat_solve_test ( )
timestamp ( )