INFORMATICA 1/1982
ALGORITHMS FOR THE SOLUTION
OF THE GENERALIZED
EIGENVALUE PROBLEM
Z. BOHTE, J. GRAD
UDK: 521.643.5 INSTITUTE OF MATHEMATICS, PHVSICS AND MECHANICS
JADRANSKA 19, LJUBLJANA, VUGOSLAVIA
ln this paper we deal with generajized eifenvaluii problem of matrix A(z) the elemcnts oi uhicli ni'
j)olynomials iu z. Tlio well known iterative metliods of Muller, lluuton and Lagaerre for finu.ini; tli^
zeros of funetion f(s') = det A(z) arv analysed. UtcomponitioiiG of matrix A(z) and its aori vat.i V'=i:
ave introducetl in ordtr to E;implify the coiuputatiouL; of the values of f(z) and its first and
seeond derivatives. Coitiparative analysi:; giveu some indicators about the rate of converiifiicu,
computer timt; and accuracy of the computed feifenvaluns for eacli of t:he mui-tiod used. Tliu cinr>lyv.'.--a
methods provcd generally to be s:;table, economical and easily applicable. FOR'1'RAII subroulines and
cjii exaniple of main calling, pro^raričme are added for Mullur*s and Lafuerre'c inothods wli.icn piovcii L
be mor'e efficient tlian Uev;ton's.
ALCORITMI ZA RnSLVANJi: TOSPLOSLNEGA PRObLJiMA I.ASTMIli VKKDNOSTI. V članku je obraviidvario miipiuriino
rvševanje posploženega problenia lastnih vruilnosti .-.a matriko A(.z) z eleiiiKfiti,- ki t;o JJO l inoini
G[<reinenl jivke z. Primurjane so dobro znane i terati vn<i metode Hullerja, Newtona in Lagur.'cra zu
računanje ničel polinoma f(z) = det A(z). Vrednosti j.iolinoma f(z) in njefovih odvoiiov ;.:(J j .•. riičunii
ne na osnovi razcepa niatriko A(z). Pi-imorjava metod dane vpof,led v hitrost konvcrfenco, [jorcib ijun
računski čas in natančnost izračunanih lastn.ih vrtdnosti za vsako metodo posobuj. Anal.i ii rauu
nit;tode so vse luimerično stabilnt-, ckonoiničnt: in enoatavno uporabne. Podprograrni v rOR'1'l'AIJ-u i n
pri.mer glavnega proj;rama so dodani za Mull^rjevo in Lague.rrovo metodo, ki :,ta Ge :izkr:/.ali aa boj i
učinkoviti od Newtonove.
Tntroduction
In this papcr three nuiuurical methods for tln;
solucion of the generalizod alj_',ehraic ei^envalue
pfoblem
A(z)x = 0 (1)
ar-e described and the programiiiefj of two of
them are prosented. Matrix A(z) is of the form
A(z) = A' •+ A.z + ... + A z
o 1 m
where all A^ a,re real squar'e
order n.
s of
The standard source of this problein is the
solution of a sy3tem of n linear homogeneous
ordinary <lifferential equations of order m.
If the matrix A^ in non-singular then tlie
problem (1) is equivalent to the classical
eigenvalu^ problom
0
0
0
Bo
I
11
0
Bl
0 ...
I ...
0 ...
B. ...
0
0
I
-1
yo
yl
ym-2
yra-l
Ci)
of order ni.n where l:. = -A A..A ;;iii,ildr
I m i
roduction is possible if A is iion-sirijni l.ir.
Although we have suveral efficient ami alai.Je
numcirical niuthodn for the solution >:>f 11it=
clnssiical cipc:nvalue j.roblem, for instaiice, <,i
alpori thni i_'lj, it is not economical to r^duoe
tlie problem (1) to the standard fonn (3).
Moreovor, niany practical pi'Cblenr; ljkc
"flutter" problem in aorodynamicE, ;iro :uiu)i
that botli A arul A are sjnpular afid ttio
o in '•
reduction to the standard forro causes fui'ttit;r
difficulties [.3].
It is therefore considered as mor>i ativantči^uu
to work directly witli the matrix (7) \n all
casi;S. Thu most naLuivil approach i:; to

44
determine the neroes of the polynomial
f(z ) = det A(z) (M)
using some of the many well known methods.
There are some efficient and stable algorithms
for tlie evaluation of the determinant of the
matrix, for instance, Gaussian elimination
with pivoting [Vf . It has been believed so far
that the evaluation of the' derivative of the
determinant involves much more work than the
evaluation of the determinant does, and
therefore the interpolation methods which
require only function values had been suggested
for the solution of the algebraic equation
f(z ) => 0 (5)
see [3J , [VJ . The most popular of these
methods is Muller's method {_HJ of successive
quadratic interpolation. The prejudice against
the methods which require the derivatives of
the function is probably based on the fact
that the derivative of the determinant is a sum
of n determinants. However, there can be found
in [l] a formula for the logarithmic derivative
of the determinant which together with its
derivative enable us to solve (5) by using
Newton's' or Laguerre's method efficiently.
In this paper algorithms for the solution of
the generalized^igenvalue problem (1), based
on Muller"s', Newton's and Laguerre"s methods
are described and compared with regard to: (i)
the number of operations (multiplications and
divisions) involved in iteration step, (ii) the
number of iterations, and (iii) the computer
time used. The computations were performed on
a number of typical exampl.es.
Muller^s method
Muller's method for the solution of the
equation (5) as described ,in [V| is the
following iterative process. From the three
previous consecutive approximations Co a root
of (5): zr_2> zr_j> zr> the parabola through
the points
<V2,f(zr_2)), 'Vi-^Vi11- (2r,f(zr>)
is formed arid its zero which is the closest fo
z^ is accepted as the next approximation, say
zr+1> to a root. The asymptotic rar.e of
convergence to a simple root in of ordtr 1.84
and at each iteration step^ only orie new
function value is required. It is necesuary to
use complex arithmetic ChroughouC.
The algorithm i3 as follows.
Fbr r=2, 3, ... we calculate
hr s 2r " zr-l' kr = hr/hr-l'
6r = fr-2kr " fr-ldr + Vkr
zr+l
- fr-ldr
(b)
The sign in the denominator of (6) is chosen ao
as to give the corr<
absolute magnitude.
as to give the correction to z the smaller
The choice of the first three approximations is
left to the user. There are no pi'oofs of the
global convergGnce of the iterative process.
In this case we have to calculate the value of
f = f(z) = det A(z)
at each iteration step what can be done most
efficiently by triangular decompouition of the
matrix A(z) using partial pivotinj; as follows.
For the given value of z we first calculate the
elements of the matrix
A s A(z) (7)
what requires m.n operations. Then we
decompose A into the forin
PA = LU (8)
where L is unit lower triangular matrix, U is
upper triangular matrix, and P is tlie
corresponding permutation matrix Q»'J . 'fhis step
requires n /3 - n/3 operations. Finally we
compute
f = ± det U = t u^u«. ... u
11 22 nn
what requires n-1 operations. Togother we must
perform
opM = m.n2 + n3/3 + 2n/3-l
operations for the calculation of f.
Newton^s methocl
The well known Newton's method for the solntion
of (5) is the iterative process'
At each iteration ntop we evaluate ttie
corriiction of the current approxiraation to a
root as the reciprocal of the logar'i thinic
derivative of the function. From |l|, | 2*1 wu
have the formula
f'(z)/f(z) = triA'* (z)l\'(z)) - tr X (0)

45
where A'(z) denotes the derivative of the
matrix A(z) with respect to z.
An efficent algorithm for calculating the
expression (9) is as follows.
First we calculate the matrices
A = A(z), B = A'(z) (10)
• 2
which requires (2m-l)n operations. Then we
decompose A as in (8) obtaining
PA = LU (11)
which requires n /3 - n/3 operations, compute
the matrices Y and X which satisfy
Vi = PB (12)
UX = Y . (13)
As we need only the trace of X, we can use
formulae
fori=n, n-1,
. j+1 ij
andk=l,2, i. These
two steps require together (n /2 - n /2) +
3 • 2 '
+ (n /6 + n /2 + n/3) bperations. So we can
calculate
f(z)/f(z) =
..
with
= (2m -
(15)
(16)
operations.
Although the asymptotic rate of convergence of
the Newton's method to a simple root is 2 there
seem to be little advantage in using Newton*s
method against Muller's..There appear to be
almost three times as many operations per
iteration step in Newton's method than in
Muller's.
S^ -- t' (zr)/f(zr)
S2 = (f (zr)2 - f(zr)f"(zr))/f(zr)2
The sign in the denominatpr is chosen so as to
give the correction to zr the smaller absolute
magnitude. The form of S^ which is suitable
for numerical evaluation iB defined by (9) and
(10) - (1Š), while the formula for S2 is
derived from (9) by differentiation, see 1.2] ,
Sj = ^((A^U^A^z^,))7) - tMA-^z^A"!^))
(17)
The algorithm for evaluation of S_ demands in
addition to the starting operations for S^,
i.e. (10) and (11), the following computations:
We must first calculate the matrix
C = A"(z) (18)
9
with additional (m-2)n operations. Then we
calculate the trace of the matrix W = A C by
steps . \ •
LZ = PC : (19)
UW = Z .
tr W
(2 0.)
(21)
3 2 '
which require (n /2 - n /2) +
3 2
+ (n /6 + n /2+n/3) + (n - 1) operations.
The first term in (17) is found by calculating
first the full matrix X by (12) and (13) with
3 •' 2
n operations and then the trace of X as
tr«') =
x, .
(22)
2
which requires additional n operations.
At each iteration step we must perform
= (3m - 3)n 2n
operations for the calculation of S- and S«.
This number is almost twice as large as (16)
for the Newton's method.
Laguerre^s method
The quite famous Laguerre's mcthod requires
apart from the function value also tlie values
of its first two derivatives. The aBymptotic
rate of convergence to a simple root is 3.
Perhaps the most attractive feature from the
users point of view is its stability in global
convergence for any value of thu firnt
approximation.
The iterative process is defined Iy [M|
zr+l = zr " n/(Sl * ((n"1)tn
for r=l,2, ..., where
'2
Numerical tests -
The programmes in FORTRAN programming language
were written for all threa algorithms and
tested on CDC CYBLR 72 computer. Single
precision complex arithmetic wds used with
t=M8, where t is the number of significant
digits in the mantissa of a binary floating-
point number. The iterative process was stopped
after either the cotnputed correction of the
computed approximation to the eigenvalue was
less than or equal to e or ^^|uiiI — e> ^or
1 < i < n, occured in (8) or (11). The value
for e was chosen as 10.1! A(z)|| „. 2~ .

The implicit deflation of f(z) was applied in
the programmes in order to prevent the
iterative process to converge to one of the
previously computed eigenvalues, say
x., x_, ..., x . This introduces the following
modifications of algorithms:
A(z) =
f(z) into f(z)/ n (z-x
i
f'(z)/f(z) into f (z)/f(z)-
into
into
m
.- l
(z-x.)
S,-
2
l (z-x.)
i=l ^
i'
m
= 1
-1
-2
(23)
.-1
and (24)
(25)
Some typioal examples of (1) were computed
with multiple eigenvalues either of finite
value or in the infinity. The values of the
parameters m and n of A(z) were on intervals
1 < m < 2 and 3 < n < 12 respectively. No
breakdown of the iterative process occured
regardless the algorithm used. All three
methods proved to be stable and accurate. It
was found that the previous computed
.eigenvalue and the values round it were the
most suitable starting values for computing
the next eigenvalue.
In the table
Muller
No T
Newton
No T
Laguerre
No T
1
2
3
4
5
6
7
8
55
0
11
0
3
0
1G
1
.57
.13
.04
.17
. 02
41
<f
19
0
4
0
41
2
.58
.28
.06
.61
.04
17
0
11
0
2
0
3
2
.45
.30
.06
.08
. 06
ALL 86 .93 107 1.57. 35 . 95
the number of iteration steps (No) and the
computer time used CT jn seo) are shqwn for
each of the computed eigenvalue z., i=l,...
of the system
.8,
+ zA. + z - 0
where I, and
[~0 -3 0
2 0 0
0 2 0
0 2
•
V
-10 0
0. .-2 0
0 0
1 0
The system has a triple eigenvalue at ± i and
a double eigenvalue at 0 [5]. Somehow similar
results were obtained with other examples.
1
15
1
-1
0
0
1
0
1
• V
-2 1 -1
3 0 1
1 0.5 0
1.00000000000
2.08866333896
-0.256555796702 ± 0.896010203022 i
-2.91609433069
11.3405425851
2) A(z) =
V
V
-1
-3
-3
-3
-3
-3
1
1
1
1
1
1
-3
-4
-3.1
-3.1
-3.1
-3.1
1
0
0
0
0
0
-3
-3.1
2.8
3.8
3.8
3.8
1
0
-1
-1
-1
-1
-3
-3.1
3.8
9.8
10.7
10.7
1
0
-1
-2
-2
-2
-3
-3.1
3.8
• 10.7
12.6
14.6
1
0
-1
-2
-3
-3
-3
-3.1
3.8
10.7
14.6
15._6_
1
0
-1
-2
-3
-2
0.908770404173 ± 1.93967680102 i
0.931536974557 ± 1.97197662562 i
4.18245919165
6.13692605089
3) A(z)
V
= Ao
~1
-1
0
0
0
0
0
+ AlZ
0
• 1
-1
0
0
0
0
+ A2z
0
0
3
-4
0
-8
0
0
0
0
7
-1.5
0
0
0
1
0
0
1
0
0
0
0
0
0
0
10
-0.4
1
0
0
0
0
0
1
Three more examples are

47
V-
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
0
0
0
0
0
0
0
0
0
0
2
0
0
-1
0
0
0
1
0
0
0.
0
0
0
0
6
0
: 0
0
0
0
0
5
. 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
9
0
0
. 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
ProRramrr.es
1. Comments
In this chapter two subroutines EIGENM and
EIGENL for computing the eigenvalues of a
general matrix (2) are presented. An example
of a possible main (calling) programme for both
subroutines is also given at the end of the
chapter. In EIGENM the eigenvalues are computed
by Muller*s method while in EIGENL by
Laguerre's method. If A(z) in (2) satisfies
m < 3 then the subroutines can be used such as
they are presented. For m > 3 some changes in
the subroutines are necessary as explained at
the end of the prograrames.
-1.09579878411
-1.08865049880 ± 1.04888469390 i
-0.541227886923 ±1.51715942168 i
z6 - z14 a» K . 107, |K| »1
The last nine eigenvalues lie in the infinity.
Muller's and Laguerre's methods were also
tested on DEC 1091 computer with t=27. In one
case w)ien some of the elements of A(z)
differed up to 10 in their absolute values
Muller's method failed to compute correctly
some of the. eigenvalues while all the
eigenvalues obtained by Laguerre's method were
correct.
References
1. BODEVJIG, E. Hatrix Calculus. North-Holland
Publishing Company, Amsterdarn, 1959.
.2. B0HTE, Z. On Algorithms for the Calculation
of Derivatives of the Determinant. Proc.
• ALGORITHMS 79 Symposium on Algorithms,
Strbske Pleso, Czechoslovakia, 1979, pp.
;nr-iio.
3. PETERS, G., AND WILKINSON, J.H. Ax =
>Bx and the Generalized Eifenproblem. SIAM
J. Numer. Anal. 7, 4(1970), 479-492.
4. W:LKINSON, J.H. The Algobraic Ligunvalue
Problem, Clarendon Pross, Oxford, 1965.
5. LAHCASTER, P. Lanibda-matriceu and Vibrating
Systems. Pergamon Press, Oxford, London,
Edinburgh, New York, Toronto, Paris,
. Braunschweig, 1966.
The subroutines LIGEHM and EIGFNL are
completely selfcontained (EIGENM is composed
of four subroutines EIGENM, F, ALAMDA, LU and
EIGENL is composed of ten subroutines EIGEML,
ALAMDA, LU, DERIVl, DERIV2, L0WER, UPPER,
UPPER2, TRACE, TRX5Q) and communication to them
is solely through their argument lists and .
C0MM0N statements. The entrace to the
subroutine EIGENM is achieved by
C0MM0N /HATRIX/ AO(ND.ND), AHND.ND), A2(ND,
1ND)
CALL EIGENM (M, N, ND, T, INDIC, LAM, PERMUT,
1A, L, U)
The entrance to the subroutine EIGENL is
achieved by
C0MM0N /MATRIX/ AOCND, ND), AKND.ND), A2CND,
1ND), A3(ND,ND)
CALL EIGINL (M, N, ND, T, INDIC, LAM, PERMUT,
1A, D, L, U) -
The nieaning of the parameters is described in>.
the following lines.
The elementE of. matrices A ,A. ,A», ..., A in
(2) are to be stored in tbe first N rowe and
columns of the two diinensional arrays A0, Al,
A2, A3,... .
M, vihere M > 1, is the degree of ž-matrix
A(z) in (2).
N, wh«:re K ; 1, is the order of matrices
Ao' V A2 AM- . •
ND, where ND > N, defines the first dimension
of the two dimensional arrays A0, Al, A2 ,..
.., A,D,L,U and the dimenaion of tlie one
dimensional array PERMUT. Dimension of the
one di.niensional arrays IMUIC and LAM must
be equal to or greatar than ND*M.
'1' is the number of biriary digits in the

