1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
{ ******************************************************************
Reduction of a square matrix to upper Hessenberg form
****************************************************************** }
unit uelmhes;
interface
uses
utypes;
procedure ElmHes(A : PMatrix;
Lb, Ub, I_low, I_igh : Integer;
I_int : PIntVector);
implementation
procedure ElmHes(A : PMatrix;
Lb, Ub, I_low, I_igh : Integer;
I_int : PIntVector);
{ ------------------------------------------------------------------
This procedure is a translation of the EISPACK subroutine Elmhes
Given a real general matrix, this procedure reduces a submatrix
situated in rows and columns I_low through I_igh to upper
Hessenberg form by stabilized elementary similarity transformations.
On input:
A contains the input matrix.
Lb, Ub are the lowest and highest indices
of the elements of A.
I_low and I_igh are integers determined by the balancing procedure
Balance. If Balance has not been used, set I_low = Lb, I_igh = Ub.
On output:
A contains the Hessenberg matrix. The multipliers which were used
in the reduction are stored in the remaining triangle under the
Hessenberg matrix.
I_int contains information on the rows and columns interchanged
in the reduction. Only elements I_low through I_igh are used.
------------------------------------------------------------------ }
var
I, J, M, La, Kp1, Mm1, Mp1 : Integer;
X, Y : Float;
begin
La := I_igh - 1;
Kp1 := I_low + 1;
if La < Kp1 then Exit;
for M := Kp1 to La do
begin
Mm1 := M - 1;
X := 0.0;
I := M;
for J := M to I_igh do
if Abs(A^[J]^[Mm1]) > Abs(X) then
begin
X := A^[J]^[Mm1];
I := J;
end;
I_int^[M] := I;
{ Interchange rows and columns of A }
if I <> M then
begin
for J := Mm1 to Ub do
begin
Y := A^[I]^[J];
A^[I]^[J] := A^[M]^[J];
A^[M]^[J] := Y;
end;
for J := Lb to I_igh do
begin
Y := A^[J]^[I];
A^[J]^[I] := A^[J]^[M];
A^[J]^[M] := Y;
end;
end;
if X <> 0.0 then
begin
Mp1 := M + 1;
for I := Mp1 to I_igh do
begin
Y := A^[I]^[Mm1];
if Y <> 0.0 then
begin
Y := Y / X;
A^[I]^[Mm1] := Y;
for J := M to Ub do
A^[I]^[J] := A^[I]^[J] - Y * A^[M]^[J];
for J := Lb to I_igh do
A^[J]^[M] := A^[J]^[M] + Y * A^[J]^[I];
end;
end;
end;
end;
end;
end.
|