DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »extorthof«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦58ca399f1⟧ »extbib« 
            └─⟦this⟧ 

TextFile

orthof=algol index.no list.yes

external
real procedure orthof(a,A,d,N,n,nf);
value N,n,nf; integer N,n,nf; array a,A,d;
begin
comment The procedure solves the least squares equation a x = b
(N linear equations with n unknowns, N>n) by orthogonalizing the
columns of the array a. b is stored in column no. n+1 of
a(1:N,1:n+1). Also the covariance matrix (i.e. the inverse of
the least squares matrix) is calculated in the first n columns
of the array A(1:n,1:n+1), while the solution x is stored in
column no. n+1 of A. In d(1:nf) an estimate of the significance
of the fixed contants is produced;
integer i,j,k,m,mt,nt;
real p,q,length,gramdet,rnorm;
m:= n+1; nt:= n+nf; mt:= nt+1; rnorm:= 0;  gramdet:= 1;
for i:=1 step 1 until N do rnorm:= rnorm+a(i,mt)**2;
rnorm:= 2*sqrt((N-n)*rnorm);
for j:=1   step  1 until nt do begin
length:= 0; for i:=1 step 1 until N do
length:= length + a(i,j)**2;
for k:=(if j>n then n else j-1) step -1 until 1 do begin q:= A(k,k);
if q>0 then begin
   p:= 0; for i:=1 step 1 until N do
   p:= p + a(i,j)*a(i,k);
   if j<=n then begin A(k,j):= p;  p:= A(j,k):= q*p;
       end else p:= p*q;
   for i:=1 step 1 until N do a(i,j):= a(i,j)-a(i,k)*p
end end k, orthogonalisation;
p:= q:= 0;
for i:=1 step 1 until N do begin
   p:= p + a(i,j)**2;
   q:= q + a(i,j)*a(i,mt)
end i;
if j>n then
   d(j-n):= if p>length*'-6 then abs q/(sqrt(p)*rnorm) else 0
else if p<length*'-8 then begin
   for k:=1 step 1 until n do A(j,k):= A(k,j):= 0;
   A(j,m):= 0
end else begin
   A(j,j):= 1/p; A(j,m):= q; gramdet:= gramdet*p/length
end end j;
for j:=n step -1 until 1 do begin q:= A(j,j);
if q>0 then begin
   p:= A(j,m); for k:=j+1 step 1 until n do
   p:= p - A(j,k)*A(k,m);
   A(j,m):= q*p
end end solution;
orthof:= gramdet;

for i:=n step -1 until 1 do begin p:= A(i,i);
if p<>0 then begin  m:= i-1;
for j:=i-1 step -1 until 1 do begin
   q:= -A(i,j); for k:=j+1 step 1 until m do
   q:= q-A(i,k)*A(k,j);
   A(i,j):= q; A(j,i):= q*p
end end end;

for i:=1 step 1 until n do if A(i,i)<>0 then begin
for j:=i step 1 until n do begin
   q:= A(i,j); for k:= j+1 step 1 until n do
   q:= q+A(i,k)*A(k,j);
   A(i,j):= A(j,i):= q
end end end;
end
▶EOF◀