Solve a system of linear equations (Intermediate to advanced level)
Contributed by: Doug Rickman, Global Hydrology and Climate Center, MSFC, NASA
This is aimed at intermediate to advanced REXX users.
The method of solution is based on the Gauss_Jordan algorithm as outlined in "Numerical Recipes." It is implemented with a suite of linear algebra routines. Rexx is very nice for this type of problem as precision becomes a very big concern with larger systems of equations. With most languages special methods must be used to maintain the precision. With Rexx one can simply increase the NUMERIC DIGITS value. This code also makes use of the "expose list" passed as a variable, which in turn containing a list of variables to be exposed.
Two test matices are provided at the end of the code.
/* Solve a square array using the Gauss-Jordan algorithm as outlined in */
/* numerical recipes. Also a collection of linear algebra routines. */
/* Implemented by Doug Rickman March 13, 1998 */
/* Input 2 matrices, a square array and an array of solution vectors. */
signal on Halt
signal on NotReady
if rxfuncquery('rexxlibregister') then do /* this will start rexxlib */
call rxfuncadd 'rexxlibregister', 'rexxlib', 'rexxlibregister'
call rexxlibregister
end
if rxfuncquery('sysloadfuncs') then do /* this will start rexxutil */
CALL RxFuncAdd 'SysLoadFuncs', 'RexxUtil', 'SysLoadFuncs'
CALL SysLoadFuncs
end
arg in1 in2
/* Read in the A. data */
i=0
do while lines(in1)>0
i=i+1
data=linein(in1)
j=0
do while data<>''
j=j+1
parse var data A.i.j data
end /* do while data ... */
end /* do while lines(in1) ... */
rc=lineout(in1)
A.0=i*j
NRowsA=i
NColsA=j
/* Read in the B. (solution) vectors */
i=0
do while lines(in2)>0
i=i+1
data=linein(in2)
j=0
do while data<>''
j=j+1
parse var data B.i.j data
end /* do while data ... */
end /* do while lines(inData) ... */
rc=lineout(in2)
B.0=i*j
NRowsB=i
NColsB=j
/* Echo the input arrays. */
elist= 'A. NRowsA NColsA'
call ShowMatrix elist
elist= 'B. NRowsB NColsB'
call ShowMatrix elist
/* Backup the original A array. */
do i=1 to NRowsA
do j=1 to NColsA
AA.i.j=A.i.j
end j
end i
/* call GaussJordan A. B. */ /* Call to un-interpreted subroutine */
/* now call the Gauss-Jordan routine */
elist= 'A. NRowsA NColsA B. NRowsB NColsB'
call GaussJordan elist
/* Critical quality check!!!!! */
/* Multiply original A. by A inverse to check for identity. */
elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
say 'Checking A x A(inverse):'
call MatrixMultiply elist
elist= 'C. NRowsC NColsC'
rc=IndentityCheck()
if rc=1 then do
say 'The result is close to the identity matrix. All is well.'
end
else do
say 'The inverse of A. may not be numerically precise enough. ',
||'You should examine the values in A. x A.inverse, in theory ',
||'this should be the identity matrix.'
elist= 'C. NRowsC NColsC'
call ShowMatrix elist
end
/* Optional quality checks. Set "QualityChecks=YES" to use. */
QualityChecks='NO'
if QualityChecks=YES then call QualityChecks
/* Print the solution matrix */
say 'Solution vectors: '
elist= 'B. NRowsB NColsB'
call ShowMatrix elist
return
/* This version uses the interpret statement, thus allowing arrays other than*/
/* A. and B. to be used. */
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - GaussJordan: -------------*/
/* Solve a square array using the Gauss-Jordan algorithm as outlined in */
/* numerical recipes. */
/* Implemented by Doug Rickman March 13, 1998 */
/* elist holds: */
/* name of the first array, */
/* the variable holding the number of rows in the first array, */
/* the variable holding the number columns in the first array. */
/* name of the second array, */
/* the variable holding the number of rows in the second array, */
/* the variable holding the number columns in the second array. */
/* */
/* A. NRowsA NColsA B. NRowsB NColsB */
GaussJordan:
procedure expose (elist)
parse var elist VArrayName1 VNRows1 VNCols1 VArrayName2 VNRows2 VNCols2
NRows1=value(VNRows1)
NCols1=value(VNCols1)
ArrayName1=strip(VArrayName1,'t','.')
NRows2=value(VNRows2)
NCols2=value(VNCols2)
ArrayName2=strip(VArrayName2,'t','.')
N=NRows1 /* number of elements and rows */
M=NCols2 /* right hand vactors is an array N by M */
NMax=50
do j=1 to N
IPIV.j=0
end j
do i=1 to N
BIG=0
do j=1 to N
if IPIV.j <> 1 then do k=1 to N
interpret 'Temp='Arrayname1'.j.k'
if IPIV.k=0 & abs(Temp)>=BIG then do
interpret 'BIG=abs('Arrayname1'.j.k)'
irow=j
icol=k
end
else if IPIV.k > 1 then do
say 'Singular matrix! Stop 1'
exit
end
end k
end j
IPIV.icol=IPIV.icol+1
if irow<>icol then do
do /*14*/ L=1 to N
interpret 'DUM='Arrayname1'.irow.L'
interpret Arrayname1'.irow.L='Arrayname1'.icol.L'
interpret Arrayname1'.icol.L=DUM'
end L
do L=1 to M
interpret 'DUM='Arrayname2'.irow.L'
interpret Arrayname2'.irow.L='Arrayname2'.icol.L'
interpret Arrayname2'.icol.L=DUM'
end L
end /* if irow<>icol then do ... */
INDXR.i=irow
INDXC.i=icol
interpret 'Temp='Arrayname1'.icol.icol'
if Temp=0 then do
say 'Singular matrix! Stop 2.'
exit
end
interpret 'PIVINV=1/'Arrayname1'.icol.icol'
interpret Arrayname1'.icol.icol=1'
do L=1 to N
interpret Arrayname1'.icol.L='Arrayname1'.icol.L*PIVINV'
end L
do L=1 to M
interpret Arrayname2'.icol.L='Arrayname2'.icol.L*PIVINV'
end L
do LL=1 to N
if LL \= icol then do
interpret 'DUM='Arrayname1'.LL.icol'
interpret Arrayname1'.LL.icol=0'
do L=1 to N
interpret Arrayname1'.LL.L='Arrayname1'.LL.L-'Arrayname1'.icol.L*DUM'
end L
do L=1 to M
interpret Arrayname2'.LL.L='Arrayname2'.LL.L-'Arrayname2'.icol.L*DUM'
end L
end /* if LL \= icol then do ... */
end LL
end i
/* Unscramble */
do L=N to 1 by -1
if INDXR.L \= INDXC.L then do K=1 to N
INDXRL=INDXR.L
INDXCL=INDXC.L
interpret 'DUM='Arrayname1'.K.INDXRL'
interpret Arrayname1'.K.INDXRL='Arrayname1'.K.INDXCL'
interpret Arrayname1'.K.INDXCL=DUM'
end K
end L
return
/* --- end subroutine - GaussJordan: -------------*/
/* --------------------------------------------------------------------------*/
/* un-interpreted subroutine. retained for clarity and comparison. */
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - GaussJordan: -------------*/
GaussJordan:
procedure expose A. B. NRowsA NColsB
N=NRowsA /* number of elements and rows */
M=NColsB /* right hand vactors is an array N by M */
NMax=50
do j=1 to N
IPIV.j=0
end j
do i=1 to N
BIG=0
do j=1 to N
if IPIV.j <> 1 then do k=1 to N
if IPIV.k=0 & abs(A.j.k)>=BIG then do
BIG=abs(A.j.k)
irow=j
icol=k
end
else if IPIV.k > 1 then do
say 'Singular matrix! Stop 1'
exit
end
end k
end j
IPIV.icol=IPIV.icol+1
if irow<>icol then do
do /*14*/ L=1 to N
DUM=A.irow.L
A.irow.L=A.icol.L
A.icol.L=DUM
end L
do L=1 to M
DUM=B.irow.L
B.irow.L=B.icol.L
B.icol.L=DUM
end L
end /* if irow<>icol then do ... */
INDXR.i=irow
INDXC.i=icol
if A.icol.icol=0 then do
say 'Singular matrix! Stop 2.'
exit
end
PIVINV=1/A.icol.icol
A.icol.icol=1
do L=1 to N
A.icol.L=A.icol.L*PIVINV
end L
do L=1 to M
B.icol.L=B.icol.L*PIVINV
end L
do LL=1 to N
if LL \= icol then do
DUM=A.LL.icol
A.LL.icol=0
do L=1 to N
A.LL.L=A.LL.L-A.icol.L*DUM
end L
do L=1 to M
B.LL.L=B.LL.L-B.icol.L*DUM
end L
end /* if LL \= icol then do ... */
end LL
end i
/* Unscramble */
do L=N to 1 by -1
if INDXR.L \= INDXC.L then do K=1 to N
INDXRL=INDXR.L
INDXCL=INDXC.L
DUM=A.K.INDXRL
A.K.INDXRL=A.K.INDXCL
A.K.INDXCL=DUM
end K
end L
return
/* --- end subroutine - GaussJordan: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - MatrixAdd: -------------*/
/* Returns a 1 if successful, result is in C. */
/* Returns a 2 if the matrices are not the same size. */
MatrixAdd:
procedure expose A. B. C. NRowsA NColsA NRowsB NColsB NColsC NRowsC
if NRowsA=NRowsB & NColsA=NColsB then do
NColsC=NColsA
NRowsA=NRowsB
do i=1 to NRows
do j=1 to NCols
C.i.j=A.i.j+B.i.j
end j
end i
return 1
end /* if ... */
else return 2
/* --- end subroutine - MatrixAdd: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - MatrixMultiply: -------------*/
/* Returns a 1 if successful, result is in C. */
/* Returns a 2 if the matrices are not the correct size. */
/* elist holds: */
/* name of the first array, */
/* the variable holding the number of rows in the first array, */
/* the variable holding the number columns in the first array. */
/* name of the second array, */
/* the variable holding the number of rows in the second array, */
/* the variable holding the number columns in the second array. */
/* name of the result array, */
/* the variable holding the number of rows in the result array, */
/* the variable holding the number columns in the result array. */
/* */
/* A. B. C. NRowsA NColsA NRowsB NColsB NColsC NRowsC */
MatrixMultiply:
procedure expose (elist)
parse var elist VArrayName1 VNRows1 VNCols1 ,
VArrayName2 VNRows2 VNCols2 ,
VArrayNameR VNRowsR VNColsR
NRows1=value(VNRows1)
NCols1=value(VNCols1)
ArrayName1=strip(VArrayName1,'t','.')
NRows2=value(VNRows2)
NCols2=value(VNCols2)
ArrayName2=strip(VArrayName2,'t','.')
ArrayNameR=strip(VArrayNameR,'t','.')
if NCols1=NRows2 then do
NRowsR=Nrows1
NColsR=NCols2
do i=1 to NRowsR
do j=1 to NColsR
interpret ArrayNameR'.i.j=0'
do k=1 to NCols1
interpret ArrayNameR'.i.j='ArraynameR'.i.j+('ArrayName1'.i.k*'ArrayName2'.k.j)'
end k
end j
end i
end /* if ... */
interpret VNRowsR'=NRowsR'
interpret VNColsR'=NColsR'
return 1
end
else return 2
/* --- end subroutine - MatrixMultiply: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - MatrixMultiplyByScalar: -------------*/
MatrixMultiplyByScalar:
procedure expose A. NRowsA NColsA Scalar
do i=1 to 3
row=''
do j=1 to 3
row=row A.i.j
end j
say row
end i
return
/* --- end subroutine - MatrixMultiplyByScalar: -------------*/
/* --------------------------------------------------------------------------*/
/* This is a start on a callable version of above. It needs to be told where*/
/* to write the result array. March 31, 1998 */
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - MatrixMultiplyByScalar: -------------*/
/* elist holds: */
/* name of the array, */
/* the variable holding the number of rows, */
/* the variable holding the number columns. */
/* the variable holding the scalar. */
MatrixMultiplyByScalar:
procedure expose procedure expose (elist)
parse var elist ArrayName1 VNRows VNCols Scalar
NRows=value(VNRows1)
NCols=value(VNCols1)
Arrayname=strip(Arrayname1,'t','.')
do i=1 to NRows
row=''
do j=1 to NCols
interpret 'row=row 'Arrayname'i.j
end j
say row
end i
return
/* --- end subroutine - MatrixMultiplyByScalar: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - IndentityCheck: -------------*/
/* Check to see if matrix is an identity matrix (diagonal=1, other=0). */
/* 1 is returned for and identity matrix, a 0 otherwise. */
/* Precision can be increased by increasing NUMERIC DIGITS for the Gauss- */
/* Jordan subroutine. Also the tolerance for error can be increased by */
/* changing the NUMERIC FUZZ value in this routine. */
/* elist holds: */
/* name of the array, */
/* the variable holding the number of rows, */
/* the variable holding the number columns. */
IndentityCheck: procedure expose (elist)
parse var elist ArrayName1 VNRows1 VNCols1
NRows=value(VNRows1)
NCols=value(VNCols1)
Arrayname=strip(Arrayname1,'t','.')
n=digits()
numeric fuzz n-3
do i=1 to NRows
do j=1 to NCols
interpret 'test=1+'Arrayname'.i.j'
if i=j & test=2 then iterate
if test=1 then iterate
else do
numeric fuzz 0
return 0
end
end j
end i
numeric fuzz 0
return 1
/* --- end subroutine - IndentityCheck: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - RowsColumns -------------*/
RowsColumns: /* Determine the number of rows and columns in a 2D matrix */
procedure expose A.
rc=cvtails(A.,tails.)
irowMax=0
jcolMax=0
do i=1 to tails.0-1
parse var tails.i irow '.'jcol
irowMax=max(irow,irowMax)
jcolMax=max(jcol,jcolMax)
end i
return irowMax jcolMax
/* --- end subroutine - RowsColumns -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - Halt: -------------*/
Halt:
say 'This is a graceful exit from a Cntl-C'
exit
/* --- end subroutine - Halt: -------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - NotReady: -------------*/
NotReady:
say 'It would seem that you are pointing at non-existant data. Oops. Bye!'
exit
/* --- end subroutine - NotReady: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - ShowMatrix: -------------*/
/* elist holds: */
/* name of the array, */
/* the variable holding the number of rows, */
/* the variable holding the number columns. */
ShowMatrix: procedure expose (elist)
parse var elist ArrayName1 VNRows1 VNCols1
NRows=value(VNRows1)
NCols=value(VNCols1)
Arrayname=strip(Arrayname1,'t','.')
say 'Array: 'ArrayName
do i=1 to NRows
row=''
do j=1 to NCols
interpret 'row=row 'ArrayName'.i.j'
end j
say row
end i
return
/* --- end subroutine - ShowMatrix: -------------*/
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
/* --- begin subroutine - -------------*/
QualityChecks:
say
/* Print the inverse of A. */
say 'Inverse of A is '
elist= 'A. NRowsA NColsA'
call ShowMatrix elist
say
/* Multiply A inverse times the solution vector. The result should equal the*/
/* original B. matrix. */
elist= 'AA. NRowsA NRowsB B. NRowsB NColsB C. NRowsC NColsC'
say 'Inverse of A times the solution B: (should equal the original matrix B.)'
call MatrixMultiply elist
elist= 'C. NRowsC NColsC'
call ShowMatrix elist
say
/* Multiply original by inverse to check for identity. */
elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
say 'Original A times the inverse of A: (should equal identity matrix)'
call MatrixMultiply elist
elist= 'C. NRowsC NColsC'
call ShowMatrix elist
return
/* --- end subroutine - -------------*/
/* --------------------------------------------------------------------------*/
/* Test values */
matrix A 9 36 204
36 204 1296
204 1296 8772
matrix B 80.5
299
1697
solution matrix 12.1848484
-1.84653670
0.182900419
/* from Numerical Recipes Examples */
Array: A1
1 2 3 4 5
2 3 4 5 1
3 4 5 1 2
4 5 1 2 3
5 1 2 3 4
Array: B1
1 1
1 2
1 3
1 4
1 5
Solution vectors:
0.0666666670 1.0
0.0666666666 0
0.0666666666 0
0.066666667 0
0.0666666656 0
*/
