DataMuseum.dk

top - metrics - download

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

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

TextFile

ortho=algol list.yes index.no

external
real procedure ortho(a,A,N,n);
value N,n; integer N,n; array a,A;
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;
integer i,j,k,m;
real p,q,length,gramdet;
m:= n+1;  gramdet:= 1;
for j:=1   step  1 until n do begin
length:= 0; for i:=1 step 1 until N do
length:= length + a(i,j)**2;
for k:=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);
   A(k,j):= p;  p:= A(j,k):= q*p;
   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,m)
end i;
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;
ortho:= 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◀