⟦cb112225b⟧

TextFile

\f



external
procedure banddiag(Nt,n,m1,m2,ev,H,U);
value Nt,n,m1,m2; integer Nt,n,m1,m2; array ev,H,U;
begin
integer i,j,k,jk,m,r,N1,N2,maxr,maxk,ugk;
real g,p,c,s,c2,s2,cs,u0,u1,lambda,norm;
array b(1:Nt,0:n), ta,tb(1:Nt);
N1:= Nt-1; N2:=Nt-2;

comment Tridiagonalization by method of
        R.H.Schwartz, Numer. Math. 12, 231 (1968).
        n is the number of bidiagonals;

for i:=1 step 1 until Nt do begin k:= if n<Nt-i then n else Nt-i;
for j:=0 step 1 until k do b(i,j):= H(i+j,j) end;
for j:=1 step 1 until N2 do begin
  maxr:= if Nt-j<n then Nt-j else n;
  for r:=maxr step -1 until 2 do begin
    for i:=j+r step n until Nt do begin
      if i=j+r then begin
        if b(j,r)=0 then goto endr;
        p:=-b(j,r-1)/b(j,r); ugk:= j;
      end else begin
        if g=0 then goto endr;
        p:=-b(i-n-1,n)/g; ugk:= i-n;
      end;
      s2:= 1/(1+p*p); s:= sqrt(s2); c:= p*s; c2:= c*c; cs:= c*s;
      u0:= b(i-1,0)*c2-b(i-1,1)*2*cs+b(i,0)*s2;
      u1:= b(i-1,0)*s2+b(i-1,1)*2*cs+b(i,0)*c2;
      b(i-1,1):= (b(i-1,0)-b(i,0))*cs+b(i-1,1)*(c2-s2);
      b(i-1,0):= u0; b(i,0):= u1;
      for k:=ugk step 1 until i-2 do begin
        u0      := b(k,i-k-1)*c-b(k,i-k)*s;
        b(k,i-k):= b(k,i-k-1)*s+b(k,i-k)*c;
        b(k,i-k-1):= u0;
      end k;
      if i<>j+r then b(i-n-1,n):= b(i-n-1,n)*c-g*s;
      maxk:= if Nt-i<n-1 then Nt-i else n-1;
      for k:=1 step 1 until maxk do begin
        u0    := b(i-1,k+1)*c-b(i,k)*s;
        b(i,k):= b(i-1,k+1)*s+b(i,k)*c;
        b(i-1,k+1):= u0;
      end k;
      if i+n<=Nt then begin
        g:= -b(i,n)*s; b(i,n):= b(i,n)*c;
      end;
    end i;
endr: end r;
end j;
b(Nt,1):= 0;
for i:=Nt step -1 until 1 do begin
   p:= abs b(i,0) + abs b(i,1);
   ta(i):= b(i,0); if i>1 then tb(i):= b(i-1,1)**2;
   if p>g then g:= p end;
c:= -g; for i:=m1 step 1 until m2 do ev(i):= g;
for m:=m1 step 1 until m2 do begin
g:= ev(m);
for lambda:= (g+c)/2 while g<>lambda do begin
  s:= 0; u1:= 0; u0:= 1;
  for i:=1 step 1 until Nt do begin
    p:= (ta(i)-lambda)*u0 - tb(i)*u1;
    u1:= u0; u0:= p;
    if u1>=0 == u0<0 then s:= s+1;
    if s=m2 then i:=Nt
  end i;
  if s>=m then begin
    g:= lambda;
    for i:=m+1 step 1 until s do begin
      if lambda<ev(i) then ev(i):= lambda end
  end else c:= lambda
end lambda;
ev(m):= lambda;
begin
array a(1:(n+1)*n//2), V(-n:Nt+n);
b(1,0):= H(1,0)-lambda; p:= abs b(1,0); s:= 1;
for i:=2 step 1 until Nt do begin
   b(i,0):= H(i,0)-lambda; u0:= abs b(i,0);
   if p>u0 then begin p:= u0; s:= i end
end;
for i:=1 step 1 until Nt do
for j:=1 step 1 until n do b(i,j):= H(i,j);
for i:=Nt+n step -1 until -n do V(i):= 0;
for j:=n step -1 until 1 do begin
   if s-j >0 then V(s-j):= -b(s,j);
   if s+j<=Nt then V(s+j-1):= -b(s+j,j) end;
for i:=s+1 step 1 until Nt do
for j:=0 step 1 until n do b(i-1,j):= b(i,j);
for i:= if s+n>Nt then N1-s else n-1
   step -1 until 0 do begin
   for j:=i+2 step 1 until n do b(s+i,j-1):= b(s+i,j);
   b(s+i,n):= 0 end;
for i:=1 step 1 until N1 do begin
   jk:= if i>n then n else i-1;
   for k:=1 step 1 until jk do a(k):= b(i,k)/b(i-k,0);
   for j:= if i+n>=Nt then N1-i else n step -1 until 0 do
      b(i+j,j):= b(i+j,j)-sum(a(k)*b(i+j,j+k),k,1,n-j);
   V(i):= (V(i)-sum(b(i,i-k)*V(k),k,i-jk,i-1))/b(i,0);
   for k:=1 step 1 until jk do b(i,k):= a(k)
end i;
norm:= 1;
for i:=N1 step -1 until 1 do begin
   p:= V(i):= V(i)-sum(b(i+k,k)*V(i+k),
                       k,1,if i+n>=Nt then N1-i else n);
   norm:= norm+p*p
end;
norm:= U(m,s):= 1/sqrt(norm);
for i:= N1 step -1 until s do U(m,i+1):= V(i)*norm;
for i:=s-1 step -1 until 1 do U(m,i):= V(i)*norm
end b,V end m end;
end
▶EOF◀
DataMuseum.dk

RC4000/8000/9000

⟦cb112225b⟧ TextFile

Derivation

TextFile
