1
{ **********************************************************************
4
* (c) J. Debord, May 2001 *
5
**********************************************************************
6
Procedures for computing eigenvalues and eigenvectors
7
**********************************************************************
9
1) Borland's Numerical Methods Toolbox : Jacobi
10
2) 'Numerical Recipes' by Press et al. : EigenVals, RootPol
11
********************************************************************** }
20
function Jacobi(A : PMatrix; Lbound, Ubound, MaxIter : Integer;
21
Tol : Float; V : PMatrix; Lambda : PVector) : Integer;
22
{ ----------------------------------------------------------------------
23
Eigenvalues and eigenvectors of a symmetric matrix by the iterative
25
----------------------------------------------------------------------
26
Input parameters : A = matrix
27
Lbound = index of first matrix element
28
Ubound = index of last matrix element
29
MaxIter = maximum number of iterations
30
Tol = required precision
31
----------------------------------------------------------------------
32
Output parameters : V = matrix of eigenvectors (stored by lines)
33
Lambda = eigenvalues in decreasing order
34
----------------------------------------------------------------------
35
Possible results : MAT_OK
37
----------------------------------------------------------------------
38
NB : 1. The eigenvectors are normalized, with their first component > 0
39
2. This procedure destroys the original matrix A
40
---------------------------------------------------------------------- }
42
function EigenVals(A : PMatrix; Lbound, Ubound : Integer;
43
Lambda_Re, Lambda_Im : PVector) : Integer;
44
{ ----------------------------------------------------------------------
45
Eigenvalues of a general square matrix
46
----------------------------------------------------------------------
47
Input parameters : A = matrix
48
Lbound = index of first matrix element
49
Ubound = index of last matrix element
50
----------------------------------------------------------------------
51
Output parameters : Lambda_Re = real part of eigenvalues
52
Lambda_Im = imaginary part of eigenvalues
53
----------------------------------------------------------------------
54
Possible results : MAT_OK
56
----------------------------------------------------------------------
57
NB : This procedure destroys the original matrix A
58
---------------------------------------------------------------------- }
60
function EigenVect(A : PMatrix; Lbound, Ubound : Integer;
61
Lambda, Tol : Float; V : PVector) : Integer;
62
{ ----------------------------------------------------------------------
63
Computes the eigenvector associated to a real eigenvalue
64
----------------------------------------------------------------------
65
Input parameters : A = matrix
66
Lbound = index of first matrix element
67
Ubound = index of last matrix element
69
Tol = required precision
70
----------------------------------------------------------------------
71
Output parameters : V = eigenvector
72
----------------------------------------------------------------------
73
Possible results : MAT_OK
75
----------------------------------------------------------------------
76
NB : 1. The eigenvector is normalized, with its first component > 0
77
2. The function returns only one eigenvector, even if the
78
eigenvalue has a multiplicity greater than 1.
79
---------------------------------------------------------------------- }
81
procedure DivLargest(V : PVector; Lbound, Ubound : Integer;
83
{ ----------------------------------------------------------------------
84
Normalizes an eigenvector V by dividing by the element with the
85
largest absolute value
86
---------------------------------------------------------------------- }
88
function RootPol(Coef : PVector; Deg : Integer;
89
X_Re, X_Im : PVector) : Integer;
90
{ ----------------------------------------------------------------------
91
Real and complex roots of a real polynomial by the method of the
93
----------------------------------------------------------------------
94
Input parameters : Coef = coefficients of polynomial
95
Deg = degree of polynomial
96
----------------------------------------------------------------------
97
Output parameters : X_Re = real parts of root (in increasing order)
98
X_Im = imaginary parts of root
99
----------------------------------------------------------------------
100
Possible results : MAT_OK
102
---------------------------------------------------------------------- }
106
function Jacobi(A : PMatrix; Lbound, Ubound, MaxIter : Integer;
107
Tol : Float; V : PMatrix; Lambda : PVector) : Integer;
109
SinTheta, CosTheta, TanTheta, Tan2Theta : Float;
110
CosSqr, SinSqr, SinCos, SumSqrDiag : Float;
111
AII, AJJ, AIJ, AIK, AJK, VIK, VJK, D : Float;
112
I, J, K, Iter : Integer;
116
for I := Lbound to Ubound do
117
for J := Lbound to Ubound do
126
for I := Lbound to Ubound do
127
SumSqrDiag := SumSqrDiag + Sqr(A^[I]^[I]);
130
for I := Lbound to Pred(Ubound) do
131
for J := Succ(I) to Ubound do
132
if Abs(A^[I]^[J]) > Tol * SumSqrDiag then
136
{ Calculate rotation }
137
D := A^[I]^[I] - A^[J]^[J];
138
if Abs(D) > MACHEP then
140
Tan2Theta := D / (2.0 * A^[I]^[J]);
141
TanTheta := - Tan2Theta + Sgn(Tan2Theta) *
142
Sqrt(1.0 + Sqr(Tan2Theta));
143
CosTheta := 1.0 / Sqrt(1.0 + Sqr(TanTheta));
144
SinTheta := CosTheta * TanTheta;
148
CosTheta := SQRT2DIV2; { Sqrt(2)/2 }
149
SinTheta := Sgn(A^[I]^[J]) * SQRT2DIV2;
153
CosSqr := Sqr(CosTheta);
154
SinSqr := Sqr(SinTheta);
155
SinCos := SinTheta * CosTheta;
156
AII := A^[I]^[I] * CosSqr + 2.0 * A^[I]^[J] * SinCos
157
+ A^[J]^[J] * SinSqr;
158
AJJ := A^[I]^[I] * SinSqr - 2.0 * A^[I]^[J] * SinCos
159
+ A^[J]^[J] * CosSqr;
160
AIJ := (A^[J]^[J] - A^[I]^[I]) * SinCos
161
+ A^[I]^[J] * (CosSqr - SinSqr);
162
for K := Lbound to Ubound do
163
if not(K in [I, J]) then
165
AIK := A^[I]^[K] * CosTheta + A^[J]^[K] * SinTheta;
166
AJK := - A^[I]^[K] * SinTheta + A^[J]^[K] * CosTheta;
177
{ Rotate eigenvectors }
178
for K := Lbound to Ubound do
180
VIK := CosTheta * V^[I]^[K] + SinTheta * V^[J]^[K];
181
VJK := - SinTheta * V^[I]^[K] + CosTheta * V^[J]^[K];
186
until Done or (Iter > MaxIter);
188
{ The diagonal terms of the transformed matrix are the eigenvalues }
189
for I := Lbound to Ubound do
190
Lambda^[I] := A^[I]^[I];
192
if Iter > MaxIter then
194
Jacobi := MAT_NON_CONV;
198
{ Sort eigenvalues and eigenvectors }
199
for I := Lbound to Pred(Ubound) do
203
for J := Succ(I) to Ubound do
204
if Lambda^[J] > D then
209
FSwap(Lambda^[I], Lambda^[K]);
210
SwapRows(I, K, V, Lbound, Ubound);
213
{ Make sure that the first component of each eigenvector is > 0 }
214
for I := Lbound to Ubound do
215
if V^[I]^[Lbound] < 0.0 then
216
for J := Lbound to Ubound do
217
V^[I]^[J] := - V^[I]^[J];
222
procedure Balance(A : PMatrix; Lbound, Ubound : Integer);
223
{ Balances the matrix, i.e. reduces norm without affecting eigenvalues }
225
RADIX = 2; { Base used for machine computations }
227
I, J, Last : Integer;
228
C, F, G, R, S, Sqrdx : Float;
233
for I := Lbound to Ubound do
237
for J := Lbound to Ubound do
240
C := C + Abs(A^[J]^[I]);
241
R := R + Abs(A^[I]^[J]);
243
if (C <> 0.0) and (R <> 0.0) then
259
if (C + R) / F < 0.95 * S then
263
for J := Lbound to Ubound do
264
A^[I]^[J] := A^[I]^[J] * G;
265
for J := Lbound to Ubound do
266
A^[J]^[I] := A^[J]^[I] * F;
273
procedure ElmHes(A : PMatrix; Lbound, Ubound : Integer);
274
{ Reduces the matrix to upper Hessenberg form by elimination }
279
for M := Succ(Lbound) to Pred(Ubound) do
283
for J := M to Ubound do
284
if Abs(A^[J]^[M - 1]) > Abs(X) then
291
for J := Pred(M) to Ubound do
292
FSwap(A^[I]^[J], A^[M]^[J]);
293
for J := Lbound to Ubound do
294
FSwap(A^[J]^[I], A^[J]^[M]);
297
for I := Succ(M) to Ubound do
304
for J := M to Ubound do
305
A^[I]^[J] := A^[I]^[J] - Y * A^[M]^[J];
306
for J := Lbound to Ubound do
307
A^[J]^[M] := A^[J]^[M] + Y * A^[J]^[I];
311
for I := (Lbound + 2) to Ubound do
312
for J := Lbound to (I - 2) do
316
function Hqr(A : PMatrix; Lbound, Ubound : Integer;
317
Lambda_Re, Lambda_Im : PVector) : Integer;
318
{ Finds the eigenvalues of an upper Hessenberg matrix }
321
I, Its, J, K, L, M, N : Integer;
322
Anorm, P, Q, R, S, T, U, V, W, X, Y, Z : Float;
324
function Sign(A, B : Float) : Float;
326
if B < 0.0 then Sign := - Abs(A) else Sign := Abs(A)
330
Anorm := Abs(A^[1]^[1]);
331
for I := Succ(Lbound) to Ubound do
332
for J := I - 1 to Ubound do
333
Anorm := Anorm + Abs(A^[I]^[J]);
339
2: for L := N downto Succ(Lbound) do
341
S := Abs(A^[L - 1]^[L - 1]) + Abs(A^[L]^[L]);
342
if S = 0.0 then S := Anorm;
343
if Abs(A^[L]^[L - 1]) <= MACHEP * S then goto 3
349
Lambda_Re^[N] := X + T;
350
Lambda_Im^[N] := 0.0;
355
Y := A^[N - 1]^[N - 1];
356
W := A^[N]^[N - 1] * A^[N - 1]^[N];
366
Lambda_Re^[N] := X + Z;
367
Lambda_Re^[N - 1] := Lambda_Re^[N];
368
if Z <> 0.0 then Lambda_Re^[N] := X - W / Z;
369
Lambda_Im^[N] := 0.0;
370
Lambda_Im^[N - 1] := 0.0
374
Lambda_Re^[N] := X + P;
375
Lambda_Re^[N - 1] := Lambda_Re^[N];
377
Lambda_Im^[N - 1] := - Z
388
if (Its = 10) or (Its = 20) then
391
for I := Lbound to N do
392
A^[I]^[I] := A^[I]^[I] - X;
393
S := Abs(A^[N]^[N - 1]) + Abs(A^[N - 1]^[N - 2]);
396
W := - 0.4375 * Sqr(S)
399
for M := N - 2 downto L do
404
P := (R * S - W) / A^[M + 1]^[M] + A^[M]^[M + 1];
405
Q := A^[M + 1]^[M + 1] - Z - R - S;
406
R := A^[M + 2]^[M + 1];
407
S := Abs(P) + Abs(Q) + Abs(R);
411
if M = L then goto 4;
412
U := Abs(A^[M]^[M - 1]) * (Abs(Q) + Abs(R));
413
V := Abs(P) * (Abs(A^[M - 1]^[M - 1]) + Abs(Z)
414
+ Abs(A^[M + 1]^[M + 1]));
415
if U <= MACHEP * V then goto 4
417
4: for I := M + 2 to N do
419
A^[I]^[I - 2] := 0.0;
420
if I <> (M + 2) then A^[I]^[I - 3] := 0.0
422
for K := M to N - 1 do
427
Q := A^[K + 1]^[K - 1];
430
R := A^[K + 2]^[K - 1];
431
X := Abs(P) + Abs(Q) + Abs(R);
439
S := Sign(Sqrt(Sqr(P) + Sqr(Q) + Sqr(R)), P);
445
A^[K]^[K - 1] := - A^[K]^[K - 1];
449
A^[K]^[K - 1] := - S * X
459
P := A^[K]^[J] + Q * A^[K + 1]^[J];
462
P := P + R * A^[K + 2]^[J];
463
A^[K + 2]^[J] := A^[K + 2]^[J] - P * Z
465
A^[K + 1]^[J] := A^[K + 1]^[J] - P * Y;
466
A^[K]^[J] := A^[K]^[J] - P * X
468
for I := L to IMin(N, K + 3) do
470
P := X * A^[I]^[K] + Y * A^[I]^[K + 1];
473
P := P + Z * A^[I]^[K + 2];
474
A^[I]^[K + 2] := A^[I]^[K + 2] - P * R
476
A^[I]^[K + 1] := A^[I]^[K + 1] - P * Q;
477
A^[I]^[K] := A^[I]^[K] - P
488
function EigenVals(A : PMatrix; Lbound, Ubound : Integer;
489
Lambda_Re, Lambda_Im : PVector) : Integer;
491
Balance(A, Lbound, Ubound);
492
ElmHes(A, Lbound, Ubound);
493
EigenVals := Hqr(A, Lbound, Ubound, Lambda_Re, Lambda_Im);
496
procedure DivLargest(V : PVector; Lbound, Ubound : Integer;
497
var Largest : Float);
501
Largest := V^[Lbound];
502
for I := Succ(Lbound) to Ubound do
503
if Abs(V^[I]) > Abs(Largest) then
505
for I := Lbound to Ubound do
506
V^[I] := V^[I] / Largest;
509
function EigenVect(A : PMatrix; Lbound, Ubound : Integer;
510
Lambda, Tol : Float; V : PVector) : Integer;
512
procedure SetMatrix(A, A1 : PMatrix; Lbound, Ubound : Integer; Lambda : Float);
513
{ Form A1 = A - Lambda * I }
517
CopyMatrix(A1, A, Lbound, Lbound, Ubound, Ubound);
518
for I := Lbound to Ubound do
519
A1^[I]^[I] := A^[I]^[I] - Lambda;
522
function Solve(A : PMatrix; Lbound, Ubound, N : Integer;
523
Tol : Float; V : PVector) : Integer;
524
{ Solve the system A*X = 0 after fixing the N-th unknown to 1 }
528
ErrCode, I, I1, J, J1, Ubound1 : Integer;
530
Ubound1 := Pred(Ubound);
532
DimMatrix(A1, Ubound1, Ubound1);
533
DimMatrix(W, Ubound1, Ubound1);
534
DimVector(B, Ubound1);
535
DimVector(S, Ubound1);
536
DimVector(X, Ubound1);
539
for I := Lbound to Ubound do
544
for J := Lbound to Ubound do
548
A1^[I1]^[J1] := A^[I]^[J];
551
B^[I1] := - A^[I]^[J];
554
ErrCode := SV_Decomp(A1, Lbound, Ubound1, Ubound1, S, W);
558
SV_SetZero(S, Lbound, Ubound1, Tol);
559
SV_Solve(A1, S, W, B, Lbound, Ubound1, Ubound1, X);
561
{ Update eigenvector }
563
for I := Lbound to Ubound do
573
DelMatrix(A1, Ubound1, Ubound1);
574
DelMatrix(W, Ubound1, Ubound1);
575
DelVector(B, Ubound1);
576
DelVector(S, Ubound1);
577
DelVector(X, Ubound1);
582
function ZeroVector(B : PVector; Lbound, Ubound : Integer; Tol : Float) : Boolean;
583
{ Check if vector B is zero }
589
for I := Lbound to Ubound do
590
Z := Z and (Abs(B^[I]) < Tol);
594
function CheckEigenVector(A1 : PMatrix; V : PVector;
595
Lbound, Ubound : Integer; Tol : Float) : Boolean;
596
{ Check if the equation A1 * V = 0 holds }
601
DimVector(B, Ubound);
604
for I := Lbound to Ubound do
605
for K := Lbound to Ubound do
606
B^[I] := B^[I] + A1^[I]^[K] * V^[K];
608
{ Check if B is zero }
609
CheckEigenVector := ZeroVector(B, Lbound, Ubound, Tol);
611
DelVector(B, Ubound);
614
procedure Normalize(V : PVector; Lbound, Ubound : Integer);
615
{ Normalize eigenvector and make sure that the first component is >= 0 }
621
for I := Lbound to Ubound do
622
Sum := Sum + Sqr(V^[I]);
624
for I := Lbound to Ubound do
625
if V^[I] <> 0.0 then V^[I] := V^[I] / Norm;
626
if V^[Lbound] < 0.0 then
627
for I := Lbound to Ubound do
628
if V^[I] <> 0.0 then V^[I] := - V^[I];
632
ErrCode, I : Integer;
636
DimMatrix(A1, Ubound, Ubound);
638
{ Form A1 = A - Lambda * I }
639
SetMatrix(A, A1, Lbound, Ubound, Lambda);
641
{ Try to solve the system A1*V=0 by eliminating 1 equation }
644
if (Solve(A1, Lbound, Ubound, I, Tol, V) = 0) and
645
CheckEigenVector(A1, V, Lbound, Ubound, Tol)
651
until (ErrCode = 0) or (I > Ubound);
655
Normalize(V, Lbound, Ubound);
659
EigenVect := MAT_NON_CONV;
661
DelMatrix(A1, Ubound, Ubound);
664
function RootPol(Coef : PVector; Deg : Integer;
665
X_Re, X_Im : PVector) : Integer;
667
A : PMatrix; { Companion matrix }
668
N : Integer; { Size of matrix }
669
I, J, K : Integer; { Loop variables }
670
ErrCode : Integer; { Error code }
676
{ Set up the companion matrix (to save space, begin at index 0) }
678
A^[0]^[J] := - Coef^[Deg - J - 1] / Coef^[Deg];
679
for J := 0 to Pred(N) do
680
A^[J + 1]^[J] := 1.0;
682
{ The roots of the polynomial are the eigenvalues of the companion matrix }
684
ErrCode := Hqr(A, 0, N, X_Re, X_Im);
686
if ErrCode = MAT_OK then
688
{ Sort roots in increasing order of real parts }
689
for I := 0 to N - 1 do
693
for J := Succ(I) to N do
694
if X_Re^[J] < Temp then
699
FSwap(X_Re^[I], X_Re^[K]);
700
FSwap(X_Im^[I], X_Im^[K]);
703
{ Transfer roots from 0..(Deg - 1) to 1..Deg }
704
for J := N downto 0 do
706
X_Re^[J + 1] := X_Re^[J];
707
X_Im^[J + 1] := X_Im^[J];
added
removed
npm/dmath/eigen.pas
