march@m.cs.uiuc.edu (11/16/89)
This is the HP28 program I wrote to emulate the RowReduce[] function of
Mathematica. The only argument is a M X N matrix A. A will then be
reduced to echelon form using Gaussian-Jordan Elimination. I have tested
this routine on everything from 2 X 2 matrices to 7 X 7 matrices with
correct results. This has been optimized (relative to the original) somewhat.
The first trial runs using a 7 X 7 matrix took close to 10 mins. Now, it
takes closer to 2 mins. Using 2 X 2 matrices usually takes about 10 seconds.
Of course, since this is one HUGE loop, optimizations tend to be multiplied
quite a bit. Please feel free to send me any and all optimizations. I hate
using STO to store variables but I am somewhat new to local variables and
don't quite know all the ways in which you can use them yet.
In addition, this algorithm is not smart in that it tries to reduce each
and every column even if it produces fractions. Hence, it is not a perfect
duplication of Mathematica's RowReduce[] function. This might be changed in
future versions after it is more optmized. I might also be tempted to
integrate this into a linear equation solver and use ISOL and others do
back substitution when full reduction has not taken place. Any thoughts?
I might also mention that this works for complex matrices as well as
real matrices.
Pascal code is also available if anyone wants it. Just email me at the address
in my .sig.
Example input: output:
1: [[ 3 4 ] 1: [[ 1 0 ]
[ 5 7 ]] [ 0 1 ]]
1: [[ 237 38 116 42 104 29 51 ] 1: [[ 1 0 0 0 0 0 0 ]
[ 290 45 142 51 132 32 65 ] [ 0 1 0 0 0 0 0 ]
[ 112 18 55 20 49 14 24 ] [ 0 0 1 0 0 0 0 ]
[ 253 40 124 45 112 30 55 ] [ 0 0 0 1 0 0 0 ]
[ 106 17 52 19 47 13 23 ] [ 0 0 0 0 1 0 0 ]
[ 257 41 126 46 114 31 56 ] [ 0 0 0 0 0 1 0 ]
[ 224 36 110 40 100 28 49 ]] [ 0 0 0 0 0 0 1 ]]
Another example:
You want to solve this linear system non-trivially ...
47x1 - 73x2 + 56x3 + 21x4 = 0
19x1 + 81x2 - 17x3 - 99x4 = 0
53x1 + 62x2 + 39x3 + 25x4 = 0
You enter: The output:
1: [[ 47 -73 56 21 ] 1: [[ 1 0 0 -4.31885078467 ]
[ 19 81 -17 -99 ] [ 0 1 0 .867684070361 ]
[ 53 62 39 25 ]] [ 0 0 1 5.13083792885 ]]
This indicates that x1 = 4.3..., x2 = -0.867... and x3 = -5.13...
From this x4 can then be determined. You will note that this system
couldn't be done using matrix division.
Well, after all that jibber-jabber ... here it is ...
ROWR [F79D]
<< DUP A STO SIZE LIST-> DROP -> m n
<< 1 1 P STO Q STO WHILE P m <= Q n <= AND REPEAT
0 'F' STO 0 'C' STO P 'K' STO
P m FOR i A i Q 2 ->LIST GET ABS -> x
<< IF x C > THEN x 'C' STO i 'K' STO END >>
NEXT
IF C .00001 < THEN 1 'F' STO Q 1 + 'Q' STO END
IF F 0 == THEN
1 n FOR j A P j -> a sr sc dr dc
<< a a sr sc 2 ->LIST GET dr dc 2 ->LIST SWAP PUT
a dr dc 2 ->LIST GET sr sc 2 ->LIST SWAP PUT >> 'A' STO
NEXT
A P Q 2 ->LIST GET 'L' STO
1 n FOR j A P j 2 ->LIST A P j 2 ->LIST GET L / PUT 'A' STO NEXT
1 m FOR i A i Q 2 ->LIST GET 'L' STO
1 n FOR j IF i P != THEN A i j 2 ->LIST A i j 2 ->LIST GET
A P j 2 ->LIST GET L * - PUT 'A' STO END
NEXT
NEXT
P 1 + 'P' STO Q 1 + 'Q' STO
END
END A
{ A C P Q L F X K } PURGE >>
>>
[I hate STO(!) ... but it works so .... I won't complain that much. :-) ]
-Steve
===============================================================================
Steve March (H) (217)328-5176/328-5230 (W) 333-7408
Department of Computer Science, University of Illinois
march@cs.uiuc.edu {uunet|convex|pur-ee}!uiucdcs!march
/* You are not expected to understand this. */ - UNIX V6 kernel source
"Time and space are modes by which we think and not conditions in which
we live." - Albert Einstein
===============================================================================