DataMuseum.dk

top - metrics - download

    Length: 3072 (0xc00)
    Types: TextFile
    Names: »extgosdiag«

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

TextFile

              
external
procedure gosdiag(n,m1,m2,c,ev,x);
value n,m1,m2; integer n,m1,m2; array c,ev,x;

begin
integer i,j,j0,j1,t,t0,t1;
real h,s,k,ct,bt,p,q;
array a,bb(1:n), b(1:n+1), r(0:n+1);
q:= 0;  t0:= n*(n+1)//2;
for i:=1 step 1 until t0 do q:= q+abs(c(i));
q:= (3*'-11)*q/t0;
for i:=n step -1 until 3 do begin
  a(i):= c(t0);  t1:= t0-1;  t0:= t0-i;  s:= 0;
  for t:=t0+1 step 1 until t1 do s:= s+c(t)**2;
  ct:= c(t1);  bt:= sqrt(s);
  bt:= b(i):= if ct>0 then -bt else bt;
  if abs(bt)>0 then begin
    h:= c(t1+1):= s-ct*bt;  c(t1):= ct-bt;
    comment the elements c(t0+1) - c(t1) equals u(1) - u(i-1).
            A*u/h is calculated in r. h was stored in c(t1+1);
    j0:= t0; k:= 0;
    for t:=i-1 step -1 until 1 do begin
      j1:= j0+t; j0:= j0-t;  s:= 0;
      for j:=1 step 1 until t do s:= s + c(j0+j)*c(t0+j);
      for j:=t+1 step 1 until i-1 do begin
        s:= s + c(j1)*c(t0+j); j1:= j1+j end;
      r(t):= s:= s/h; k:= k+c(t0+t)*s
    end t;
    k:= k*0.5/h;  j0:= t0;
    comment the elements of r - k*u is calculated in r, and 
            finally c is transformed;
    for t:=i-1 step -1 until 1 do r(t):= r(t) - c(t0+t)*k;
    for t:=i-1 step -1 until 1 do begin
      j0:= j0-t;
      for j:=1 step 1 until t do
        c(j0+j):= c(j0+j) - c(t0+t)*r(j) - c(t0+j)*r(t)
    end t
  end abs(bt)>0
end i, tridiagonalization;

if n>1 then begin a(2):= c(3); b(2):= c(2) end;
a(1):= c(1);  b(1):= b(n+1):= r(0):= r(n+1):= 0;
k:= abs(a(n)) + abs(b(n));
for i:=n-1 step -1 until 1 do begin
  h:= abs(a(i)) + abs(b(i)) + abs(b(i+1));
  if h>k then k:= h
end;
for i:=1 step 1 until n do bb(i):= b(i)**2;
for i:=m1 step 1 until m2 do ev(i):= k;
h:= -k;
for j0:=m1 step 1 until m2 do begin
  k:= ev(j0);
  for ct:= (h+k)/2 while k<>ct do begin
    p:= 0; q:= 1; j:= 0;
    for i:=1 step 1 until n do begin
      s:= (a(i)-ct)*q - bb(i)*p;
      p:= q; q:= s;
      if q<0 == p>=0 then j:= j+1
    end sturm;
    if j>=j0 then begin
      k:= ct; j1:= if j<m2 then j else m2;
      for i:=j0+1 step 1 until j1 do begin
        if ct<ev(i) then ev(i):= ct end
    end else h:= ct
  end bisection;
  ev(j0):= ct;  h:= abs(a(1)-ct); t:= 1;
  for i:=2 step 1 until n do begin
    k:= abs(a(i)-ct);
    if k<h then begin h:= k; t:= i end
  end;
  r(1):= r(n):= q:= 1;  i:= 0;
  for i:=i+1 while i<t do
    r(i+1):= -((a(i)-ct)*r(i)+b(i)*r(i-1))/b(i+1);
  s:= 1/r(t);
  for i:=t-1 step -1 until 1 do begin
    r(i):= p:= r(i)*s;  q:= q+p*p end;
  i:= n+1;
  for i:=i-1 while i>t do
    r(i-1):= -((a(i)-ct)*r(i)+b(i+1)*r(i+1))/b(i);
  s:= 1/r(t);
  for i:=t+1 step 1 until n do begin
    r(i):= p:= r(i)*s;  q:= q+p*p end;
  r(t):= 1; q:= 1/sqrt(q);
  for i:=1 step 1 until n do r(i):= r(i)*q;
  t1:= n-1; t0:= 1;
  for t:=2 step 1 until t1 do begin
    t0:= t0+t;  h:= c(t0+t+1);
    if h<>0 then begin
      s:= 0;
      for i:=1 step 1 until t do s:= s + c(t0+i)*r(i);
      s:= s/h;
      for i:=1 step 1 until t do r(i):= r(i) - c(t0+i)*s
    end h<>0
  end t;
  h:= ct;  for i:=1 step 1 until n do x(j0,i):= r(i)
end j0
end;
end
▶EOF◀