?? unit2.pas
字號:
unit Unit2;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
Procedure BALANC(var A:matrx2; N:integer);
Procedure ELMHES(var A:matrx2; N:integer);
Procedure HQR(A:matrx2;N:integer;var WR,WI:array of real);
implementation
Procedure BALANC(var A:matrx2; N:integer);
Label 1,2,3;
var
I,J:integer;
RADIX,SQRDX,C,R,LAST,G,F,S:real;
begin
RADIX:=2 ;
SQRDX:=Sqr(RADIX);
1: Last:=1;
For I:=1 To N do
begin
C:=0;
R:=0;
For J:=1 To N do
begin
If J <> I Then
begin
C:=C + Abs(A[J, I]);
R:=R + Abs(A[I, J]);
end;
End;
If (C <> 0) And (R <> 0) Then
begin
G:=R / RADIX;
F:=1;
S:=C + R;
2: If C < G Then
begin
F:=F * RADIX;
C:=C * SQRDX;
GoTo 2;
end;
G:=R * RADIX;
3: If C > G Then
begin
F:=F / RADIX;
C:=C / SQRDX;
GoTo 3;
end;
If (C + R) / F < 0.95 * S Then
begin
Last:=0;
G:=1 / F;
For J:=1 To N do
A[I, J]:=A[I, J] * G;
For J:=1 To N do
A[J, I]:=A[J, I] * F;
end;
end;
end;
If Last = 0 Then GoTo 1;
End;
Procedure ELMHES(var A:matrx2; N:integer);
var
M,J,I:integer; X,Y:real;
begin
If N > 2 Then
begin
For M:=2 To N - 1 do
begin
X:=0 ;
I:=M;
For J:=M To N do
begin
If Abs(A[J, M - 1]) > Abs(X) Then
begin
X:=A[J, M - 1];
I:=J;
end;
end;
If I <> M Then
begin
For J:=M - 1 To N do
begin
Y:=A[I, J];
A[I, J]:=A[M, J];
A[M, J]:=Y;
end;
For J:=1 To N do
begin
Y:=A[J, I];
A[J, I]:=A[J, M];
A[J, M]:=Y;
end;
end;
If X <> 0 Then
begin
For I:=M + 1 To N do
begin
Y:=A[I, M - 1];
If Y <> 0 Then
begin
Y:=Y / X;
A[I, M - 1]:=Y;
For J:=M To N do
A[I, J]:=A[I, J] - Y * A[M, J];
For J:=1 To N do
A[J, M]:=A[J, M] + Y * A[J, I];
end;
end;
end;
end;
end;
end;
Procedure HQR(A:matrx2;N:integer;var WR,WI:array of real);
Label 1,2,3,4;
var
I,J,NN,III,K,ITS,L,M:integer;
ANORM,T,S,X,Y,W,R,Q,P,U,ZZ,Z,AAA,BBB,V:real;
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;
1: If NN >= 1 Then
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 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;
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:=P*P + W;
Z:=Sqrt(Abs(Q));
X:=X + T;
If Q >= 0 Then
begin
if P>=0 then
ZZ:=1
else
ZZ:=-1;
Z:=P + Abs(Z) * ZZ;
WR[NN]:=Z + X;
WR[NN - 1]:=WR[NN];
If Z <> 0 Then WR[NN]:=X - W / Z;
WI[NN]:=0;
WI[NN - 1]:=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 ShowMessage(' too many iterations ');
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 * S * 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));
BBB:=Abs(A[M + 1, M + 1]);
AAA:=Abs(A[M - 1, M - 1]) + Abs(Z) + BBB;
V:=Abs(P) * AAA;
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;
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];
R:=0;
If K <> NN - 1 Then R:=A[K + 2, K - 1];
X:=Abs(P) + Abs(Q) + Abs(R);
If X <> 0 Then
begin
P:=P / X;
Q:=Q / X;
R:=R / X;
end;
end;
if P>=0 then
ZZ:=1
else
ZZ:=-1;
S:=Sqrt(P*P + Q*Q + R*R) * ZZ;
If S <> 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;
If NN > K + 3 Then
III:=K + 3
Else
III:=NN;
For I:=L To III 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;
GoTo 1;
end;
end;
end.
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -