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