Flexowriter SFD Algol-60 example

This is an example of an Algol-60 program for the Electrologica X8, as it must have appeared in 1968 on the Flexowriter SFD, when it was typed (and simulaneously punched in papertape) by our late colleague Walter Hoffmann (1944-2013).
The various non-ASCII symbols used by the Flexowriter are available in Unicode. The easiest way to render Unicode is to embed it in an HTML document, as we have done here.

begin comment testprogramma voor qrivalsym2, 080268, R1089, W. Hoffmann;
integer i, n;
real eps;
array em[0:5];
comment mca 2011;

procedure elmveccol(l, u, i, a, b, x); value l, u, i, x;
integer l, u, i; real x; array a, b;
for l:= l step 1 until u do a[l]:= a[l] + b[l,i] × x;

comment mca 2012;

procedure elmcolvec(l, u, i, a, b, x); value l, u, i, x;
integer l, u, i; real x; array a, b;
for l:= l step 1 until u do a[l,i]:= a[l,i] + b[l] × x;

comment mca 2013;

procedure elmcol(l, u, i, j, a, b, x); value l, u, i, j, x;
integer l, u, i, j; real x; array a, b;
for l:= l step 1 until u do a[l,i]:= a[l,i] + b[l,j] × x;

comment mca 2300;

procedure tfmsymtri2(a, n, d, b, bb, em); value n; integer n;
array a, b, bb, d, em;
begin integer i, j, r, r1;
real w, x, a1, b0, bb0, d0, machtol, norm;
norm:= 0;
for j:= 1 step 1 until n do
begin w:= 0;
for i:= 1 step 1 until j do w:= abs(a[i,j]) + w;
for i:= j + 1 step 1 until n do w:= abs(a[j,i]) + w;
if w > norm then norm:= w
end;
machtol:= em[0] × norm; em[1]:= norm; r:= n;
for r1:= n - 1 step - 1 until 1 do
begin d[r]:= a[r,r]; x:= tammat(1, r - 2, r, r, a, a);
a1:= a[r1,r]; if sqrt(x) ≤ machtol then
begin b0:= b[r1]:= a1; bb[r1]:= b0 × b0; a[r,r]:= 1 end
else
begin bb0:= bb[r1]:= a1 × a1 + x;
b0:= if a1 > 0 then - sqrt(bb0) else sqrt(bb0);
a1:= a[r1,r]:= a1 - b0; w:= a[r,r]:= 1 / (a1 × b0);
for j:= 1 step 1 until r1 do b[j]:= (tammat(1, j, j, r,
a, a) + matmat(j + 1, r1, j, r, a, a)) × w;
elmveccol(1, r1, r, b, a, tamvec(1, r1, r, a, b) × w ×
.5);
for j:= 1 step 1 until r1 do
begin elmcol(1, j, j, r, a, a, b[j]);
elmcolvec(1, j, j, a, b, a[j,r])
end;
b[r1]:= b0
end;
r:= r1
end;
d[1]:= a[1,1]; a[1,1]:= 1; b[n]:= bb[n]:= 0
end tfmsymtri2;

comment mca 2320;

integer procedure qrivalsymtri(d, bb, n, em); value n; integer n;
array d, bb, em;
begin integer count, max, n1, k, k1;
real tol, tol2, macheps2, bbmax, t, r, w, c, s, shift, u, g, t2,
w2, p2, dk, cos2, sin2, oldcos2;
tol:= em[2] × em[1]; tol2:= tol × tol; macheps2:= em[0] ↑ 2;
count:= 0; bbmax:= 0; max:= em[4];
in: k:= n; n1:= n - 1;
next: k:= k - 1; if k > 0 then
begin if bb[k] ≥ tol2 then goto next;
if bb[k] > bbmax then bbmax:= bb[k]
end;
if k = n1 then n:= n1 else
begin
twoby2: t:= d[n] - d[n1]; r:= bb[n1]; if k < n - 2 then
begin w:= bb[n - 2]; c:= t × t; s:= r / (c + w);
g:= (w + s × c) × s; if g < tol2 then
begin n:= n - 1; n1:= n - 1; if g > bbmax then bbmax:= g;
goto twoby2
end
end negligible bb;
if abs(t) < tol then s:= sqrt(r) else
begin w:= 2 / t; s:= w × r / (sqrt(w × w × r + 1) + 1) end;
if k = n - 2 then
begin d[n]:= d[n] + s; d[n1]:= d[n1] - s; n:= n - 2 end
else
begin count:= count + 1; if count > max then goto end;
shift:= d[n] + s; if abs(t) < tol then
begin w:= d[n1] - s;
if abs(w) < abs(shift) then shift:= w
end;
k:= k + 1; g:= d[k] - shift; t2:= g × g; w2:= bb[k];
p2:= t2 + w2; oldcos2:= 1;
for k:= k + 1 step 1 until n do
begin k1:= k - 1; sin2:= w2 / p2; cos2:= t2 / p2;
dk:= d[k] - shift; u:= (g + dk) × sin2;
d[k1]:= g + u + shift; g:= dk - u;
t2:= if cos2 < macheps2 then w2 × oldcos2 else g × g
/ cos2; w2:= bb[k]; p2:= w2 + t2; bb[k1]:= sin2 × p2;
oldcos2:= cos2
end;
d[n]:= g + shift
end sweep
end;
if n > 0 then goto in;
end: em[3]:= sqrt(bbmax); em[5]:= count; qrivalsymtri:= n
end qrivalsymtri;

