MARST-compatible Algol-60 example
Expectedly, the
MARST translator, when fed with the Algol text below, complains about non-declared identifiers.
Capitalized identifiers like PRINTTEXT, NLCR, REMAINDER were defined in the X8 system environment; 'tammat', 'matmat', 'tamvec', 'count', 'qrivalsymtri' etc, were defined in a subroutine library which had to be loaded separately (also from papertape, probably as pre-compiled code).
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