DataMuseum.dk

top - metrics - download

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

└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
    └─⟦b2ec5d50f⟧ 
        └─⟦09b4e9619⟧ »thcømat« 
            └─⟦this⟧ 

TextFile

\f


message diag

diag=algol message.no
external
boolean procedure diag(n,z,d);
value n; integer n; array z,d;
begin
  integer i,j,k,l,m;
  real f,g,h,hh,b,c,p,r,s;
  array e(1:n);
  for i:= n step -1 until 2 do 
  begin
    l:= i-2; f:= z(i,i-1); g:= 0;
    for k:= 1 step 1 until l do g:= g+z(i,k)**2;
    h:= g+f*f;
    if g<=2**(-2014) then 
    begin
      e(i):= f; h:= 0 
    end else
    begin
      l:= l+1;
      g:= e(i):= if f>=0 then -sqrt(h) else sqrt(h);
      h:= h-f*g; z(i,i-1):= f-g; f:= 0;
      for j:= 1 step 1 until l do 
      begin
        z(j,i):= z(i,j)/h; g:= 0;
        for k:= 1 step 1 until j do g:= g+z(j,k)*z(i,k);
        for k:= j+1 step 1 until l do g:= g+z(k,j)*z(i,k);
        e(j):= g/h; f:= f+g*z(j,i) 
      end;
      hh:= f/(h+h);
      for j:= 1 step 1 until l do 
      begin
        f:= z(i,j); g:= e(j):= e(j)-hh*f;
        for k:= 1 step 1 until j do
        z(j,k):= z(j,k)-f*e(k)-g*z(i,k) 
      end
    end;
    d(i):= h 
  end;
  d(1):= e(1):= 0;
  for i:= 1 step 1 until n do 
  begin
    l:= i-1;
    if d(i)<>0 then for j:= 1 step 1 until l do 
    begin
      g:= 0;
      for k:= 1 step 1 until l do g:= g+z(i,k)*z(k,j);
      for k:= 1 step 1 until l do z(k,j):= z(k,j)-g*z(k,i) 
    end;
    d(i):= z(i,i); z(i,i):= 1;
    for j:= 1 step 1 until l do z(i,j):= z(j,i):= 0 
  end;

  for i:= 2 step 1 until n do e(i-1):= e(i);
  e(n):= b:= f:= 0;
  for l:= 1 step 1 until n do 
  begin
    h:= 2**(-35)*(abs(d(l))+abs(e(l)));
    if b<h then b:= h;
    for m:= l step 1 until n do
    if abs(e(m))<=b then goto cont;
cont:
    if m<>l then 
    begin
      for j:= 0,j+1 while abs(e(l))>b do
      begin
        if j=30 then 
        begin
          diag:= false; goto error 
        end;
        p:= (d(l+1)-d(l))/(2*e(l)); r:= sqrt(p**2+1);
        h:= d(l)-e(l)/(if p<0 then p-r else p+r);
        for i:= l step 1 until n do d(i):= d(i)-h;
        f:= f+h;
        p:= d(m); c:= 1; s:= 0;
        for i:= m-1 step -1 until l do 
        begin
          g:= c*e(i); h:= c*p;
          if abs(p)>=abs(e(i)) then
          begin
            c:= e(i)/p; r:= sqrt(c**2+1);
            e(i+1):= s*p*r; s:= c/r; c:= 1/r 
          end
          else 
          begin
            c:= p/e(i); r:= sqrt(c**2+1);
            e(i+1):= s*e(i)*r; s:= 1/r; c:= c/r 
          end;
          p:= c*d(i)-s*g;
          d(i+1):= h+s*(c*g+s*d(i));
          for k:= 1 step 1 until n do 
          begin
            h:= z(k,i+1);
            z(k,i+1):= s*z(k,i)+c*h;
            z(k,i):= c*z(k,i)-s*h 
          end 
        end;
        e(l):= s*p; d(l):= c*p;
      end
    end;
    d(l):= d(l)+f 
  end;
  for i:= 1 step 1 until n do 
  begin
    k:= i; p:= d(i);
    for j:= i+1 step 1 until n do if d(j)<p then 
    begin
      k:= j; p:= d(j)  
    end;
    if k<>i then 
    begin
      d(k):= d(i); d(i):= p;
      for j:= 1 step 1 until n do 
      begin
        p:= z(j,i); z(j,i):= z(j,k); z(j,k):= p 
      end 
    end 
  end;
  diag:= true;
error:
end diag
; end
▶EOF◀