48
mantissa of a single precision floating-
point number. T is of integer type.
The array INDIC indicates the sucoess of
the subroutine as:follows
value of INDIC(I) e.igenvalue (I)
.0 not found after 100
iteration steps
1 found
LAM is one dimensional complex array. The
computed eigenvalues will be found in the
first N*M places.
The arrays PERMUT, A, D, L, U is the working
storage.
The main (calling) programme should contain
the following declaration statements
INTEGER, M,N,T,INDIC,PERMUT '
REAL A0,Al,A2,A3
C0MPLEX LAM,A,L,U
C0MM0H /MATRIX/ AO(ND.ND), AKND.ND),
1A2(HD,ND),'A3(ND,ND).
DIMENSI0N INDIC (ND*M),LAM(ND*M).PERMUT
l(fft)),A(ND,ND+l),L<ND,ND+l),U(ND,ND+l)
EQUIVALENCE
and the following two additional declaration
statements for the subroutine EIGENL
C0MPLEX D
DIMENSI0N DCND.ND)
Note that.the dimensions of the arrays in
DIMENSI0N and C0MM0N statements of the main
programme and in C0MM0N sjatements within
subroutines must be constants of integer type.
See the example in 5 where ND was replaced by
12, ND+1 by 13 and ND*M by 36.
2. Subroutine EIGENM - Muller's method
The subroutine is oomposed of four subroutines
EIGENM, F, ALAMDA and LU.
SUBR0UTINE EIGENM(M,N,ND,T,INDIC,LAM,
1PERMUT,A,L,U)
INTEGER M,N,ND,T,INDIC,fERMUT,I,IND,
1ITERST, NM, R
REAL EPS.EPS100, IM, RE
REAL AO,A1,A2,A3
C0MPLEX LAM,A5L,U,DEN0M1,DEN0M2,DR,FR,
C0MM0N /MATRIX/ A0(ND,ND),Al(ND.ND),
1A2(ND,ND),A3(ND,ND)
DIMENSI0N INDIC(1),LAM(1).PERMUT(1),
1A(ND,1),L(ND,1),U(ND,1)
NM = N*M
ITERST = 102
D0 1 1 = 1, NM
INDIC(I) = 0
1 C0NTINUE
1 = 1
C THE MULLER ITERATIVE PR0CESS.
C DETERMINATI0N 0F THREE STARTING
C C0NSECUTIVE APPR0XIMATI0NS T(3 THE
C C0MPUTED EIGENVALUE LAM(I),F0R
C 1=1,2,...NM.
2 ZR = (l.O,O.O)*FL0AT(I)
3 ZRMl = 0.5*ZR
ZRM2 = 0.1*ZR
C THE INITIAL C0MPUTATI0NS 0F THE NEXT
C APPR0XIMATI0N Z(R+1) T0 LAM(I),F0R R=2.
R = 2
CALL FCZRM^.FRM^.EPS.I.IND.M.N.ND.T.LAM,,
1PERMUT,A,L,U)
IF(IND.EQ.1)G0 T0 5
CALL F(ZRM1,FRM1,EPS,I,IND,M,N,ND,T,LAM,
1PERMUT,A,L,U)
IF(IND.EQ.1)G0 T0 6
QALL F(ZR,FR,EPS,I,IND,M,N,ND,T,LAH,
1PERMUT,A,L,U)
IF(IND.EQ.1)G0 T0 7
C
HRMl = ZRM1-ZRM2
HR = ZR -ZRMl
KR = HR/HRMl
C THE MAIN L00P F0R THE C0MPUTATI0N 0F
C THE NEXT APPR0XIMATI0N Z(R+l) T0 LAM(I),
C F0R R=2,3,....ITERST.
i* DR=1.0 + KR
GR=FRM2*KR**2 - FRM1*DR**2 + FR*(KR+DR)
DEN0M2 = CSQRT(GR**2-H.*FR*DR*KR*
1 <FRM2''KR-FRM1*DR+FR))
DEN0M1 = GR+DEN0M2
DEN0M2 = GR-DEN0M2
IF(CABS C DEN0M2).GT.CABS(DEN0M1>)DEN0M1
1 = DEN0M2
KR = -2.O*FR*DR/DEN0M1
HR = HR*KR
GR = ZR
ZR = ZR+HR
FRM2 = FRMl
FRHl = FR
CALL F(ZR,FR,EPS,I,IND,M,N,ND,T,LAM,
1 PLRMUT.A.L.U)
C
IF(IMD.EQ.1)G0 T0 7
IF(CABS(ZR-GR).LE.1O.*EPS)G0 T0 7

