TAOCP Vol 1, Section 2.2.6, Exercise 15. MMIX version. 
(Pivoting of a matrix implemented  as a two dimensional linked list.)

Conversion from MIX by Ladislav Sladecek sladecek@users.sourceforge.net


** FILES

README			this file

pivot_book.mms		The MMIX code (just the procedure PivotStep 
			which is printed in the book)

pivot_run.mms 		Complete runnable MMIX code with main() and 
			routines for creating, loading and dumping
			matrices and node allocation. PivotStep is tested
			indirectly by running MatrixInversion (a quick
			and dirty procedure for matrix inversion)	
			which calls the PivotStep.

pivot_test.pl		A perl wrapper which generates a sequence of random
			matrices, compiles pivot_run.mms, and runs it. The
			resulting matrix is multiplied with the original matrix
			and the difference is evaluated.

mmix-sim.ch		The testing routines depend on a kludge for floating 
			point I/O which is in the file mmix-sim.ch.  
			The pivot_book code does not depend on this kludge in 
			any way. 


** HOWTO run this


1) Tested on RedHat Linux 7.2; should work on most other systems.

2) Prerequisites: 

a) mmix source code (I used version of October 14, 2001).

b) perl module Math::MatrixReal, can be downloaded from www.cpan.org
or better can be installed directly by running:

perl -MCPAN -e shell;	
install Math::MatrixReal
	
3) Unpack the mmix source in
~/mmixmasters/support/mmix_with_printfloat. Unpack this archive in
~/mmixmasters/chapter02/pivot

4) Copy mmix-sim.ch to ~/mmixmasters/support/mmix_with_printfloat

5) cd  ~/mmixmasters/support/mmix_with_printfloat

6) make 

7) mv mmix mmix_with_printfloat

7) cd ~/mmixmasters/chapter02/pivot 

8) perl pivot_test.pl 
(Read the source of pivot_test.pl for explanation of arguments.)


** DISCUSSION

1) The original MIX code used ROW(basecol[]) for PTR[]. My MMIX code
uses a separate array ptr[]  because the values of ROW(basecol[])
should be negative to function as a sentinel. The algorithm depends on
the sentinel. I suspect there are few bugs in the original MIX code.

2) I used a slightly more complicated allocation scheme to avoid necessity
to create a linked list of free nodes in advance.

3) The most complicated was the choice of pointer size (UP and LEFT pointers).
Assuming that the values are 64 bit floats, all  nodes must be 
aligned  toward octabyte. Thus three options make sense:

a) both pointers and indexes are octabytes. This is the simplest 
option allowing for almost unlimited matrix size, but is rather 
ineffective for small size matrices. We will waste 400% of memory 
for pointers. This option was chosen here for its simplicity.

b) pointers and indexes are tetrabytes which count nodes starting 
from the  base address :nodebase. The overhead is still 200% and 
the matrices  are limited to about 4 billions  nodes. 

c) pointers and indexes are wydes; the overhead is just 100% but the
 number of nodes is limited by 65536 - 2n. The effective limit for n
is just few hundreds.


** RANDOM IDEAS

1) Can someone explain to me why the pointers are oriented in the
UP-LEFT direction instead of the more natural DOWN-RIGHT direction?
There must be a reason for this, I am just to stupid to find it
myself. The UP-RIGHT and DOWN-LEFT combinations do not make sense
because they will make matrix multiplication a nightmare, but the
DOWN-RIGHT direction would allow to print the matrix in the natural order. 

2) It would be nice to have few extensions to mmix and mmixal:

a) floating point constants in mmixal.

b) a library of the most useful functions like integer and float IO.

c) a mmixal lint which will print warnings (see MY MOST STUPID ERRORS).


3) Why the matrix inversion works: 

(Mathematicians are kindly asked to SKIP THIS!)

The simplest method to invert a matrix $A$ by Gauss elimination is to
write an unary matrix $B$ of the same dimension and perform the same
operations in parallel on both the matrices. While we transform $A$ to unary
form by Gauss elimination; we get the inverse of $A$ in $B$.

We can do the algorithm in-place storing only one matrix because we know
that half of the elements are always zeroes and ones.


** MY MOST STUPID ERRORS

1) UP IS VAL + 1 

(Any non-alphanumeric character starts a comment. Hence any
additional space (a typo) makes an error which is not detected by the 
compiler.)

2) POP 1

POP 1 is interpreted as POP 0,1 instead of POP 1,0.


Ladislav Sladecek