# This is a blind linear equation solver.
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 :) #
(INT n = 3;
MODE R = LONG LONG REAL;
MODE M = [n, n+1]R;
MODE A = [n]R;
PROC pr = (M m)VOID: (FOR i TO n DO FOR j TO n+1 DO
print((SHORTEN SHORTEN m[i, j], " ")) OD; print(newline) OD; print(newline));
PROC slv = (REF M m)A: (
FOR j TO n DO # select diagonal element #
R p = 1.0/m[j, j];
FOR k FROM j+1 TO n+1 DO m[j, k] *:= p OD;
FOR k FROM j+1 TO n DO
FOR h FROM j+1 TO n+1 DO
m[k, h] -:= m[k, j]*m[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 #
m[k, n+1] -:= m[k, j]*m[j, n+1] OD OD;
m[,n+1]);
# posit secret answer: (3, -3, 7) #
M x := (
(2, 4, 6, 36),
(4, 7, 6, 33),
(7, 4, 3, 30));
pr(x);
print(slv(x));
SKIP)