# 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)))