|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 33792 (0x8400)
Types: TextFile
Notes: RCSL-44-RT-900
Names: »nllprtx«
└─⟦667bb35d6⟧ Bits:30007480 RC8000 Dump tape fra HCØ.
└─⟦4334b4c0b⟧
└─⟦this⟧ »nllprtx«
; gi nll procedures * page 1 5 09 77, 15.13;
; nll_pr
; ******
if listing.yes
char 10 12 10
nll_pr = set 1
nll_pr = algol
external integer procedure nll_pr
_________________________________
_ (z, t, exist);
zone z;
integer t;
boolean exist;
begin
integer q_proc;
real str;
nll_pr :=
q_proc := 10;
write(z, nl, 3, <:; :>, q_proc, <:_procedures:>, nl, 2,
_ <:; proc_name__________version:>, nl, 1);
for t := 1 step 1 until q_proc do
if exist then
begin
str := real (case t of (
<* 1 *> <:nll<95>pr:>,
<* 2 *> <:otn<95>gier:>,
<* 3 *> <:nll<95>gier:>,
<* 4 *> <:nll:>,
<* 5 *> <:nll1:>,
<* 6 *> <:nll<95>epu:>,
<* 7 *> <:otn<95>epu:>,
<* 8 *> <:epu<95>block:>,
<* 9 *> <:epu<95>print:>,
<* 10 *> <:nll<95>block:>,
<* *> <::>));
wr_date_time(z, proc_transla(string str)
+ 0*write(z, sp, 16 - write(z, nl, 1, string str, <:, :>)));
end;
end nll_pr;
end
if warning.yes
(mode 0.yes
message nll_pr not ok
lookup nll_pr)
\f
; gi nll procedures * page 2 5 09 77, 15.13;
; otn_gier
***********
if listing.yes
char 10 12 10
otn_gier = set 1
otngier = algol
external procedure otn_gier
__________________________
(B, NL, n);
value n;
array B, NL;
integer n;
comment
otn_gier (call) no type procedure
the procedure computes the contributions from
the product B**T times B to the triangular matrix
NL mapped as a linear array columnwize.
B (call) real array
a column of observation equations (for one observ.).
NL (call and return) real array
a triangular array of normals mapped on NL with
nl(1, 1)=NL(1), nl(1, 2)=NL(2), nl(2, 2)=NL(3), etc.
n (call) integer
the number of columns in B and the triagular NL;
begin integer I, r, c;
I:=0;
for r:=1 step 1 until n do
for c:=1 step 1 until r do
begin I:=I+1; NL(I):=NL(I)+B(r)*B(c) end
end otn_gier;
end
if warning.yes
(mode 0.yes
message otn_gier not ok
lookup otn_gier)
\f
; gi nll procedures * page 3 5 09 77, 15.13;
; nllgier
; *******
if listing.yes
char 10 12 10
nllgier = set 1 disc
nllgier = algol
external
real procedure nll(N, n, fr, fc, BS);
value n, fr, fc, BS; integer n, fr, fc; boolean BS; array N;
begin integer r, c, p, I, Ir, Ic; real sum;
for r:=fr+1 step 1 until n do
begin Ir:=r*(r-1)//2; I:=Ir+fc;
for c:=fc+1 step 1 until r do
begin sum:=0; Ic:=c*(c-1)//2; I:=I+1;
for p:=fc+1 step 1 until c-1 do
sum:=sum+N(Ir+p)*N(Ic+p);
N(I):=N(I)-sum; if r<>c then N(I):=N(I)/N(Ic+c)
end;
if r<>n then N(I):=sqrt(N(I))
else
begin nll:=N(I);
if BS then for c:=n-1 step -1 until 1 do
begin Ir:=I:=I-1; Ic:=c*(c+1)//2; N(I):=N(I)/N(Ic);
for p:=c-1 step -1 until 1 do
begin Ir:=Ir-1; Ic:=Ic-1;
N(Ir):=N(Ir)-N(I)*N(Ic)
end
end
end
end
end nll;
end
if warning.yes
(mode 0.yes
message nll_gier not ok
lookup nll_gier)
\f
; gi nll procedures * page 4 5 09 77, 15.13;
; nll
; ***
if listing.yes
char 10 12 10
nll = set 1 disc
nll = algol
external
real procedure nll(N, n, explim, maxlim, SING, KT);
value n, explim; integer n, explim;
array N; long array SING; real KT, maxlim;
begin integer r, t, tt, c, p, I, Ir, Ic, expNI, maxlim1, maxtest;
real sum, NI; boolean sing, BS;
tt:=0; t:=1; maxlim1:=-8338607;
explim:=explim shift 12;
BS:=n>0;
for r:=1 step 1 until n do
begin
I:=Ir:=r*(r-1)//2;
for c:=1 step 1 until r do
begin
sum:=0; Ic:=c*(c-1)//2; I:=I+1;
for p:=1 step 1 until c-1 do
sum:=sum-N(Ir+p)*N(Ic+p);
NI:=N(I);
if r<>c then N(I):=(NI+sum)/N(Ic+c)
else
if r<>n then
begin
tt:=tt+1; if tt=49 then begin t:=t+1; tt:=1; end;
SING(t):=SING(t) shift 1;
expNI:=(NI extract 12) shift 12; NI:=NI+sum;
maxtest:=expNI-((NI extract 12) shift 12);
sing:=maxtest>explim;
if -, sing & NI>0.0 then N(I):= sqrt(NI)
else begin N(I):='600; SING(t):=SING(t)+1; end;
if maxtest>maxlim1 then maxlim1:=maxtest;
end;
end;
end;
nll:=NI:=NI+sum;
KT:=mu*ln(-NI/sum)/2; maxlim:=maxlim1/4096;
if BS then for c:=n-1 step -1 until 1 do
begin
Ir:=I:=I-1; Ic:=c*(c+1)//2; N(I):=N(I)/N(Ic);
for p:=c-1 step -1 until 1 do
begin
Ir:=Ir-1; Ic:=Ic-1;
N(Ir):=N(Ir)-N(I)*N(Ic);
end;
end;
end nll;
end
if warning.yes
(mode 0.yes
message nll not ok
lookup nll)
\f
; gi nll procedures * page 5 5 09 77, 15.13;
; nll1
; ****
if listing.yes
char 10 12 10
nll1 = set 1 disc
nll1 = algol
external
real procedure nll(N, n, explim, maxloss, SING, KT);
value n, explim; integer n, explim, maxloss;
array N; long array SING; real KT;
begin integer r, t, tt, c, p, I, Ir, Ic, expNI, maxloss1, maxtest;
real sum, NI; boolean sing, BS;
tt:=0; t:=1; maxloss1:=0;
explim:=explim shift 12;
BS:=n>0; n := abs n;
for r:=1 step 1 until n do
begin
I:=Ir:=r*(r-1)//2;
for c:=1 step 1 until r do
begin
sum:=0; Ic:=c*(c-1)//2; I:=I+1;
for p:=1 step 1 until c-1 do
sum:=sum-N(Ir+p)*N(Ic+p);
NI:=N(I);
if r<>c then N(I):=(NI+sum)/N(Ic+c)
else
if r<>n then
begin
tt:=tt+1; if tt=49 then begin t:=t+1; tt:=1; end;
SING(t):=SING(t) shift 1;
expNI:=(NI extract 12) shift 12; NI:=NI+sum;
maxtest:=expNI-((NI extract 12) shift 12);
sing:=maxtest>explim;
if -,sing and NI>0.0 then N(I):= sqrt(NI)
else begin N(I):='600; SING(t):=SING(t)+1; end;
if maxtest>maxloss1 then maxloss1:=maxtest;
end;
end;
end;
comment TEST: :
write(out, <:<10>nll1:NI,sum= :>,
<< -d.ddd'-d>, NI, sum);
nll:=NI:=NI+sum;
KT:= if -sum>'-300 and NI>0 then mu*ln(-NI/sum)/2 else '-50;
maxloss := maxloss1 shift (-12);
if BS then for c:=n-1 step -1 until 1 do
begin
Ir:=I:=I-1; Ic:=c*(c+1)//2; N(I):=N(I)/N(Ic);
for p:=c-1 step -1 until 1 do
begin
Ir:=Ir-1; Ic:=Ic-1;
N(Ir):=N(Ir)-N(I)*N(Ic);
end;
end;
end nll;
end
if warning.yes
(mode 0.yes
message nll1 not ok
lookup nll1)
\f
; gi nll procedures * page 6 5 09 77, 15.13;
;
; n l l e p u
; ____________
;
if listing.yes
char 10 12 10
nllepu=set 36
nllepu=algol message.no
external
real procedure NLL(N, n, explim, maxlim, SING, KT);
___________________________________________________
value n, explim;
integer n, explim, maxlim;
array N;
long array SING;
real KT;
begin
comment
Solution of a system of linear equations with
positiv definit coefficient matrix. Choleski
decomposition and the aritmetic unit epu are
used.
Parameters.
NLL (real, return value). The pvv-sum
_ minus the square-sum of the inter-
_ mediate right side. The pvv-sum is
_ the square - sum of the original
_ right side.
N (real array, call and return value)
_ At call the coefficients and the
_ right side. Only the upper half is
_ stored, which is done column by
_ column, i.e. (1, 1), (1, 2), (2, 2),
_ (1, 3), (2, 3), (3, 3), (1, 4), ...,
_ (4, 4), ..., (1, n), ..., (n, n) is
_ stored in N(1) to N(n*(n+1)/2). The
_ last column is the right side and
_ the very last element (which is not
_ changed) ougth to be the pvv-sum.At
_ return the Choleski factor and the
_ solution is stored in the same way.
_ Double precision is used at call-
_ time. All the tails are stored in
_ the last part of N, that is from
_ element (n*(n+1)/2 + 1) to element
_ (n*(n+1)). At return single preci-
_ sion is used, i.e. the tails are un-
_ defined, exept the diagonal elements
_ which holds the status word and the
_ exponent loss from the epu-diagonal-
_ operation.
;
\f
comment gi nll procedures * page 7 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
comment
n (integer, call value). The order of
_ the system plus one. If sign is ne-
_ gativ the absolute value is used
_ and no back-substition is preformed.
explim (integer, call value). Maximum al-
_ lowed loss of digits in the subtrac-
_ tion which forms the square of the
_ new diagonal element. A value of 30
_ is a reasonable standard.
maxlim (integer, return value). The actual
_ maximum digit loss.
SING (long array, return value). Bit pat-
_ tern giving the singularities of
_ the matrix.
KT (real, return value). The Krarup
_ number 0.5*log(NLL/pvv).
External referances.
integer procedure abs_addr.
real procedure sub_double.
long prvers.
procedure epu_block (GI nr 75004).
process ympe.
GI program index nr 75006.
January 1975, WLW.
updated feb 77, WLW.
- -> o 0 o <- -
;
\f
comment gi nll procedures * page 8 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
integer r, t, tt, c, p, I, Ir, Ic, expNI, maxlim1,
_ add1, div1, str1, mla1,
_ zmls, stp, zadd1, mls, dia,
_ max_test, addr_sum, addr_sum1,
_ baseN, baseH, last_p, code_length,
_ n1;
real sum, sum1, NI, op_code;
zone epu((3*6*100+3)//4, 1, epu_block);
integer array field word;
boolean field byte;
integer field ord;
boolean sing, BS, epu_test, bs_test;
procedure code_dump;
begin
write(out, <:<10>epu code dump.:>);
for byte:= 1 step 1 until word do
write(out, false add 10,
if byte mod 6 = 1 then 1 else 0,
<< dddd>, epu.byte extract 12);
for ord:= 2 step 2 until word do
begin
byte:= ord -1;
if ord mod 6 = 2 then
begin
write(out, false add 10, 1,
if (epu.byte shift (-11) extract 1) = 1
then <:z :> else <: :>);
op_code:= real (
if (epu.byte extract 11 > 7)
then <:ill:> else ( case (epu.byte extract 7+1) of
(<:add:>, <:sub:>, <:mla:>, <:mls:>, <:str:>,
_<:div:>, <:sqr:>, <:dia:>, <:xxx:>)));
write(out, string op_code,
<< dddd>, epu.ord extract 12 <*count*>);
end
else write(out,
<<-ddd ddd ddd>,
(if (ord mod 6 = 4 and epu.ord <> addr_sum )
then ((epu.ord - baseN) // 4)
else (if (ord mod 6 = 0 and epu.ord <> addr_sum1)
_ then ((epu.ord - baseH) // 4)
_ else epu.ord)));
end for ord;
end epu_dump;
BS := n>0;
n := abs(n);
n1 := n - 1;
prvers:= 5 09 77 15 13;
comment index check in N;
r:= system(3)array_bounds:(c, N);
if r>1 or c<n*(n+1) then
begin
write(out, <:<10>***Coefficient matrix:>,
_ <:<10> in call of nllepu:>,
_ <:<10> is too small.:>,
_ <:<10> Bad bounds::>, r, c);
system(9)run_time_alarm:(n, <:<10>order___:>);
end if r>1;
word:= 0;
code_length:= 3*6*100;
open(epu, 14, <:ympe:>, 31 shift 18);
out_rec6(epu, code_length);
\f
comment gi nll procedures * page 9 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
add1 := 0 shift 23 + 1;
zadd1:= 1 shift 23 + 0 shift 12 + 1;
mla1 := 2 shift 12 + 1;
mls := 3 shift 12;
zmls := 1 shift 23 + 3 shift 12;
str1 := 4 shift 12 + 1;
div1 := 5 shift 12 + 1;
dia := 7 shift 12 + explim extract 12;
stp := 4095 shift 12 + 1;
addr_sum := abs_addr(sum);
addr_sum1:= abs_addr(sum1);
baseN := abs_addr(N);
baseH := baseN + (n*(n+1)*2);
epu_test := false;
bs_test := false;
if epu_test then
write(out, <:<10>adresser N, H:>, baseN, baseH,
<:<10>sum, sum1:>, addr_sum, addr_sum1);
for r:=1 step 1 until n do
begin
I:=Ir:=r*(r-1)*2;
for c:=1 step 1 until r do
begin
I:= I+4;
Ic:= c*(c-1)*2;
comment generate epu-code for
N(c, r) - summa (N(p, r)*N(p, c)),
p= 1, ..., c-1;
if epu_test then write(out, <:<10>mls r, c, word:>, r, c, word);
lastp:= if c=1 then 3 else 6;
for p:= 1 step 1 until lastp do
epu.word(p):= case p+6-lastp of (
zmls + (c-1),
baseN + 4 + Ir,
baseN + 4 + Ic,
(if c=1 then 1 shift 23 else 0) + add1,
baseN + I,
baseH + I);
word:= word + 2*last_p;
comment if c=1 the mls will work as a no-op, therefor
it is omitted and the clearing of the epu-accumulator
R is made by the 1 shift 23 in the add1.
;
\f
comment gi nll procedures * page 10 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
if word+24 >= code_length then
begin
comment start and wait for the epu-calculations. As
the epu-unit can't be reserved the contens of the
accumulator R must be saved before stop of the epu;
for p:= 1 step 1 until 6 do
epu.word(p):= case p of (
str1,
addr_sum,
addr_sum1,
stp,
0,
0);
if epu_test then
write(out, <:<10>fuld<10>:>);
word:= word + 12;
changerec6(epu, word);
outrec6(epu, code_length);
if epu_test then code_dump;
if epu_test then
for p:= 1 step 1 until I//4 do
write(out, nl, 1, <<dd>, p, << -d.ddd ddd'+d>, N(p));
word:= 0;
comment reload R from sum;
epu.word(1):= zadd1;
epu.word(2):= addr_sum;
epu.word(3):= addr_sum1;
word:= word+6;
end if word+18>= code_length;
\f
comment gi nll procedures * page 11 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
if r<>c then
begin
comment non-diagonal element;
epu.word(1):= div1;
epu.word(2) := baseN + I;
epu.word(3) := baseN + Ic + 4*c;
word:= word + 6;
comment Pay attention to the fact that both divisor and
the result are single precision. The divisor is a diagonal
element , i.e. the result of a sqrt . This operation must
therefor deliver a single precision result;
end non-diag
else
begin
comment diagonal element;
comment TEST: :
write(out, <:<10>nll epu dia::>, I);
if r = n then
begin
comment diagonal element of right hand side.
Start and wait for epu with save of sum;
for p:= 1 step 1 until 6 do
epu.word(p):= case p of (
str1,
addr_sum,
addr_sum1,
stp,
0,
0);
word:= word + 12;
comment start and wait for the epu-calculation;
if epu_test then write(out, <:<10>diag<10>:>, r, c, <:<10>:>);
if epu_test then
for p:= 1 step 1 until n*(n+1)//2 do
write(out, <:<10>:>, <<dd>, p, << -d.ddd ddd'+d>, N(p),
N(p + n*(n+1)//2) );
changerec6(epu, word);
if eputest then write(out, <:<10>changerec done, word:>, word);
if epu_test then code_dump;
outrec6(epu, code_length);
word:= 0;
end r = n
\f
<* gi nll procedures * page 12 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 *>
else
begin
comment computation of square root;
epu.word(1):= dia;
epu.word(2):= baseN + I <* sqrt *>;
epu.word(3):= baseH + I <* status shift 24 + exp_loss *>;
word := word + 6;
end if r<>n;
end diagonal element;
end for c, i.e. row no r clompleted;
end;
\f
comment gi nll procedures * page 13 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
comment the triangular choleski-factor is now
stored in N and the forward substitutio is
preformed. The new intermediate rigth sides,
which is the result of the forwart substitu-
tion, is stored in the last part of N;
if epu_test then
write(out, <:<10>Forlængs slut, I, sum, N(I):>,
<< -d.dd'd>, I, sum, N(I//4));
NLL:= NI:= sum + sum1 <* rounding, nll:= pvv - pww *>;
sub_double(sum, sum1, N(I//4), N(I//4 + n*(n+1)/2), sum, sum1);
comment TEST: :
write(out, <:<10>pww = :>, << -d.ddd'-ddd>,
sum, sum1, sum + sum1);
sum1:= sum + sum1 <* rounding, sum1 := -pww *>;
if abs(sum1) < '-100 then KT := - '300 else
if (-NI/sum1) < '-10 then KT := - '300 else
KT:=mu*ln(-NI/sum1)/2;
comment test of exponent losses;
max_lim1:= 0;
tt := 0 <* bit pointer *>;
t := 1 <* double word pointer *>;
c := n * (n + 1) // 2 <* base index of tail *>;
for p:= 1 step 1 until n1 do
begin
tt:= tt + 1;
if tt = 49 then
begin
tt:= 1;
t := t + 1;
end;
SING(t):= SING(t) shift 1;
r:= c + p * (p + 1) // 2 <* index of tail of diagonal element *>;
NI:= N(r);
if NI shift (-24) extract 24 <* status *> > 1 then
SING(t):= SING(t) add 1;
if NI extract 24 > max_lim1 then
max_lim1:= NI extract 24;
end for p;
max_lim:= max_lim_1;
\f
comment gi nll procedures * page 14 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
comment clear 2. part of rigth side;
c:= n*(n+1)//2 + n*(n-1)//2;
for p:= n-1 step -1 until 1 do
N(c+p):= 0;
if bs_test then
begin
write(out, <:<10>nll I:>, I);
write(out, <:<10>nll højre<10>:>);
for c:= n-1 step -1 until 1 do
write(out, false add 10,
if c mod 4 = 0 then 1 else 0,
<< -d.ddd ddd'+dd>, N(n*(n-1)//2 + c) );
end if bs_test;
if BS then
for c:= n-1 step -1 until 1 do
begin
comment compute solution no c and subtract
the effect on the rigth side;
Ir:= I:= I-4;
Ic:= c*(c+1)*2;
comment generate code for final solution
for unknown no c, x(c):= x(c)/diagonal(c);
if bs_test then
write(out, <:<10>nll div:>, << -d.ddd ddd'+dd>, N(I//4), N(Ic//4));
for p:= 1 step 1 until 6 do
epu.word(p):= case p of (
zadd1,
baseN + I,
baseH + I,
div1,
baseN + I,
baseN + Ic);
word:= word + 12;
\f
comment gi nll procedures * page 15 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
for r:= c-1 step -1 until 1 do
begin
comment subtract the influence of unknown no c
from the upper c-1 elements of the rigth side;
Ir:= Ir-4;
Ic:= Ic-4;
if word + 36 >= code_length then
begin comment start and wait for the epu;
epu.word(1):= stp;
epu.word(2):= 0;
epu.word(3):= 0;
word:= word + 6;
changerec6(epu, word);
out_rec6(epu, code_length);
word:= 0;
end if word + 30;
for p:= 1 step 1 until 9 do
epu.word(p):= case p of (
zadd1,
baseN + Ir,
baseH + Ir,
mls + 1,
baseN + I,
baseN + Ic,
str1,
baseN + Ir,
baseH + Ir);
word:= word + 18;
end for r;
end for c;
epu.word(1):= stp;
epu.word(2):= 0;
epu.word(3):= 0;
word:= word + 6;
changerec6(epu, word);
close(epu, true);
if bs_test then
begin
write(out, false add 10, 2);
for c:= n-1 step -1 until 1 do
write(out, false add 10,
if c mod 4 = 0 then 1 else 0,
<< -d.ddd ddd'+dd>, N(n*(n-1)//2 + c));
end if bs_test;
end NLL;
end
if warning.yes
(mode 0.yes
message nll_epu not ok
lookup nll_epu)
\f
; gi nll procedures * page 16 5 09 77, 15.13;
if listing.yes
char 10 12 10
otnepu=set 1
otnepu=algol
external procedure otn_epu(B, NL, n);
_____________________________________
value n;
integer n;
array B, NL;
begin
comment
The procedure add the contribution to the
normal equations from a single observation
equation. The arithemetic unit epu is used,
which give double precision results.
Parameters.
B (real array, call value). One row in the
_ observation equations, stored in the
_ elements B(1), ..., B(n).
NL (real array, call and return value). The
_ dobbel precision normal equations, in-
_ cluding the right-hand side. Only the
_ upper triangular part is stored, which
_ is done column by column. In the ele-
_ ments NL(1) to NL(n*(n+1)/2) the high
_ part of the coefficients are stored,
_ and in the last n*(n+1)/2 elements the
_ low part of the coefficients are stored.
n (integer, call value). The number of un-
_ knowns plus one.
External refrences.
integer procedure abs_addr
long prvers
procedure epu_block
process ympe
GI-program no 75005,
WW, jan. 1975.
;
\f
comment gi nll procedures * page 17 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
comment declarations;
integer z_add1, mla_1, str_1, stp_1,
_ head, tail, base, r, c, code_length;
integer array field word;
zone epu((100*18 + 6 + 3)//4, 1, epu_block);
prvers:= 5 09 77 15 13;
comment address constants for epu-calculations;
base:= abs_addr(B);
head:= abs_addr(NL);
tail:= head + (n*(n+1))*2;
n := 4 * n;
comment initialize epu-instruction nmenonics;
z_add_1:= 1 shift 23 + 0 shift 12 + 1;
_ mla_1:= 2 shift 12 + 1;
_ str_1:= 4 shift 12 + 1;
_ stp_1:= 4095 shift 12 + 1;
comment initialize the zone epu
and the associated variables;
word:= 0;
code_length:= 100*18;
open(epu, 14, <:ympe:>, 63 shift 17);
out_rec_6(epu, code_length + 6);
\f
comment gi nll procedures * page 18 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
for r:= 4 step 4 until n do
for c:= 4 step 4 until r do
begin
comment generate code for calculation of
the influence on element (c, r) in the
normal equations;
head:= head + 4;
tail:= tail + 4;
if word = code_length then
begin
comment empty the code-buffer;
epu.word(1):= stp_1;
epu.word(2):=
epu.word(3):= 0;
out_rec_6(epu, code_length + 6);
word:= 0;
end;
epu.word(1):= z_add_1;
epu.word(2):= head;
epu.word(3):= tail;
epu.word(4):= mla_1;
epu.word(5):= base + r;
epu.word(6):= base + c;
epu.word(7):= str_1;
epu.word(8):= head;
epu.word(9):= tail;
word:= word + 18;
end for r and c;
epu.word(1):= stp_1;
epu.word(2):=
epu.word(3):= 0;
word:= word + 6;
change_rec_6(epu, word);
close(epu, true);
end otn_epu;
end
if warning.yes
(mode 0.yes
message otn_epu not ok
lookup otn_epu)
\f
; gi nll procedures * page 19 5 09 77, 15.13;
if listing.yes
char 10 12 10
epublock=set 1
epublock=algol list.no
external procedure epublock(z, s, b);
_____________________________________
zone z;
integer s, b;
begin
comment
Block procedure for the arithmetic unit epu.
The status bit 6 = 1 shift 17 = underflow
is handeled by applying the z-modif to the
following instruction and restarting
the unit.
If the procedure is called with b < 0 it will
return imediately with b=no of underflows
and this underflow counter will
be cleared.
GI program no 75004,
WW, jan 1975.
;
integer bit, share, record_base,
_ result, buffer_address,
_ new_status;
integer array descriptor(1:20),
_ answer(1:8);
own integer underflows;
boolean first, test, old_z;
boolean field byte;
\f
comment gi nll procedures * page 20 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
if b >= 0 then
begin
comment block procedure part;
first:= true;
test := false;
byte := 1;
getzone6(z, descriptor);
record_base:= descriptor(14);
share:= descriptor(17);
getshare6(z, descriptor, share);
if test then
begin
write(out, <:<10>epu block proc:>, s, s extract 12, b,
<:<10>share no :>, share,
<:<10>record base :>, record_base,
<:<10>share state:>, descriptor(1),
<:<10>first shared:>, descriptor(2),
<:<10>last shared:>, descriptor(3),
<:<10>first instru:>, descriptor(5),
<:<10>last instru:>, descriptor(6),
<:<10>R1, R2 :>, << -d.ddd ddd ddd'-dd>,
descriptor(7), descriptor(8));
end if test;
for bit:= 16 step 1 until 22 do
if s shift (-bit) extract 1 = 1 then
begin
comment status bit no bit is on;
if bit <> 17 then
begin
write(out, if first then
<:<10>***epu status::> else
<:<10>______________:>,
case bit-15 of (
<:16:>,
<:underflow:>,
<:overflow:>,
<:address fault:>,
<:20:>,
<:21:>,
<:illegal instruction:>));
first:= false;
end bit <>17;
\f
comment gi nll procedures * page 21 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
case bit - 15 of
begin
;
begin
comment 1 shift 17: underflow;
if byte > 1 and test then
write(out, <:<10>epu-extra:>, byte, b);
byte:= b + byte;
underflows:= underflows + 1;
if test then write(out, <:<10>uf:>, byte, z.byte extract 12);
comment modify the operation and start for the rest;
old_z := z.byte;
z.byte:= false add (old_z extract 11) add (1 shift 11);
descriptor(5):= descriptor(5) + b;
descriptor(4):= 5 shift 12;
setshare6(z, descriptor, share);
buffer_address:= monitor(16)send_message:(z, share, descriptor);
if buffer_address = 0 then
system(9)run_time_alarm:(0, <:<10>buf_claim:>);
result:= monitor(18)wait_answer:(z, share, answer);
if result <> 1 then
system(9)run_time_alarm:(result, <:<10>epu-ans.:>);
comment reset the modified instruction and
the status word;
z.byte:= old_z;
new_status:= answer(1);
b:= answer(2);
s:= logand(exor(-1, 1 shift 17 + 1 shift 8 + 1), s) extract 24;
s:= logor(s, new_status shift (-12) shift 12) extract 24;
bit:= 16;
if test then write(out, <:<10>***epu-new status::>,
new_status shift (-12), new_status extract 12,
s shift (-12), s extract 12);
end case underflow;
;
;
;
;
;
end case bit of;
end if s;
\f
comment gi nll procedures * page 22 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
if -, first then write(out, <:<10>:>);
if s <> 1 shift 1 then stderror(z, s, b);
end block-procedure part
else
begin
comment return the underflowcounter;
b:= underflows;
underflows:= 0;
end;
end;
end
if warning.yes
(mode 0.yes
message epu_block not ok
lookup epu_block)
\f
; gi nll procedures * page 23 5 09 77, 15.13;
; epu_print
; *********
if listing.yes
char 10 12 10
epuprint=set 1 disc
epuprint=algol
external procedure epu_print
____________________________
_ (epu, code_length);
value code_length;
zone epu;
integer code_length;
begin
integer instr;
boolean z;
integer array field word;
codelength:= (codelength//6)*6 - 1;
write(out, <:<10><10>epu-instruktioner<10>:>);
for word := 0 step 6 until codelength do
begin
z := epu.word(1)<0;
write(out, false add 10, if z then 2 else 1);
instr := epu.word(1) shift (-12) extract 11;
if z and instr=2047 then
write(out, <:__stp:>)
else
begin
write(out, if z then <:z:> else <:_:>);
if instr < 9 then
write(out, <:_:>,
case instr + 1 of (
<:add:>, <:sub:>, <:mla:>, <:mls:>, <:str:>,
<:div:>, <:sqr:>, <:dia:>, <:chl:>))
else write(out, <<dddd>, instr);
end;
write(out, <:__:>, <<dddd>, epu.word(1) extract 12);
z := z and instr=2047;
if -, z then
write(out, <:, ___:>, <<-d_ddd_ddd>, epu.word(2));
if instr <> 6 and -, z then
write(out, <:, ___:>, <<-d_ddd_ddd>, epu.word(3));
end word-loop;
write(out, false add 10, 2);
end;
end
if warning.yes
(mode 0.yes
message epu_print not ok
lookup epu_print)
\f
; gi nll procedures * page 24 5 09 77, 15.13;
; nll_block
; *********
if listing.yes
char 12 10
nllblock=set 1
nllblock=algol list.no
external procedure nllblock(z, s, b);
_____________________________________
zone z;
integer s, b;
begin
comment
Simple block procedure for the arithmetic unit epu (RC 4090).
The procedure write on current output the name (ref a, page 11)
or the no (algol equivalent, ref b, p 6-11) of the 12 hardware
status bits of the epu status-word s.
References:
(a) P. Koch Andersson:
Functional Description for the RC4090,
External Processing Unit.
RCSL 44-RT 900.
(b) Hans Dinesen Hansen:
Algol 6 User's Manual.
RCSL 31-D 322.
nll_pr procedure lib no 10.
Willy Lehmann Weng, nov 1976.
;
\f
comment gi nll procedures * page 25 5 09 77, 15.13
0 1 2 3 4 5 6 7 8 9 ;
integer bit;
boolean first;
first:= true;
for bit:= 12 step 1 until 23 do
if s shift (-bit) extract 1 = 1 then
begin
comment status bit no bit is on;
write(out, if first then
<:<10>***epu status::> else
<:<10>______________:>,
case bit-11 of (
<:12:>,
<:13:>,
<:14:>,
<:15:>,
<:16:>,
<*17*> <:underflow:>,
<*18*> <:overflow:>,
<*19*> <:address fault:>,
<:20:>,
<*21*> <:time out:>,
<*22*> <:illegal instruction:>,
<:23:>));
first:= false;
end if s;
if -, first then
begin
write(out, <:<10>:>);
epu_print(z, b);
end;
if s <> 1 shift 1 then
begin
case epumode+1 of
begin
stderror(z, s, b);
system(9)alarm:(s, <:nllblock:>);
system(9)alarm:(s, <:nllblock:>);
end;
end;
end;
end
if warning.yes
(
mode 0.yes
message syntax fejl i nllblock
)
\f
; gi nll procedures * page 26 5 09 77, 15.13;
if 0.no
(
if listing.yes
char 10 12 10
nll_pr = compresslib ,
otn_gier,
nll_gier,
nll,
nll1,
nll_epu,
otn_epu,
epu_block,
epu_print,
nll_block,
if 2.no
scope project ,
nll_pr,
otn_gier,
nll_gier,
nll,
nll1,
nll_epu,
otn_epu,
epu_block,
epu_print,
nll_block,
lookup ,
nll_pr,
otn_gier,
nll_gier,
nll,
nll1,
nll_epu,
otn_epu,
epu_block,
epu_print,
nll_block,
)
if 10.no
end
if 10.no
finis
\f
; gi nll procedures * page 27 5 09 77, 15.13;
; test of nll_1 and nll_epu
; *************************
if 0.yes
end
if 0.yes
finis
mode 0.no ; clear error flag
pie_test = set 1
pie_test = algol list.no message.no
begin
comment
Test progran for the procedures nll_1 and nll_epu.
The matrix of case 4 from appendix C of Westlake
'A Handbook of Numerical Matrix Inversion
and Solutions of Linear Equations'
is used as test matrix.
The order of the matrix is read from current input.
A negativ value means single precision solution and
a positiv value means double precison sulution.
This order reading is repeated until order zero
is read.
Willy Weng, jan 1977.
;
integer i, j, n1, n, max, exp_lim, max_los;
real dia, non_dia, kt;
boolean double;
exp_lim:= 30;
for i:= read(in, n) while n <> 0 do
begin
double:= n > 0;
n := abs(n);
n1 := n - 1;
max := if double then n*(n+1) else n*(n+1)/2;
begin
long array sing(1:(n+47)//48);
real array A(1:max);
write(out, nl, 2, <:Pie-test of order :>, n1);
dia := 2/(2*n1-1);
non_dia:= - 1/((n1 - 1)*(2*n1 - 1));
for i := 1 step 1 until n1 do
for j := 1 step 1 until i do
A(i*(i-1)/2 + j) :=
if i=j then dia else non_dia;
i := n * (n - 1) / 2 + 1;
A(i) := 1;
for i := i + 1 step 1 until max do
A(i) := 0;
A(n*(n+1)/2) := 1;
if double then
write(out, <:<10>Double nll = :>,
nll_epu(A, n, exp_lim, max_los, sing, kt))
else
write(out, <:<10>Single nll = :>,
nll_1(A, n, exp_lim, max_los, sing, kt));
write(out, <:<10>kt = :>, kt,
_ <:<10>los= :>, max_los,
_ <:<10>:>);
max := n*(n+1)/2;
j := n*(n-1)/2 + 1;
for i:= j step 1 until max do
write(out, << -d.ddd'-dd>, nl,
if (i-j) mod 5 = 0 then 1 else 0,
if i=j then ((A(i) - n1) / n1) else (A(i) - 1));
end array decl block;
end read block;
end
if warning.yes
(mode 0.yes
o c
message pie test not ok
end)
if 0.yes
finis
pie_test
4 10 -4 -10 0
; 3 ; 10 ; 20 ; 50 ;
; -3 ; -10 ; -20 ; -50 ; 0 ;
end
finis
▶EOF◀