DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »tminn2d«

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

TextFile

\f


message minn2d

minn2d=algol message.no
external
real procedure minn2d(n,f,X,eps,minimum);
value n; integer n; real procedure f;
array X; real eps; boolean minimum;
begin
  array A(1:n*(n+1)/2),B,Y,PHI(1:n);
  real q,q0,x,y,z,b,e,del;
  integer i,i1,j,k,ti,tj;
  eps:=eps**2; del:=1;
  q0:=f(X,B,A);
NEW:
  minimum:=true;
  comment Cholesky factorization;
  ti:=0;
  for i:=1 step 1 until n do 
  begin
    i1:=i-1; tj:=0;
    for j:=1 step 1 until i1 do 
    begin
      x:=A(ti+j);
      for k:=j-1 step -1 until 1 do
      x:=-A(ti+k)*A(tj+k)+x;
      A(ti+j):=x;
      tj:=tj+j 
    end;
    x:=A(ti+i);
    for k:=i1 step -1 until 1 do 
    begin
      y:=A(ti+k); z:=A(ti+k):=A(k*(k+1)/2)*y;
      x:=-y*z+x 
    end;
    if x>0 then 
    begin
      ti:=ti+i; A(ti):=1/x 
    end
    else 
    begin
      minimum:=false; i:=n 
    end
  end;
  comment check norm of gradient;
  b:=0;
  for i:=1 step 1 until n do b:=B(i)**2+b;
  if b<=eps and minimum then goto END;
  comment compute phi-vector;
  if minimum then 
  begin
    j:=0;
    for i:=1 step 1 until n do 
    begin
      y:=B(i);
      for k:=i-1 step -1 until 1 do y:=y-A(j+k)*PHI(k);
      PHI(i):=y;
      j:=j+i 
    end;
    z:=0;
    for i:=n step -1 until 1 do 
    begin
      j:=i*(i+1)/2;
      y:=PHI(i)*A(j);
      for k:=i+1 step 1 until n do 
      begin
        y:=-A(j+i)*PHI(k)+y;
        j:=j+k 
      end;
      PHI(i):=y; z:=B(i)*y+z 
    end
  end
  else 
  begin
    for i:=1 step 1 until n do PHI(i):=B(i);
    z:=b 
  end;
  comment find minimum in phi direction;
  del:=if del>=.5 then 1 else 2*del;
ITER:
  e:=0;
  for i:=1 step 1 until n do 
  begin
    x:=X(i); y:=Y(i):=x-PHI(i)*del;
    e:=abs(x-y)+e 
  end;
  comment if X<>Y then continue search;
  if e>0 then 
  begin
    q:=f(Y,B,A);
    x:=(q-q0)/del/z;
    if x>-.25 then 
    begin
      if x<0 then 
      begin
        y:=0;
        for i:=1 step 1 until n do y:=PHI(i)*B(i)+y
      end
      else y:=1;
      if y>=0 then 
      begin
        del:=.5*del/(x+1);
        goto ITER 
      end 
    end;
    q0:=q;
    for i:=1 step 1 until n do X(i):=Y(i);
    goto NEW 
  end;
END:
  eps:=sqrt(b); minn2d:=q0
end minn2d
; end
▶EOF◀