# Algol 68 Matrix Inversion #
MODE R = REAL;
PROC inv = ([,]R m) [,]R: (INT n = 1 UPB m;
   IF n /= 2 UPB m THEN print("not square") FI;
   [n, 2*n]R am; FOR i TO n DO FOR j TO n DO
      am[i, j] := m[i, j]; am[i, j+n] := (i=j|1|0) OD OD;
   FOR i TO n DO INT im := 0; R mv := ABS am[i, i];
     FOR j FROM i+1 TO n DO R ae = ABS am[j, i];
       IF mv < ae THEN mv := ae; im := j FI OD;
     IF mv = 0.0 THEN print("Singular") FI;
     IF im > i THEN FOR n FROM i TO 2*n DO
        R s = am[i, n]; am[i, n] := am[im, n]; am[im, n] := s OD FI;
   (R r = 1.0 / am[i, i]; FOR j FROM i TO 2*n DO am[i, j] *:= r OD);
   FOR k FROM i+1 TO n DO R f = am[k, i]; 
     FOR j FROM i+1 TO 2*n DO am[k, j] -:= f * am[i, j] OD OD OD;
   FOR i TO n DO FOR j FROM i+1 TO n DO R p = am[i, j];
     FOR k FROM i+1 TO 2*n DO am[i, k] -:= am[j, k] * p OD OD OD;
   am[,n+1:2*n]);

(1|(INT n = 10; [n, n]R hil; FOR i TO n DO FOR j TO n DO hil[i, j] := R(1.0)/(i+j-1) OD OD;
  ([,]R i = inv(hil); print((i, newline, inv(i))))),
print(inv(inv([,]R((3, 6, 5, 2), (3, -4, 6, 2), (7, 5, 4, 2), (3, 6, 7, 1))))),
print(inv([,]R((0,0,0,1), (0,0,1,0), (1,0,0,0), (0,1,0,0))))
)