C
C
c
PA = LU
WHERE L IS UNIT L0WER TRIANGULAR MATRIX,
. U IS UPPER TRIAMGULAR MATRIX AHD P IS
THE C0RRESP0NDING PERMUTATI0N MATRIX
<SEE WILKINS0N).
CALCULATI0N 0F A(Z.).
CALL ALAMDA (Z,M,N,ND,A)
CALCULATI0N 0F THE SMALL P0SITIVE
NUMBER EPS.
C = (0. ,0.)
c
c
c
. R = R+l
. IF(R.LT.ITERST)G0 T8 t
LAMCI) = ZR
G0 T0 9
5 2R = ZRM2
G0 T0 7 C
6 ZR = ZRMl C
7 LAM(I) = ZR
INDIC(I) =1
1=1+1 C
C0MPUTE THE C0MPLEX C0NJUGATE VALUE 0F
LAM(I).
RE = REAL(ZR)
IM = AIMAG(ZR)
EPSIOO = 10O.O*EPS .
IF<ABS(I;M).LE.EPS1OO)G0 T0 8
LAM(I) = C0NJG(ZR)
INDIC(I) =1
I =1 + 1
8 IFCNM.LT.DG0 T0' 9
DETERMINATI0N 0F THE STARTING APPR0XIM
ATI0N T0 THE C0MPUTED EIGENVALUE LAM
(1+1). .
ZR^ = 0.9*ZR
. RE = ABS(RE)+ABS(IM)
IF(RE.LE.EPS1OO.0R.1./RE.LE.EPS1OO)G0
1T0 2 .
G0 T0 3
9 C0NTINUE
. RETURN
END . •
SUBR0UTINE F<Z,FUNC,EPS,1,IND,M,N,ND,T,
1LAM,PERMUT,A,L,U)
INTEGER I,IND,M,N,ND,T,PERMUT,IM1,J,K
REAL EPS.FK C
C0MPLEX Z,FUNC,LAM,A,L,U,C C
DIMENSI0N LAM(1),PERMUT(1),A(ND,1), C
1L(ND,1),U(ND,1) C
C0MPUTATI0N 0F FUNC = FCZ) = DET A(Z) »
DET A. MATRIX A IS DEC0MP0SED INT0 THE
D0 1 J=1,N .
D0 1 K=1,N
C = C + A(J,K)**2
1 C0NTINUE
EPS = CABS(CSQRT(C))*2.0**(-T)
C0MPUTATI0N 0F THE MATRICES L AND U,
WHERE PA = L*U.
CALL LU(EPS, IND.N^ND.PERMUT.A.L.U)
IF<IND.EQ.1)G0 T0 4
J = -IND
K = (-1)**J
FK = FL0ATCK)
FUNC = CMPLX(FK,0.0)
D0 2 K=1,N
FUNC = FUNC*U(K,K)
2 C0NTINUE
M0DIFICATI0N 0F.FUNC = F(Z) INT0 F(Z)/
((Z-X(1))(Z-X(2))...(Z-X(I-D) ACC0RDING
T0 EQUATI0N (23) 0F.THIS PAPER.
IMl = 1-1
IF(IM1.EQ.O)G0 T0 4
D0 3 J=1,IM1 '
FUNC = FUNC/(Z-LAM(J))
3 C0NTINUE . . •
4 C0NTINUE
RETURN
END '.'.;•'
SUBRflUTINE ALAMDA(Z,M,N,ND,A)
INTEGER M,N,ND,I,J,K
REAL A0,-Al,A2,A3
C0MPLEX Z.A.AIJ
C0MPLEX Z2..Z3
C0MM0N /MATRIX/ AO(ND.ND),A1(NU,ND),
1A2(HD,ND) ,A3(tlD,ND) . •
DIMENSI0N A(MD,l)
THIS SUBR0UTINE CALCULATES THE ELEMENTS
0F MATRIX A(Z), VIHICH ARE FUNCTI0NS 0F Z,
F0R A GIVEN VALUE 0F Z.
Z2 = Z**2
Z3 = Z**3
C0NTINUE . , • • '
D0 4"» 1 = 1,N
D0 4t J=1,M . "
AIJ = AO(I.J)
D0 3 3 K=1,M . . . :
G0 T0 <1,2,3),K
1 AIJ = AIJ+Ai(I,J)*Z
G0 T0 3 3
2 AIJ = AIJ+A2(I,J)*Z2
G0 T0 3 3
3 AIJ = AIJ+A3(I,J)*Z3
33 C0NTINUE
A(I,J) = AIJ .

50
•m C0NTIHUE
RETURN
END.
SUBR0UTINE LU (EPS,IND,N,ND,PERMUT,A,L,
1U)
INTEGER IND,N,ND,PERMUT,I,K,NP1,R,RC,
lRMl.RPi.T
REAL EPS.SRC.STABS
C0MPLEX A,L,U,ST,URR
DIMENSI0N PERMUT(1),A(ND,1),L(ND,1 ),
1U(ND,1)
C
C THIS SUBR0UTINE PERF0RMS THE
C DEC0MP0SITION 0F MATRIX (PA) IHT0 L*U,
C WHERE L IS UNIT L0WER TRIANGULAR MATRIX,
C U IS UPPER TRIANGULAR MATRIX AND P IS
C THE C0RRESP0NDING PERMUTATI0N MATRIX.
C PARTIAL PIV0TING IS INTR0DUCED (SEE J.H.
C WILKINSON, PAGE 225).
IND = 0
D0 1 I = 1.,N
PERMUT(I) = I
1 C0NTINUE
C
NPl = N+l
D0 11 R=1,N
RMl = R-l
RPl = R+l
RC = R
SRC = 0
D0 4 T=R,N
ST s A(T,R)
IF(RM1.EQ.O)G0 T0 3
C
D0 2 K=1,RM1
ST = ST-L(T,K)*U(K,R)
2 C0NTINUE
3 A(T,NP1) s ST
STABS s CABS(ST)
IF(SRC.GL.STABS)G0 T0 U
SRC = STABS
RC » T
H C0NTINUE
IF(SRC.GT.I:PS)G0 T0 s
IND = 1
RETURN
C
5 IF(R.EQ.RC)G0 T0 7
T = PERMUT(R)
PERMUT(R) = PLFHUT(RC)
PERMUT(RC) = T
C THE F0LL01V1NG STATl.Mi:NT( IHD» I MD-1 ) MAV
C BE TAKEN 0UT IN EIGEML. '
IND = IMD-1
D0 6 Ksl,NP1
ST = A(R,K)
A(R,K) = A(RC.K)
A(RC,K) = ST
6 C0NTINUE
7 URR = A(R,MP1)
U(R,R) = URR
IF(R.EQ.N)G0 T0 11
D0 10 T=RP1,N
L(T,R) s A(T,NP1)/URR
ST = A(R,T)
IF(RM1.EQ.O)G0 T0 9
8
9
10
11
D0 8 K=l
ST =
C0NTINUK
U(R,T) =
C0NTINUE
C0NTINUE
RETURN
END
,RM1
ST - L(R,K)*U(K,T)
ST
3. Subroutine EIGENL - Laguerre's method
The subroutine ic composed of ten subroutines
EIGENL.ALAMDA (see 2.), LU(see 2.), UERIVl,
DERIV2, L0WER, UPPER, UPPER2, TRACE and TRXSQ.
SUBR0UTINE EIGENL(M,N,ND,T,IND1C,LAM,
1PERMUT,A,D,L,U)
INTECER M,N,ND,T,INDlC,PtRMUT,I,IM1,1ND,
1ITERST,J,K,NM,NM1,R
REAL EPS,LPS100,IM,RE
REAL A0.A1,A2,A3
C0MPLEX LAM , A , D, L ,U, CEPS , DEN0M1, D1.N0M2 ,
1S1,S2,TRAC1;X,TRACEW,TRACX2,ZR,ZRP1
C0MM0N /MATRIX/ A0(ND,ND),A1(NU.ND)
1A2(ND,ND),A3(ND,ND)
DIMENSI0N ItlDIC(l) ,L,AM(1) .PERMUTCl),
1A(ND,1),D(ND,1),L(ND,1),U(ND,1)
C
NM s N*M
ITERST =101
00 1 1=1,NM
INDIC(I)=0
1 C0NTINUE
1 = 1
THE LAGUERRE ITERATIVE PR0CLSS.
DETF.RMIHATI0N 0T TME STARTING
APFR0XIMATI0N T0 TIIE FIRST C»MPUTi;U
EIGENVALUl: l,AM(l).
ZR = (1.,1.)
2 NMl = NH-I
C0MPUTAT10U 0F TIIE NEXT LIGEMVALUt LAM(I).
R = 1
C
C TIIE HMM L00P F0R Thl. C0MPUTATI0II TliL

