code 31363;
procedure LUPZERORTPOL (N, M, B, C, ZER, EM);
value N, M; integer N, M; array B, C, ZER, EM;
begin 
procedure RATQR(N,M,POSDEF,DLAM,EPS)TRANS:(D,B2);
  value N,M,POSDEF,DLAM,EPS;
  integer N,M;
  boolean POSDEF;
  real DLAM,EPS;
  array D,B2;
comment QR ALGORITHM FOR THE COMPUTATION OF THE LOWEST EIGENVALUES
  OF A SYMMETRIC TRIDIAGONAL MATRIX. A RATIONAL VARIANT OF THE
  QR TRANSFORMATION IS USED, CONSISTING OF TWO SUCCESSIVE QD STEPS
  PER ITERATION.
  A SHIFT OF THE SPECTRUM AFTER EACH ITERATION GIVES AN ACCELERATED
  RATE OF CONVERGENCE. A NEWTON CORRECTION,DERIVED FROM THE
  CHARACTERISTIC POLYNOMIAL,IS USED AS SHIFT.
   RATQR IS IMPLEMENTED BY REINSCH AND BAUER, SEE WILKINSON AND REINSCH
   ,1971, CONTR. II-6.
  FORMATS: D,B2[1:N];
begin 
  integer I,J,K,T; real DELTA,E,EP,ERR,P,Q,QP,R,S,TOT;
comment LOWER BOUND FOR EIGENVALUES FROM GERSHGORIN, INITIAL SHIFT;
  B2[1]:= ERR:= Q:= S:= 0; TOT:= D[1];
  for I:= N step -1   until 1 do 
  begin 
    P:= Q; Q:= SQRT(B2[I]); E:= D[I]-P-Q;
    if E < TOT then TOT:= E
  end I;
  if POSDEF & TOT < 0 then TOT:= 0 else 
  for I:= 1 step 1 until N do D[I]:= D[I]-TOT;
  T:= 0;
  for K:= 1 step 1 until M do 
  begin 
NEXT QR TRANSFORMATION:  T:= T + 1;
  TOT:= TOT + S; DELTA:= D[N]-S; I:= N;
  E:= ABS(EPS*TOT); if DLAM < E then DLAM:= E;
  if DELTA < = DLAM then goto CONVERGENCE;
  E:= B2[N]/DELTA; QP:= DELTA+E; P:= 1;
  for I:= N-1 step -1 until K do 
  begin 
    Q:= D[I]-S-E; R:= Q/QP; P:= P*R+1;
    EP:= E*R; D[I+1]:= QP+EP; DELTA:= Q-EP;
    if DELTA < = DLAM then goto CONVERGENCE;
    E:= B2[I]/Q;QP:= DELTA+E; B2[I+1]:= QP*EP
  end I;
  D[K]:= QP; S:= QP/P;
  if TOT+S > TOT then goto NEXT QR TRANSFORMATION;
comment IRREGULAR END OF ITERATION,
  DEFLATE MINIMUM DIAGONAL ELEMENT;
  S:= 0; I:= K; DELTA:= QP;
  for J:= K+1 step 1 until N do 
    if D[J] < DELTA then 
    begin I:= J; DELTA:= D[J] end;
CONVERGENCE:
  if I < N then B2[I+1]:= B2[I]*E/QP;
    for J:= I-1 step -1 until K do 
      begin D[J+1]:= D[J]-S; B2[J+1]:= B2[J] end J;
    D[K]:= TOT; B2[K]:= ERR:= ERR+ABS(DELTA)
  end K;
EM[5]:=T;EM[3]:=INFNRMVEC(1,M,T,B2);
end RATQR;
    procedure DUPCEV (L, U, SHIFT, A, B);
    value L,U,SHIFT;integer L,U,SHIFT;array A,B;
    for U:=U step -1 until L do A[U]:=B[U+SHIFT];
    integer I;real NRM;
    NRM:=ABS(B[0]);
    for I:=1step 1until N-2do if C[I]+ABS(B[I])>NRM"THEN"
         NRM:=C[I]+ABS(B[I]);
     if N>1then NRM:=if NRM+1>=C[N-1]+ABS(B[N-1])then NRM+1else 
                       C[N-1]+ABS(B[N-1]);
     EM[1]:=NRM;
    DUPCEV(1,N,-1,B,B);
    DUPCEV(2,N,-1,C,C);
    RATQR (N, M, EM[6] = 1, EM[2], EM[0], B, C);
    DUPVEC (1, M, 0, ZER, B)
end LUPZERORTPOL;
        eop