⟦0a000276f⟧

TextFile

cptridql=algol list.yes index.no

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

comment The procedure diagonalizes a hermitian matrix of 
        order n. The upper triangular part of this is stored
        with the real part on and below the diagonal and the
        imaginary part above the diagonal in the quadratic
        array c(1:n,1:n). The eigenvalues in ascending order
        are formed in a(1:n), and the unitary transformation
        matrix ( a = x c x' ) is formed with the real part
        in c and the imaginary part in the array x(1:n,1:n).
        First part of the procedure is a modified Householder
        tridiagonalization;

begin
integer i,j,k,l,m;
real f,g,h,br,bi,qr,qi,q,s,t,macheps;
array b,pr,pi(1:n);
macheps:= 3'-11; q:= 0;
for i:=1 step 1 until n do
for j:=1 step 1 until n do q:= q+abs(c(i,j));
q:= macheps*q/(n*n);
for i:=n step -1 until 3 do begin
  a(i):= c(i,i); l:= i-1; qr:= c(i,l); qi:= c(l,i);
  s:= f:= qr*qr+qi*qi;
  for j:=i-2 step -1 until 1 do s:= s+c(i,j)**2+c(j,i)**2;
  b(l):= s:= sqrt(s);
  if s>q then begin
    f:= sqrt(f);  c(i,i):= h:= (s+f)*s;
    if f>0 then begin
      pr(i):= br:= -qr/f; pi(i):= bi:= qi/f;
      c(i,l):= qr-s*br;  c(l,i):= qi+s*bi
    end else begin
      c(i,l):= s; pi(i):= c(l,i):= bi:= 0; pr(i):= br:= -1 end;
    comment the elements c(i,1) - c(i,i-1) and c(1,i) - c(i-1,i)
            contain real and imaginary parts of u(1) - u(i-1).
            h was stored in c(i,i). p = c*u/h is calculated
            in pr and pi. g = u'*p/2h is formed as well;
    g:= 0;
    for j:=1 step 1 until l do begin
      qr:= c(j,j)*c(i,j); qi:= c(j,j)*c(j,i);
      for k:=j-1 step -1 until 1 do begin
        qr:= qr+c(j,k)*c(i,k)+c(k,j)*c(k,i);
        qi:= qi-c(k,j)*c(i,k)+c(j,k)*c(k,i)
      end;
      for k:= j+1 step 1 until l do begin
        qr:= qr+c(k,j)*c(i,k)-c(j,k)*c(k,i);
        qi:= qi+c(j,k)*c(i,k)+c(k,j)*c(k,i)
      end;
      pr(j):= qr:= qr/h;  pi(j):= qi:= qi/h;
      g:= g+qr*c(i,j)+qi*c(j,i);
    end j;
    g:= g/(h+h);
    comment the elements of p - g*u is calculated in p, and 
            finally c is transformed;
    for j:=1 step 1 until l do begin
      pr(j):= pr(j)-c(i,j)*g;
      pi(j):= pi(j)-c(j,i)*g
    end;
    for j:=1 step 1 until l do
    for k:=1 step 1 until j do begin
      c(j,k):= c(j,k)-pr(j)*c(i,k)-pi(j)*c(k,i)
                     -pr(k)*c(i,j)-pi(k)*c(j,i);
      if j>k then
      c(k,j):= c(k,j)-pr(j)*c(k,i)+pi(j)*c(i,k)
                     +pr(k)*c(j,i)-pi(k)*c(i,j)
    end;
    comment exp(i*delta(i)) used to make b(i-1) positive
            was evaluated and stored in p(i). The elements of
            c with index (i-1) are transformed accordingly;
    for j:=l-1 step -1 until 1 do begin
      qr:= c(l,j); qi:= c(j,l);
      c(l,j):= br*qr+bi*qi; c(j,l):= br*qi-bi*qr
end end sigma>q
end i, tridiagonalization;

a(1):= c(1,1); c(1,1):= 1; x(1,1):= b(n):= 0;
if n>1 then begin
  a(2):= c(2,2); br:= c(2,1); bi:= c(1,2);
  b(1):= f:= sqrt(br*br+bi*bi);
  if f>0 then begin c(1,1):= br/f; x(1,1):= -bi/f
end end;
m:= n+1;
for i:=3 step 1 until m do begin
  l:= i-1;
  if b(l)<=q then begin
    c(l,l):= 1; x(l,l):= 0;
    for j:=l-1 step -1 until 1 do
      c(j,l):= c(l,j):= x(j,l):= x(l,j):= 0
  end else begin
  h:= c(i,i);
  comment u is transferred to p;
  for j:=1 step 1 until l do begin
    pr(j):= c(i,j); pi(j):= c(j,i)
  end;
  g:= 1/b(l); g:= h*g*g; qr:= pr(i); qi:= pi(i);
  c(l,l):= qr-qr*g;
  x(l,l):= qi-qi*g;
  br:= (qr*pr(l)-qi*pi(l))/h;
  bi:= (qi*pr(l)+qr*pi(l))/h;
  for j:=i-2 step -1 until 1 do begin
    qr:= qi:= 0;
    for k:=i-2 step -1 until 1 do begin
      qr:= qr+c(j,k)*pr(k)-x(j,k)*pi(k);
      qi:= qi+c(j,k)*pi(k)+x(j,k)*pr(k)
    end;
    c(j,l):= x(j,l):= 0; qr:= qr/h; qi:= qi/h;
    for k:=1 step 1 until l do begin
      c(j,k):= c(j,k)-qr*pr(k)-qi*pi(k);
      x(j,k):= x(j,k)-qi*pr(k)+qr*pi(k)
    end;
    c(l,j):= -br*pr(j)-bi*pi(j);
    x(l,j):= -bi*pr(j)+br*pi(j)
end end j end i;

comment Now we find the eigenvalues and eigenvectors of the
        tridiagonal 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 ascending order. The eigenvectors are
        formed in the array x(1:n,1:n), overwritting the
        accumulated transformations supplied by the tridia-
        gonalization above (Numer. Math. 11, 293-306 (1968));
br:= f:= 0;
for l:=1 step 1 until n do begin
  j:= 0; h:= macheps*(abs a(l) + abs b(l));
  if br<h then br:= h;
  comment look for small sub-diagonal elements;
  for m:= l step 1 until n do
  if abs b(m)<=br then goto cont;
cont:
  if m=l then goto root;
nextit:
  if j=30 then goto fail;
  j:= j+1;
  comment form shift;
  t:= (a(l+1)-a(l))/(2*b(l)); q:= sqrt(t*t+1);
  h:= a(l)-b(l)/(if t<0 then t-q else t+q);
  for i:=l step 1 until n do a(i):= a(i)-h;
  f:= f+h;
  comment QL transformation;
  t:= a(m); bi:= 1; s:= 0;
  for i:= m-1 step -1 until l do begin
    g:= bi*b(i); h:= bi*t;
    if abs t>=abs b(i) then begin
      bi:= b(i)/t; q:= sqrt(bi*bi+1);
      b(i+1):= s*t*q; s:= bi/q; bi:= 1/q
    end else begin
      bi:= t/b(i); q:= sqrt(bi*bi+1);
      b(i+1):= s*b(i)*q; s:= 1/q; bi:= bi/q
    end;
    t:= bi*a(i)-s*g;
    a(i+1):= h+s*(bi*g+s*a(i));
    comment form vector;
    for k:=1 step 1 until n do begin
      qr:= x(i,k); qi:= x(i+1,k);
      x(i+1,k):= s*qr+bi*qi;
      x(i,k)  := bi*qr-s*qi;
      qr:= c(i,k); qi:= c(i+1,k);
      c(i+1,k):= s*qr+bi*qi;
      c(i,k)  := bi*qr-s*qi
    end k
  end i;
  b(l):= s*t; a(l):= bi*t;
  if abs b(l)>br then goto nextit;
root:
  a(l):= a(l)+f
end l;
comment order eigenvalues and eigenvectors;
for i:=1 step 1 until n do begin
  k:= i; t:= a(i);
  for j:=i+1 step 1 until n do if a(j)<t then
  begin k:= j; t:= a(j) end;
  if k<>i then begin
    a(k):= a(i); a(i):= t;
    for j:=1 step 1 until n do begin
      t:= x(i,j); x(i,j):= x(k,j); x(k,j):= t;
      t:= c(i,j); c(i,j):= c(k,j); c(k,j):= t
    end j
end end;
fail:
end;
end
▶EOF◀
DataMuseum.dk

RC4000/8000/9000

⟦0a000276f⟧ TextFile

Derivation

TextFile
