# This is a blind matrix inverter.
 It is oblivious to ill conditioning.
 Note faith in calculation of p. 
 It does not even avoid introducing ill conditioning! 
 It almost always muddles thru :)
 It is simpler than it is fast. #
(INT n = 5;
MODE R = LONG LONG REAL;
MODE RF = PROC(R)R;
MODE M = [n, n]R;
MODE A = [n]R;
PROC pr = (M m)VOID: (FOR i TO n DO FOR j TO n DO
   print((SHORTEN SHORTEN m[i, j], " ")) OD; print(newline) OD; print(newline));
PROC inv = (M m)M: ([n, 2*n]R x; x[ , 1:n] := m;
  FOR i TO n DO FOR j TO n DO x[i, j+n] := (i=j|1|0) OD OD;
  FOR j TO n DO # select diagonal element #
    (R p = 1.0/x[j, j];
     FOR k FROM j+1 TO 2*n DO x[j, k] *:= p OD);
     FOR k FROM j+1 TO n DO
       FOR h FROM j+1 TO 2*n DO
         x[k, h] -:= x[k, j]*x[j, h] OD OD OD;
# At this point the entries on the diagonal are noise and logically 1.
  The entries below the diagonal are noise and logically 0. #
  FOR j FROM n DOWNTO 1 DO # selects diagonal element #
    FOR k TO j-1 DO # selects row #
      FOR h FROM n+1 TO 2*n DO
        x[k, h] -:= x[k, j]*x[j, h] OD OD OD;
  x[,n+1:2*n]);

# Matrix multiply #
OP * = (M x, M y)M: ([n, n]R a;
  FOR i TO n DO
    FOR j TO n DO R t := 0;
      FOR h TO n DO t +:= x[i, h]*y[h, j] OD;
      a[i, j] := t OD OD;
  a);

M x := (
 (2, 4, 6, 43.31, 8), 
  (4, 7, 6, -2, 1.1),
  (7, 4, 3, 4, 2),
  (-2, 4.4, 2, -9.4, 0),
  (3, 5, 3.2, 6.1, 4));

pr(x);
(M y = inv(x); pr(y); pr(x*y)))