DataMuseum.dk

top - metrics - download

    Length: 5376 (0x1500)
    Types: TextFile
    Names: »exttridql«

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

TextFile

tridql=algol list.yes index.no

external
procedure tridql3(n,a,x);
value n; integer n; array a,x;

comment The procedure diagonalizes a symmetric matrix of order n.
        The lower triangular part of this is stored in the array
        x(1:n,1:n). The eigenvalues in ascending order are formed
        in a(1:n), and the orthonormal transformation matrix u
        ( a = u c utr ) is overwritting c in the array x. First
        part of the procedure is a Householder tridiagonali-
        zation (Num. Math. 11,181-195(1968));

begin
integer i,j,k,l,m;
boolean sorter;
real f,g,h,p,q,c,s,t,tol;
integer array pm(1:n);
array b(1:n);
for i:=1 step 1 until n do begin
  pm(i):= i; b(i):= x(i,i);
  for j:=i+1 step 1 until n do x(i,j):= x(j,i)
end;
k:= 1; m:= n-1; sorter:= false;
for i:=2 step 1 until m do begin
  g:= 0; l:= i;
  for j:=i step 1 until n do
  if g<abs x(k,pm(j)) then begin
    g:= abs x(k,pm(j)); l:= j
  end;
  k:= pm(l);
  for j:=l-1 step -1 until i do pm(j+1):= pm(j);
  pm(i):= k;  sorter:= sorter or k<>i
end;
if sorter then begin
  for i:=1 step 1 until n do begin
    k:= pm(i); x(i,i):= b(k);
    for j:=i+1 step 1 until n do begin
      l:= pm(j);
      x(j,i):= if k<l then x(k,l) else x(l,k)
end end end;
tol:= 0.0 shift 35 add (-2012 extract 12);
for i:=n step -1 until 2 do begin
  l:= i-1; f:= x(i,l); g:= f*f; s:= 0;
  for j:=i-2 step -1 until 1 do s:= x(i,j)*x(i,j)+s;
  if s>tol then begin
    s:= s+g; g:= sqrt(s);
    g:= b(l):= if f>0 then -g else g;
    a(i):= h:= -g*f+s; x(i,l):= f-g; f:= 0;
    for j:=1 step 1 until l do begin
      s:= 0;
      for k:=1   step 1 until j do s:= x(j,k)*x(i,k)+s;
      for k:=j+1 step 1 until l do s:= x(k,j)*x(i,k)+s;
      a(j):= s:= s/h; f:= x(i,j)*s+f;
    end j;
    s:= f/(h+h);
    for j:=1 step 1 until l do begin
      f:= x(i,j); a(j):= g:= a(j)-s*f;
      for k:=1 step 1 until j do x(j,k):= -a(k)*f-x(i,k)*g+x(j,k)
    end j;
  end s>tol else begin
    b(l):= f; a(i):= 0
end end i;
a(1):= b(n):= 0;

comment accumulation of transformation matrices;
for i:=1 step 1 until n do begin
  l:= i-1; h:= a(i);
  if h<>0 then begin
    for j:=1 step 1 until l do begin
      s:= 0;
      for k:=1 step 1 until l do s:= x(i,k)*x(j,k)+s;
      s:= s/h;
      for k:=1 step 1 until l do x(j,k):= -x(i,k)*s+x(j,k)
  end j end;
  a(i):= x(i,i); x(i,i):= 1;
  for j:=1 step 1 until l do x(i,j):= x(j,i):= 0
end i;
comment Now we find the eigenvalues and eigenvectors of the tri-
        diagonal matrix given with its diagonal elements in the
        array a(1:n) and its subdiagonal elements in the first
        n-1 stores of the array b(1:n), using QL transformations.
        The eigenvalues are overwritten in the array a in ascen-
        ding order. The eigenvectors are formed in x(1:n,1:n),
        overwritting the accumulated transformations supplied
        by the tridiagonalization (Num. Math. 15, 450 (1970));
tol:= 0.0 shift 35 add (-34 extract 12);
for l:=1 step 1 until n do begin
  j:= 0;
  comment look for small sub-diagonal elements;
nextit:
  m:= l-1;
  for m:=m+1 while m<n and abs b(m)>(abs a(m)+abs a(m+1))*tol do;
  if m=l then goto root;
  if j=30 then goto fail;
  j:= j+1; q:= a(l);
  comment form shift;
  g:= (a(l+1)-q)/(2*b(l)); f:= sqrt(g*g+1);
  q:= a(m)-q+b(l)/(if g<0 then g-f else g+f);
  comment QL transformation;
  s:= 1; f:= c:= t:= 0;
  for i:= m-1 step -1 until l do begin
    g:= b(i);
    if c=0 then begin p:= g*s; h:= -t*p+g end
           else begin h:= g*c; p:= -t*h+g end;
    if abs p<=abs q then begin
      c:= p/q; p:= sqrt(c*c+1);
      b(i+1):= p*q; s:= c/p; t:= c/(1+p); c:= 0;
      g:= a(i+1)-f; q:= a(i)-g; p:= h+h;
      f:= ((-t*p+q)*s+p)*s;  a(i+1):= g+f;
      q:= ((-q*t-p)*s+q)*s+h;
      for k:=1 step 1 until n do begin
        g:= x(i,k); p:= x(i+1,k);
        x(i  ,k):= -(g*t+p)*s+g;
        x(i+1,k):= -(p*t-g)*s+p
    end end else begin
      c:= q/p; q:= sqrt(c*c+1);
      b(i+1):= p*q; t:= c/(1+q); c:= c/q;
      g:= a(i+1)-f; q:= a(i)-g; p:= h+h;
      f:= ((-t*p-q)*c+p)*c; a(i+1):= a(i)+f; f:= f+q;
      q:= ((-t*q+p)*c+q)*c-h;
      for k:=1 step 1 until n do begin
        g:= x(i,k); p:= x(i+1,k);
        x(i  ,k):=  (p*t+g)*c-p;
        x(i+1,k):= -(g*t-p)*c+g
      end k
  end end i;
  b(l):= q; a(l):= a(l)-f; b(m):= 0;
  goto nextit;
root:
end l;
comment order eigenvalues and eigenvectors;
for i:=1 step 1 until n do begin
  k:= i; p:= a(i);
  for j:=i+1 step 1 until n do if a(j)<p then
  begin k:= j; p:= a(j) end;
  if k<>i then begin
    a(k):= a(i); a(i):= p;
    for j:=1 step 1 until n do begin
      p:= x(i,j); x(i,j):= x(k,j); x(k,j):= p
    end j
end end;
if sorter then begin
  for i:=1 step 1 until n do begin
    for j:=2 step 1 until n do b(pm(j)):= x(i,j);
    for j:=2 step 1 until n do x(i,j):= b(j)
end end;
fail:
end;
end
▶EOF◀