procedure hqr(var a:TArray2D;n:integer;var wr,wi: TArray1D);
label 2,3,4;
var
nn,m,l,k,j,its,i,mmin:integer;
z,y,x,w,v,u,t,s,r,q,p,anorm:TElem;
begin
anorm := abs(a[1,1]);
for i := 2 to n do
for j := i-1 to n do
anorm := anorm + abs(a[i,j]);
nn := n;
t := 0;
while nn >= 1 do
begin
its := 0;
2: for l := nn downto 2 do
begin
s:= abs(a[l-1,l-1])+abs(a[l,l]);
if s = 0.0 then
s := anorm;
if abs(a[l,l-1])+s = s then
goto 3
end;
l := 1;
3: x := a[nn,nn];
if l = nn then
begin
wr[nn] := x + t;
wi[nn] := 0.0;
nn := nn-1
end
else
begin
y := a[nn-1,nn-1];
w := a[nn,nn-1]*a[nn-1,nn];
if l = nn-1 then
begin
p := 0.5*(y-x);
q := sqr(p)+w;
z := sqrt(abs(q));
x := x + t;
if q >= 0.0 then
begin
z := p + sign(z,p);
wr[nn] := x + z;
wr[nn-1] := wr[nn];
if z <> 0.0 then
wr[nn] := x - w/z;
wi[nn] := 0.0;
wi[nn] := 0.0
end
else
begin
wr[nn] := x + p;
wr[nn-1] := wr[nn];
wi[nn] := z;
wi[nn-1] := -z
end;
nn := nn-2
end
else
begin
if its = 30 then
begin
writeln('pause in routine HQR');
writeln('too many iterations') ;
readln;
end;
if (its = 10) or (its = 20) then
begin
t := t + x;
for i := 1 to nn do
a[i,i] := a[i,i] - x;
s := abs(a[nn,nn-1])+abs(a[nn-1,nn-2]);
x := 0.75*s;
y := x;
w := -0.4375*sqr(s)
end;
its := its + 1;
for m := nn-2 downto l do
begin
z := a[m,m];
r := x - z;
s := y - z;
p := (r*s-w)/a[m+1,m]+a[m,m+1];
q := a[m+1,m+1] - z - r - s;
r := a[m+2,m+1];
s := abs(p) + abs(q) + abs(r);
p := p/s;
q := q/s;
r := r/s;
if m = l then goto 4;
u := abs(a[m,m-1])*(abs(q)+abs(r));
v := abs(p) * (abs(a[m-1,m-1]) + abs(z) + abs(a[m+1,m+1]));
if u+v = v then goto 4
end;
4: for i := m+2 to nn do
begin
a[i,i-2] := 0;
if i <> m+2 then
a[i,i-3] := 0.0
end;
for k := m to nn-1 do
begin
if k <> m then
begin
p := a[k,k-1];
q := a[k+1,k-1];
if k <> nn-1 then
r := a[k+2,k-1]
else
r := 0.0;
x := abs(p) + abs(q) + abs(r);
if x <> 0.0 then
begin
p := p/x;
q := q/x;
r := r/x;
end
end;
s := sign(sqrt(sqr(p) + sqr(q) + sqr(r)),p);
if s <> 0.0 then
begin
if k = m then
begin
if l <> m then
a[k,k-1] := -a[k,k-1];
end
else
a[k,k-1] := -s*x;
p := p+s;
x := p/s;
y := q/s;
z := r/s;
q := q/p;
r := r/p;
for j := k to nn do
begin
p := a[k,j]+q*a[k+1,j];
if k <> nn-1 then
begin
p := p + r * a[k+2,j];
a[k+2,j] := a[k+2,j] - p*z
end;
a[k+1,j] := a[k+1,j] - p*y;
a[k,j] := a[k,j] -p*x
end;
mmin := min(nn,k+3);
for i := l to mmin do
begin
p := x * a[i,k]+y*a[i,k+1];
if k <> nn-1 then
begin
p := p + z *a[i,k+2];
a[i,k+2] := a[i,k+2] - p*r
end;
a[i,k+1] := a[i,k+1] - p*q;
a[i,k] := a[i,k] - p
end
end
end;
goto 2
end
end
end
end;
Recent