comment mca2322;

integer procedure qrivalsym2(a, n, val, em); value n; integer n;
array a, val, em;
begin array b, bb[1:n];
tfmsymtri2(a, n, val, b, bb, em);
qrivalsym2:= qrivalsymtri(val, bb, n, em)
end qrivalsym2;

comment Van een matrix a worden de elementen aldus bepaald: a[i,
j]:= SUM(k, 1, n, ar(i, k, n) × labda[k] × ar(k, j, n))/n. De
eigenwaarden van a zijn de elementen van een willekeurig te vullen
array labda[1 : n] en de elementen van de matrix van
eigenvectoren X zijn aldus bepaald: X[i, j] = ar(i, j, n)/sqrt(n).
n mag alleen een macht van 2 zijn, bij voorkeur een macht van 4;

integer procedure ar(i, j, n); value i, j, n; integer i, j, n;
comment ar(i, j, n) krijgt het teken van eigenvectorelement X[i, j];
begin integer sig;
sig:= 1; goto in;
dev: if i > n ∧ j > n then sig:= - sig; if n = 1 then goto end;
i:= REMAINDER(i, n); j:= REMAINDER(j, n); if i = 0 then i:= n;
if j = 0 then j:= n;
in: n:= n / 2; goto dev;
end: ar:= sig
end ar;

comment Procedure om een symmetrische matrix te genereren met
eigenwaarde mu[1] t/m mu[n], de procedure gebruikt integer
procedure ar(i, j, n);

procedure symtestmx(a, mu, n); value n; array a, mu; integer n;
begin integer i, j, k;
real s;

integer procedure col(i, j, n); value i, j, n; integer i, j, n;
begin integer k;
k:= j;
down: n:= n / 2; if i > n then
begin i:= i - n; j:= if k > n then j - n else j + n end;
if k > n then k:= k - n; if i > 1 then goto down; col:= j
end col;

for j:= 1 step 1 until n do
begin s:= 0;
for k:= 1 step 1 until n do s:= ar(1, k, n) × ar(k, j, n) ×
mu[k] / n + s; a[1,j]:= s; a[j,j]:= a[1,1]
end;
for j:= 3 step 1 until n do
for i:= j - 1 step - 1 until 2 do a[i,j]:= a[1,col(i,j,n)]
end symtestmx;

procedure primat(A, labda, n, numval, p, s); value n, numval;
integer n, numval, p; array A, labda; string s;
begin integer i, j, a, b;
NLCR; PRINTTEXT(s); a:= 1; b:= (if numval > 6 then 6 else numval);
NLCR;
AA: NLCR; NLCR; NLCR; if p = 1 then
begin for j:= a step 1 until b do print(labda[j]); NLCR end;
for i:= 1 step 1 until n do
begin NLCR;
for j:= a step 1 until b do print(A[i,j])
end;
if b ≠ numval then
begin a:= b + 1; b:= (if (b + 6) > numval then numval else b + 6);
goto AA
end;
NLCR
end primat;

start: NLCR; PRINTTEXT(
≮ mach. precisie qr precisie max. iter≯); NLCR;
for i:= 0, 2, 4 do
begin eps:= em[i]:= read; PRINT(eps) end;
n:= read;
next: begin integer l, j;
real t, sqinvn;
array ar, vec[1:n,1:n], val[1:n];
NLCR; PRINTTEXT(≮order≯); ABSFIXT(2, 0, n); NLCR;
for i:= 1 step 1 until n do val[i]:= read; symtestmx(vec, val, n);
for j:= 1 step 1 until n do
for i:= j step 1 until n do ar[i,j]:= ar[j,i]:= vec[j,i];
sqinvn:= 1 / sqrt(n); primat(ar, val, n, n, 0, ≮input≯); NLCR;
NLCR; PRINTTEXT(≮time≯); t:= time;
l:= qrivalsym2(vec, n, val, em); ABSFIXT(3, 2, time - t); NLCR;
NLCR; PRINTTEXT(≮number of calculated eigenvalues =≯);
ABSFIXT(2, 0, n - l); NLCR; NLCR; PRINTTEXT(≮eigenvalues are≯);
NLCR;
for i:= 1 step 1 until n do
begin PRINT(val[i]); NLCR end;
NLCR; PRINTTEXT(
≮ norm max negl. codiag. element number of qr steps≯
); NLCR;
for i:= 1, 3, 5 do PRINT(em[i]); NLCR; n:= read; if n ≠ 0 then
begin NEWPAGE; goto next end
end
end