51
C . NEXT APPR0XIMATI0N Z(R+l5 T0 LAM(I).
C CALCULATI0N 0F THE EQUATI0NS (10) - (22)
C 0F THIŠ PAPER.
C CALCULATI0N 0F A(Ž). ' C
3 CALL ALAMDA(ZR,M,N,ND,A)
C CALCULATI0N 0F THE SMALL P0SITIVE NUMBER
C EPS.
CEPS =0 C
D0 t J=1,N C
DO.it K=1,N . C
CEPS = CEPS+A(J,K)**2 C
'"» C0NTINUE •
EPS = CABS(CSQRT(CEPS))*2.0**(-T)
EPS100 =.1OO.*EPS
C C0MPUTATI0N 0F THE MATRICES L AND U,
C WHERE PA '= L*U. ; .
CALL LU(EPS,IND,N,ND,PERMUT,A,L;U)
IF(IND.EQ.1)G0 T0 8
C C0MPUTATI0N 0FTHE TRAGE 0F MATRIX X
C AND THE TRACE 0F X**2.
CALL DERIV1(ZR,M,N,ND,PERMUT,D)
CALL L0WER(N,ND,D,L)
CALL UPPER(N,ND,D,U)
• CALL TRACE(N,ND,D,TRACEX)
CALL TRX9Q(N,ND,D,TRACX2) C
C C0MPUTATI0N 0F THE TRACE 0F MATRIX W.
TRACEW =0.0
1F(M.LT.2)G0-T0 5
CALL DERIV2(ZR,M,N,ND,PERMUT,D)
CALL L0WER(N,ND,D,L)
CALL UPPER2 (N.ND.D.U)
CALL'TRACE(N,ND,D,TRACEW)
C ••• • •
5 Sl = TRACEX
S2 = TRACX2-TRACEV/
C ' M0DIFICATI0MS 0F Sl AND S2 ACC0RDING
,C T0 THE EQUATI0NS (2^) AND (2S) 0F THIS
C . PAPER.
' • IMl - 1-1 c
IF(IM1.EQ.O)G0 T0 7 C
D0 6 J=1,IH1 . • c
• CEPS = l.d/(ZP-LAM(J)) C
. Sl ;=•• Sl-CEPS • " .
S2 = S2-CEPS**2
6 C0NTINUE
C ' ' '
7 CEPS = CSQRT(MM1*((HM1+1)*S2-S1**2))
DEN0M1- = Sl + CEPS
•; DEN0M2 = Sl-CEPS
IF(CABS(DEM0Ml).GT.CAbS(DEK0H2))Dl:M0M2= .
1DEN0M1
ZRPl = ZR-(NH1+1)/DEK0H2
C
IF(CABS(ZRP1-ŽR).L1:.1O.*LPS)C0 Td) 9
ZR = ZPPl
R = R + l
IF(R.LT.ITERST)G0 T0 3
LAM(I) - ZR
G0 T0 11 • • ' • .•
S ZRPl = ZR
9 LAM(I) = ZRPl
INDIC(I) .= 1 ' '
DLTERMINATI0N 0F THE STARTIiJG
APPR0XIMATI0N T0 THE NEXT C0MPUTED
EIGENVALUE LAM(I.+ 1) AND C0MPUTATI0N 0F
THE C0MPLEX C0NJUGATE VALUL 0F LAM(I).
ZR = O.9*ZRP1
1 = 1+1
RE = REAL(ZRPl)
IM =AIMAG(ZRP1)
IF(ABS(IM).LE.1.0E-3*ABS(RE))ZR=0.5*
1CMPLX(RE,RE)
RE = ABS(RE)+ABS(IM)
IF(RE.LE.EPS1OO.0R.1./RE.LE.EPS1OO)
1ZR=(1.,1.)*FL0AT(I)
IF(ABS(IM).LE.EPS1OO)G0 T0 10
LAM(I) r C0NJG(ZRP1)
INDIC(I) - 1
I = 1+1
10 IF(I.LE.NH)G0 T0.2
11 C0NTINUE
RETURN
END
SUBR0UTINE DXRIVKZ,M,N,ND,PERMUT,A)
INTEGER M,N,ND,PERMUT,I,II,J,K,MM1
REAL AO,A1,A2,A3
C0MPLEX Z.A.AIJ
C0HPLEX Z2.Z3
C0MM0N /MATRIX/ AO(HD.ND),A1(ND,ND)
DIMENSI0N A(ND,l),PERMUT(l)
THIS SUBR0UTINE CALCULATES THE ELEMENTS
0F THL FIRST DERIVATIVE 0F MATRIX A(Z),
F0R A GIVEN VALUE 0F Z.
Z2 = 2.*Z ' . .
Z3 = 3.*Z**2
C0NTINUE
MMl = M-l
D0 5S 11=1,N ,
I = PERMUT(II)
D0 55 J=1,N
AIJ a -AKI.J)
IF(M.LT.2)G0,T0 HM
D0 3 3 K=1,MM1
G0 T0 (2,3),K
AIJ = AIJ+A2(I,J)*Z2
G0 T0 3 3
AIJ = AIJ^A3(I,J)*Z3

52
c
c
c
c
C
C
C
C
C
33 C0NTINUE
44 A(II,J)=AIJ
55 C0NTINUE ' • '
RETURN
END • . • '
SUBR0UTINE DERIV2(Z,M,N,ND,PERMUT,A)
INTEGER M,N,ND,PERMUT,I,II,J,K,MM2
REAL AO,A1,A2,A3 . -
C0MPLEX Z.A.AIJ
• . C0MPLEX Z2.Z3
C0MM0N /MATRIX/ AO(ND.ND) jAKND.ND) ,
1A2(ND,ND),A3(ND,ND) :
DIMENSI0N A(ND,1),PERMUT.(1)
THIS SUB'R0UTINE CALCULATES THE ELEMENTS
0F THE SEC0ND DERIVATIVE 0F MATRIX A(Z),
F0R A GIVEN VALUE 0F Z. .
• Z3 = 6.*Z.
C0NTINUE
MM2 = M-2
D0 55 11=1,N
I = PERMUTUI)
D0 55 Jsl.N
AIJ = 2.*A2(I,J)
IF(M.LT.3)G0 T0 44 . .
Vti 33 K=1,MM2
G0 T0 (3,3),K
3 AIJ = AIJ+A3(I,J)*Z3
33 C0NTINUE
44 A(II,J) = AIJ
55 C0NTINUE
RETURN
END
SUBR0UTINE L0WER(N,ND,A,L)
INTEGER N,ND,I,IM1,J,K . ;
C0MPLEX,.A,L,AIJ
DIMENSI0N A(ND.l),L(ND,l)
THIS SUBR0UTJNE CALCULATES THE ELEMENTS
0F MATRIX M WHICH SATISFIES THE RELATI0N
L*M=A, WHERE L IS UNIT L0WER TRIANGULAR
MATRIX. M IS ST0RED INT0. ARRAV A.
D0 2 1 = 2,N . .
IMl s 1-1
D0 2 J=1,N
AIJ = (0. ,0.)
D0 1 K=l,tMl
AIJ = AIJ+L(l',K)*A(K,J)
1 C0NTINUE
A(I,J)=A(I,J)-AIJ
2 C0NTINUE
RETURN
END
C
C
C
C
C
SUBR0UTINE UPPER(N,ND,A,U)
INTEGER N,ND,I,II,IP1,J,K
C0MPLEX A.U.AIJ
DIMENSI0N A(ND,l),U(ND,D
• THIS SUBR0UTINE CALCULATES THE ELEMENTS
0F MATRIX M WHICH SATISFIES THE RELATI0N
U*M=A, WHERE U IS UPPER TRIANGULAR
MATRIX. M IS ST0RED INT0 ARRAY A.
D0 3 11=1,N
I = N+l-II
IPl = 1+1
D0 3 J=1,N
AIJ = (0. ,0.)
IFU.EQ.N)G0 T0 2
D0 1 K=IP1,N
AIJ = AIJ+U(I,K)*A(K,J)
1 C0NTINUE
2 A(I,J) = (A(I,J)-AIJ)/U(I,I)
3 C0NTINUE
RETURN
END
SUBR0UTINE TRACE(N,ND,A,TR)
INTEGER N.ND.I
C0MPLEX A,TR
DIMENSI0N A(ND,D
THIS SUBR0UTINE CALCULATES THE TRACE
0F MATRIX A.
TR = (0. ,0. )
D0 1 1=1,N
TR = TR+A(I,I)
1 C0NTINUE
RETURN
END
SUBR0UTINE TRXSQ(N,ND,A,TR)
INTEGER N,ND,I,K
C0MPLEX A,TR
DIMENSI0N A(HD,D
THIS SUBR0UTINE CALCULATES THE TRACE 0F
A**2.
TR = (0.,0.)
D0 1 1=1,N
D0 1 K=1,N
TR = TR+A(I,K)*A(K,I)
1 C0NTINUE
RETURN
END
SUBP.flUTIIir. UPPEP2 (N,ND,A,U)
IHTEGER N,ND,I,II,IP1,J,K
C0MPI,CX A,U,AI.I
DIHEMSI0H A(ND,l),U(ND,l)

53
C
C THISSUBR0UTINE CALCULATES THE
C ELEMENTS OF THE L0WER TRIANGULAR 0F
C. '. MATRIX M. M SATISFIES THE RELATI0N U*M=A,
C ,WHERE U IS UPPER TRIANGULAR MATRIX. M
C IS ST0RED IHT0 ARRAY A.
D0 3 11=1,N
I = N+.l-II
IPl =1+1
D0 3 J=1,I
AIJ = (0.,0. )
IFU.EQ.N)G0 T0 2
D0 1 K=IP1,N
AIJ = AIJ+U(I,K)*A(K,J)
1 . C0NTINUE
2 A(I,J)=(A(I,J)-AIJ)/U(I,I)
3 C0NTINUE
RETURN
END
4. Modifications
For matrix A(z) satisfying (2), where m>3,
some of the statements must be replaced with
other stat.enients as shown below for the case
1) Statements
REAL AO,A1,A2,A3
C0MM0N /MATRIX/ A0(ND.ND),Al(ND.ND),
1A2(ND,HD),A3(ND,MD)
in the main programme, EIGKNM, ALAMDA, LIGLNL,
DERIVl and DERIV2 by
REAL AO,A1,A2,A3,A4,A5
C0MM0H /MATRIX/ A0(ND,WD),A1(ND,ND),
; 1A2(ND,HD),A3(ND,ND),A4(ND.HD),A5(ND.ND)
2) Statement C0HPLi:X Z2.Z3 in ALAMDA.DERIVl
and DERIV2 by
C0MPLEX 7.2,2 3 ',?M ,Z5
3) Statemcnt C0NTINUL
(i) in ALAMDA by
ZM = Z**4
25 = 7,**b
(ii) in DERIVl by ••
ZH = H.*2**3
Z5 = 5.*Z**4
(iii) in DERIV2 by
ZM = 4 . *3 . *Z**2'
Z5 = 5. *i\. ftZ**3
^) Statement G0 T0(1,2,3),K in ALAMDA by
G0 T0(1,2.,3,M ,5) ,K
5) Statement G0 T0(2,3),K in DLKIVl by
G0 T0C2,3,4,5),K
6) Statement G0 T0(3,3),K in DLRIV2 by
G0 T0(3,4,5),K
7) Statement 33 C0IJTINUE in ALAMDA, DERIVl aiul
DERIV2 by
G0 T0 3 3
. i* AIJ = AIJ+A4(IjJ)*Z4
G0 T0 3 3
5 AIJ = AIJ+A5CI,J)ftZ5
33 C0NTINUE
5. Main - calling programme
In the following liries an example of a
possible calling programme for subroutines
EIGENH and EIGENL is presented. It finds the
eigenvalues of A(z) in (2) with ms3 and nsl2.
The programme stops after reading value zero
for n.
INTEGER INDIC,PERMUT,M,N,T,I,J,K,MN,MP1
REAL AO,A1,A2,A3
C0MPLEX LAM,A,L,U
C0MM0N /MATRIX/ A0(12,12),A1(12,12) ,
1A2(12,12),A3(12,12)
EQUIVALEMCE (A(l),L(1),U(1))
DIMENSI0N INDIC(36),LAM(36),PERMUT(12) ,
1A(12,13),L(12,13),U(12,13)
C0MPLEX D
• DIMENSI0N D(12,12)
C
1 READ 2,N,M,T
2 F0RHATOI5)
3 F0RMAT(8F1O.O)
IF(N.EQ.O)G0 T0 12
C THE PR0GRAMME ST0PS AFTER READING ZER0
C F0R N.
MIJ = M*N
MPl = M+l
D0 B 1=1,MPl
G0 T0 (4,5,6,7),I
H READ 3,((A0(K,J),J=l,tO,K=l,N)
G0 T0 8
5 READ 3,((A1(K,J),J=1,N),K=1,M)
G0 T0 8
6 REAO 3,((A2(K,J),J=1,N),K=1,N)
G0 T0 8
7 RLAD 3 ,((A3(K,J),J = 1,H),K = 1,N)
8 C0MTIKU1:

54
PRINT 9,((A0(K,J),J=l,N),K=l,N)
9 F0RMAT (1H .6E16.7)
PRINT 9,((A1(K,J),J = 1,N)-,K=1,N)
IF(M.GT.1)PRINT 9,((A2<K,J),J=l,N),
1K=1,N)
IF(M,6T.2)PRINT 9,<(A3(K,J) ,J = 1,N),
1K=1,N>
CALL EIGENM(M,N,12,T,INDIC,LAM,Pi:RMUT,
1A,L,U)
PRINT 1
10 F0RMAT(1HO,1OX,6H INDIC,10X,9HREAL PART,
112X,10H IMAG,PART,///(1H ,12X,I2,5X,
lEl-9.12,fX,E19.12))
D0 11 1=1,MN
LAM(I) =0.0
11 C0NTINUE
CALL EIGENL(M,N,12,T,INDIC,LAM,PERMUT,
1A,D,L,U)
PRINT10,(INDIC(I),LAM(I),I=l,MN)
G0 T0 1
12 C0NTINUE
END
CENIK OGLASOV
Ovitek - notranja stran (za lethik 1982)
2 stran 28 .000 din
3 stran 21.000 din
Vmesne strani (za letnik 1982)
1/1 stran 13.000 din
1/2 strani T 9.000 din
Vmesne strani za posamezno številko
1/1 stran 5.000 din
1/2 strani 3.300 dln
Oglasi o potrebah po kadrih (za posamezno številko)
2.000 din
Razen oglasov v klasičnl obliM so zaželjene tudi krajše
poslovne, strokovne in propagandne informacije in članki.
Cene objave tovrstnega materiala se bodo določale spo-
razumno.
ADVERTIZING RATES
Cover page (for all issues of 1982)
2nd page 1300 %
3rd page 1000 %
Inside pages (for alt issues of 1982)
1/1 page 790 $
1/2 page 520 $
Inside pages (Individual issues)
1/1 page 260 $
1/2 page 200 %
Rates for classified advertlzing:
each ad • 66 $
In addition to advertisment, we welcome short business
or product news, notes and artlcles. The related cliarges
are negotiable.