- Email: [email protected]
An algorithm for the generalized eigenproblem
Comput. & Elect. Engng Vol. 5, pp. 311-320 © PergamonPress I.td., 1978. Printed in Great Britain
0045-7906/7811201--0311/$02.00/0
AN A L G O R I T H M FOR T H E G E N E R A L I Z E D
EIGENPROBLEMt
Louis F. GODBOUT, JR. Department of Electrical Engineering, University of Hartford, 200 Bloomfield Avenue, West Hartford, CT 06117, U.S.A. and DAVID JORDAN Department of Electrical Engineering & Computer Science, Box U-157, University of Connecticut, Storrs, CT 06268, U.S.A.
(Received 22 May 1978; received[or publication 19 September 1978) Abstract--A new method is presented that determines the eigenvalues and corresponding eigenvectors for
the most general form of the eigenproblem. An example illustrating the technique is included. The appendix contains a listing of the Fortran program developed for use on a Data General Eclipse/S200 computing
system. BACKGROUND
The generalized eigenvalue problem is concerned with the nontrivial solutions A and x of Ax = ABx
(1)
where A and B are n × n constant matrices. The complex values A are called eigenvalues while the complex vectors x are referred to as eigenvectors. When B has an inverse (1) reduces to the standard eigenvalue problem
B-JAx = Ax
(2)
which can be solved using Francis' QR algorithm[l] together with Sridhar and Jordan's method [2] to find the n eigenvalues and eigenvectors respectively. Whenever B does not have an inverse the reduction (2) cannot be accomplished, However, (1) can be rewritten as (AB-A)x = 0.
(3)
[AB-AI = 0
(4)
The A's are given by
where I'l represents the determinant of a matrix. Expansion of (4) yields
(7) (~) (n~-,) Inl A~ + ~ IM,"'IA~-' + ~ IMi(2)IA"-2+ " " + ~ IM~¢"-"IA+1- AI : 0 i=1
i=I
i=1
(5)
where (:) is the notation for a binomial coefficient. M,.°) is the ith possible matrix that can be formed by replacing j rows of B by the corresponding j rows of - A. Notice that AB-A is a regular pencil of matrices[3]. Also note from (5) that there are less than n solutions A to (4) when B is singular. Thus for this case one must also identify the number of solutions, r~ < n, along with the r~ solutions for A and x. tThis work was supported in part by NSF Grants GK-37841 and ENG-75-17962. 311
312
[.. F. GODBOUT, JR. and D, JORDAN
Previous work has been done with regard to the generalized eigenproblem. Moler and Stewart[4] have introduced the QZ algorithm (modified QR algorithm) to solve the problem. Peters and Wilkinson[5] also have a method. The following discussion presents an alternative approach to the problem using the regular pencil properties of hB-A. ALGORITHM Consider the generalized eigenproblem in the form of eqn (3) where B is singular. Since AB-A is a regular pencil of matrices it can be transformed by appropriate, square, nonsingular matrices P and Q to a strictly equivalent quasi-diagonal form where: (i) P(AB-A)Q = diag {AB~- A~, AN-I} (ii) The diagonal block A N - I is (n - fi) x (n - fi) where I is the identity matrix and N is a quasi-diagonal nilpotent matrix of the form 0.
1
"1 "0 0
1 (6)
N= '1 0 0.
1
l 0 (iii) the diagonal block AB~-A~ is ti × fi (r~ being the number of solutions to the problem) where fi = n - £ # ,
(7)
I=1
with the gi's being the v block sizes of N, and (iv) B~ is nonsingular. The proof of these results and an algorithm for generating the canonical pencil along with the P and Q transformations is given in Godbout and Jordan [6]. The application of the canonical form technique to (3) yields P(AB - A)Qz = 0
(8a)
or
{:;J: {:]
['LA where z = Q Ix.
(9)
partitioned vectors zl and z2 of z are ri and (n - ti) dimensional respectively. The partition XBt - At is ti x ri while the A N - I p a r t i t i o n is ( n - ri) X (n - ti). The upper partition of (8b) yields the equation
The
(ABI - A t ) z l = 0.
(10)
313
An algorithm for the generalized eigenproblem
However since B1-t exists one gets (ll)
(AI - Bt-lAl)Zl = 0
which is the standard eigenvalue problem. Francis' OR method[l] can be employed to find the ri eigenvalues A while Sridhar and Jordan's algorithm[2] will generate t h e corresponding ti eigenvectors z~. The lower partition of (8b) gives (12)
(AN - / ) z 2 = 0.
Because of the form of N the [AN- I[ is unity (plus or minus) and no new ,~ values will be found, Also note that ( A N - / ) - I exists forcing z2 = 0.
(13)
With the knowledge of z~ and z2 for each A as found by (11) and (13) the desired eigenvectors x are given by x=Q
[z,] 0 •
(14)
EXAMPLE
Consider the generalized eigenvalue problem (AB - A)x
=0
where
B=
ii0 0, 1 Ii°°°°l i I'212l°°1 0 1 0 0 0 0 -l 0
0 1 0 -1
0 -1 0 2
00 0
and
0 0 -1 0
A=
.
-l
0 0 1 0 0 • 0 0 0 0 0 0
Using Godbout and Jordan's algorithm[6] P and Q transformations can be found to bring ;tB - A to canonical form. They are
I~ 1/2 0 1/2 P =
0 112 0 0
0 0 l 0
0 0 0 - 112 0 0 0 0 -l
and
-1/2
0 1/2 112 0
Q=
o
0 1/2 112 0
0 1 1 0
0 1 0 0 • 0 0 l 0
The result of applying these transformations to AB - A is P ( A B - A) Qz = 0
0 0 01
o o t A I--101/2 -11/2 0 blo 0 0 0 0 olo
o[
VI[2
[. o
ol-
O_l
-1/2
o
0 0 0
zl
o o
1 0
z2
o
0 0
0 1
0 0
=
1
The upper partition of the above yields (AI - B t - I A t ) z l = 0
or
[ -1/2
~ + 1/2'
- 1/2
~+l/2jz~=o.
[:1
314
L . F . (JODBOUT.JR. and D. JORDAN
The eigenvalues a and corresponding eigenvectors z~ are
A=O
with
z~-[{
with
z,=
1
and
a=-I
['1 -I
"
The desired eigenvector x for a = 0 is given by (14) as
X=
]/2 0 I/2 1/2 0
I~
I 0 0 0 1 0 0 I 0 0 0 I 0
-1/2
0 I/2 - 1/2
0
In the same manner the eigenvector x for a =
x =
=
0 ~) .
'1
0
1 is
I
.
o
0
CONCI.USIf)NS
An algorithm has been outlined to find the eigenvalue and eigenvectors of the most general
eigenproblem. The properties of regular matrix pencils were used to transform the problem to a form in which proven methods (QR method, Sridhar and Jordan method) could be applied to obtain the solution. The appeal of the technique lies in its simplicity and use of existing processes. An example has been given to illustrate the computational procedure. REFERENCES 1. J. A. F. Francis, The QR transformation, Part 1 and I1. Comput. J. No. 4. 265-271 and 332-345 (i962L 2. B. Sridhar and D. Jordan, An algorithm for calculation of the Jordan canonical form of a matrix. Comput. Elect. Engng 1,237-254 (1973). 3. F. R. Gantmacher, The Theory o[ Matrices, Vol. 2, Chelsea, New York (1959). 4. C. B. Moler and G. W. Stewart, An algorithm forthe generalized matrix ¢igenvalue problem. Rep. CNA32, Center for Numerical Analysis, Univ. of Texas at Austin. 5. G. Peters and J. H. Wilkinson, Ax = )tBx and the generalized eigenproblem. SIAM ,L Numer Anal. 7, 479-492 (t970). 6. L. F. Godbout and D, Jordan, On state equation descriptions of linear differential systems. Trans. ASME Dynamic Systems Meas. and Contr. 97, 333-344 (1975). APPENDIX The following Fortran program successfully implements the solution to the generalized eigenproblem discussed in this paper. The program has been written to run on a Data General Ecfipse/S200 computer. Only slight modifications are required to have it run on other systems. The program uses as input the A and B matrices. If it is determined that B is nonsingutar the output generated is simply the B-'A product. However, if B is singular the canonical form indicated in the paper is generated. For this case the B, ~A. product together with the Q transformation is the outpuL In either case the user should implement the resultant output to solve the remaining standard eigenproblem. COMPILER DOUBLE PRECISION •DIMENSION A ( | 0 o I 0 ) , B ( I 0 . I 0 % . P t I O , 1 0 I p B I N V A ( I O , 1 0 )
DIMENSION
G(IO,IO)#U(IO,IO),T|f10,IO)
DIMENBION O U M M Y t I 0 , 1 0 ) , S C R A T C H f I O . 1 0 ) DIMENSION X ( 1 0 ) ° X I ( I 0 ) , X 2 ( I ~ } o Y f I 0 ) , Y ! ( 1 0 ) DIMENSION I N C O L t I 0 } , I C L ~ ( I O ] LOGICAL ONEMAT,NOGO INTEGER ~ANK,COLM,ETA,ETANFW C • C
OPENS CARD READER AND LINF PRINTER AS INPUT ANn OUTPUT
An algorithmfor the generalized eigenproblem C C c
DEVICES R E S P E C T I V E L Y . THIS lS A NECESSITY ON A DATA GENERAL E C L I P S F / S ~ 0 0 SYSTE M.
WHEN RUNNING SATCH
CALL O P E N ( | , ' T M P ~ , t e I E R ) CALL OPEN(2,'SYSOUTW,3,1ER) 19 w R I T E ( 2 , t ) ] FORMAT(/,SX,'IHIS ROUTINE DFTERMINES THE EIGENVALUES,LAMRDAo IAND THE EIGENVECTORS,X, THAT SATISFY AtX=LAM~DAtR*X.') C C C C C
READ IN THE OATA. N=SYSTE~ ORDER. EPS=ALLDWABLE ERROR. IF IEND>O ANOTMER SYSTEM DATA SET FOLLOWS T H I S SET. READ(I,2)N,IENO,EPS FOR~AT(2110,D?O.IO) REAO(I,~)((8(I,J),J=I,N),I=I,N] HEAD(I,3]((A(I,J],J=t,N),I=toN) FORMAT(~020.10)
C C C
PPINT THE DATA.
8 q 5
6 7 C C C C
~RITE(2,a)N FORMATf/,~X,'IHE SYSTEM OWnFR TS ° , 1 2 , ' . ' ) ~ITE(2,8] FURMATf/,SX,'THE B ~ATRIX T S ' ) DO G 1 = 1 , ~ , ~RIT~(2,7)(~fI,J),J=I,~) ^wIIE(2,5) FO~HAT(I,~X,'IHF A ~ATRIx I S ' ) DO ~ I = 1 , ~ ~ITE[2,/i(A(I,J),JmloN} FbR~AT(bD25.15) PE~FOW~ G~USS ELIMIwATIhN ON ~. 5Av@ ~ T~ OU~'~Y I~ CASE R INVFRSE EXISTS. I F ( ~ . F : } . I I G O T{I 10 O~F~AI~.FALSE. CALL GAUSS(A,~,P,]HCOL,N,UNFMAT,FPS]
C C C C C C C C
NIiTF:~ A~O A A ~ TRANSF(IR~FD MATRICES OF THF ORIGINALS. Pl AND q l a ~ IRE T R A N S F O P ~ A T I O N S . ~1 I~F(,q~AIIu.~ II~ ]NCOL. Pl TN P. A M A ~ I~ A. H~AR IN ~. F]~,O ] ~
RA~K OF H ANi) OUTPUT
IT.
II=i~-I+1 I F { ~ a ~ S E P t I I , I I ) ) . G I . ~ P S ) G P TO 12 11 CONTINUE 12 WRIIE(2,1]]RANK I ) FORMATt/,SX,'THE RANK OF B TS ' , I 2 , ' . ' ) C C C
LOOK FOR AN ERROR. IF(RANK.NE.O]GO TO tA ~RITE(2,15) 15 FORMAT(/,SX,'THE 8 MATRIX GO TO 5 2
C C C C
IS
NULL.')
IF 8 INVERSE EXISTS THEN DETFRMINF I T . OUTPUT THE 8 INVERSE*A PRODUCT. t~
IF(RANK.NE.N)GO TO 16 wRITE(2,SB) 58 FORaAT(/,$X,'B INVERSE ~ X I S T S . ' ) CALL INVERSE(B,P,INCUL,SCRATCH°N) CALL ~MULI(SCRATCH,DUMMYoMTNVA°N,N,N) ~R]TE(2,59) ~g FOWMATf/,SX~'IHE 8 INVERSE*A MATRIX IS ~) DO bO I = I , N bO ~ R I I E ( 2 , 7 ] ( ~ I N V A ( I , J ) , J = I , N ) GO TO ~2 C C C C E C C C C C
~ OOES Nor HAVE AN INVERSF. T H E R E F O R E GO T H R O U G H THF R E D U C T I O N ~1
PROCESS.
UPPER LEFT PARTITION OF BRAR.
F O ~ M 0 2 ANO P~, 92 IN (~. P2 INVERSE IN P. ETA=MAXIMUM DIAGONAL HLOCK S I I F
COLM=I O,~E~AI=.TwUE. ETA=! O~ ~ d ]=IRAN~,N
OF THE NILPOTFNT m A T R I x .
315
L. F. GODBOUT, JR. and D. JORDAN
316 E T~'~E ,~:O C C C
1;~ TI-~', I : , , f
g
COLu~'r~
IIF
,,2.
]F ( I - j ) 2 ~ , 85,,~,4 2g
~;[J,COL~,~] =1. ohO :{ { J ) : 1 . (h")O G(! 1 0 # ~
#~1
,~, ( j , C OL ~'l ) = 0 . b4)(~ x[J]:O.OnO CI}~ 1 l ~ . u ;
8~
P6
;(# ( J 1 : - ' ~ (
CALL
J,
l )
~ A C S u H ( ~ , X~, x l ,RANK)
DI, ~ 7
J=I,~4~'JK ,.(j, C,3L~.') =X 1 [ J ) X{J)=,l[j]
27 C C
r.f/~,,,~
vF['Ta)R
a
UF
PP" TNVF~SF.
C t3~'
CAt L k A I v [ " C ( A , X , Y,N) r [,(;n=.l- aLSE . I)tJ c~'4 J : l p~i
P(j,COLM]=Y[J) IF ( J - R A N K ] 3 0 , $ O , 3 1 Y]{J}--Y(J)
30
GO
31
TO
29
IF { {OABS{Y ( J ] ) . G T . E P S ) . A N D . [ . N O T . COtx, T INUE
2Q
C C C
C~ECK FOR A SINGULAR PENCIl_.
NOGOI )NOGO=.TRUE.
USE U FOR STOPAGF.
O0 3~ J = I , N DU 3l W : I , N ! F [ K - C O L M ] 38# 3 8 , 3 Q U(j,K]=P{j,K)
38
GO TO 3 7 U(j,K)=O.ODO COP,'T l RuE CONTINUE CALL GAUSS [ DUMMy ~ U, SCRATCH, ]'CLMo N, ONEMAT, EPS) I F ( O A ~ S [ U [ C O L m , C O L N ' I } . G T . F P S ) C , O TO ~0 K.#ITE(2,~I) "l F O R M A T ( I , S X , ' L A ~ B D A * B - A I S A SINGULAR P E N C I L . ' , / , 1SXo'IHERFFO~'E THERE AreE NO SOLUTIONS TO T H I S P P O R L E F . * ) GO TO .~2
39 37 ].~}
C C
SE' FAP TF~E PENCIL fS NOT STNG!JLAP. vECTORS L)E P2 AND 0 ~ .
C
CONTINUE
FINDING
C O L U MN
C COL ~:COL ~", 1 F I ANF ¢,=F f AfgE"~+ 1 IF(,'~OGO)GO TO 21 CALL ~ / ~ C S u H [ R , Y I , X I o RAN~ ) Ou ~ J=I,N
qO
,~5
G( J,CUL~')=O.UUO
~a
x(Jj:o.nno r:-L. [ r , s ' , ~(J,EOL,~')mxI
[d)
X ( J 1::'(1 [,J} CU~' T I ,U~" Gr TO d 'a ]F(EIA,\,FW . ( ; [ .
3~
2.1 C C
]r,,)ILAIF
IHF
ETA)FTA:FTANF~"
NU~"FJlr~ Ok S('ILL,TI[JNS THAT
EXISI.
C
CuL-~:C CI[. '~'" 1 •:'
I IJ
~ ;
C
C C
~=ILL ~1
I),'
(:UT
u~
K:I,>
: i~
,a q
.~:
,J I =hi*
j j:
T~F P 2 ] N V
C~'
I #l,I
} -
)
,I-CI'~L a+J ] I F [ J l -COLM)
~. ~9
~D
~18, ~le~, q 7
IF ( K - I ) ~ 8 , ~ l g , a 8 U(I,J1)~Q(I,,JI) Q(I,JJ}=Q(I,JI} GO TO ~ 5
T~A~;SFFII4~AT~ONS.
An algorithmfor the generalizedeigen~oblem ~8
U{IwJl]wPfIoJl} PtloJJ)sP(l,Jl) GO TO ~5
a7 ~5 Q~
U(IoJI)=0.000 CONTINUE CONTINUE CALL GAUSS(OUMMYeU°Tt°ICLM°N,ONEMAT,EP$)
CALL GAU88(DUMMYwTI.SCRATCH,ICLMeNeONEMAT~EP$) CALL INVERSE{TI,SCRATCM,ICLM,OUMMYeN) DO 50
I=I,N
ICM=COLM+I 00 %1 J = I C M , N JJ= J- COLM
IF(K-I)S2wS~,S2 $3
52 51 50 ~3 C C C C C C
~(I,JJ]=SCRATCH(I.J) GO TO 51 P(I,JJ)=SCRATCMII,J] CONTINUE CONTINUE CONTINUE
OETEwMI~E P~ FROM ITS INVFRSF. T l USED TO STORE P 2 . HwIGGLE IN H. A~IGGLE IN A.
CALL GAUSSCOUM~YeP,SCRATCM, ICLM,N°ONEMAT,EPS) CALL INVERSE(R,SCRATCMoICLM,TI,N) C C C
DETERMINE H~IGGLE AND AWTGGLE. CALL CALL CALL CALL
C C C
MVULTCTI,R,DUMMY,NoN,N) M~ULT(DUMMY,Q°d°NoN°N) ~UULTCTI°A,OUMMY,N,~,N} MMULT(OU~Yo~,A,N°N°N)
FOR~' Q=UI*02
hl
C C C C C C
O0 bl I = l , ~ ICLM(]}=I CALL ~SOwIIO,N,Ir~COL,ICL~)
EXTRACT T~F ~At~Ix PARTITIONS. ~ l ]~] ~. A[ IN A. ~ IN H . - i ~ In T I ° f i l I:~VERSF*AI IN wII~vA. N TN P. TME~. SFI l ~ h LO~E~ LEFT PAWTTTTONS TO ZERO. ~U=N-~EIG DO 71 I = l , , ~ U TI:MFIG+i O0 ~ J=T,MFIG U(I,J)=H(II,J) ~(IIoJ)=O.ODU IT(]°J)=-A(II,J] 6~ A(II,,I)=O.OPn ~U b5 J=t°~U JJ=~EIG+J b5 RtI,JI=H(IIoJJ) 7] CO~TINu~
C C C
P~I~'T THE FINAL TkANSFOR~FD MATRICES, nWIIF(2,b~) 5~ FN~AI(IoSXo'T~E T~ANSFO~F~ ~ ~ATRIX I S ' ) O0 b~ I = I , ' ~ 55 ~ I I ~ ( ~ , / J [ R ( | , J J , j = I ° ~ ) %0 F O ~ V A T ( / , ~ X , ' T ~ F T~A~SFORMF~ A ~ATRIx IS ° ) DO 57 ]=1,~', ~7 ~ | T ~ ( ~ , ? ) C A C I , J } , J : ] ° ~ ]
C C C
FOO~ ~1
I'~VFkSE*A].
CALL GAUS$(DUM~YoH°SCkAtCM°TNC[}L,MEIG#ONEMAToFPS} CALL I'~VFkSE[~,SCWAICM,T.~COI°CiIMMY°MEIG) CALL u~ULT(DIJ~Y,A,MI~A°~FT~°~FTG,~FIG}
C C C C
Ft.~v APHAI:DLIMY~=M~*RINVA-A~. AGHAI l~'J SCRATCH. CALL MMUL T (LJe~[~JVAeDt~Y,NIJ,~FTG°~EIG) CALL ~AOO(OU~MY, TI ° SCRATEN.~II°~F I~)
~UILO Q~ ] ~
~
LN~ER LEFT P A R T I I I N N OF THE TRANSFORMATION Q3. 1 1 . LURER LEFT P A R T I T I O N OF 03 IN U.
CALL E Q u A T E ( S C W A T C H p u , n U , M F I G } I F ( ~ T A . E Q . 1}GI) TO b7
317
L. F. GODBOUT. JR. and D. JORDAN
318
CALL r ~ U L I ( P , S E R A I C M , D u ~ Y , N U , N U , ~ E I G ) CALL ~ U L T I D u ~ Y , ~ I N V A , ~ C R A T ~ , N U , ~ E I G ° ~ E I G ) CALL r~AD0(LI~SC~AIC~,u°NU,~FT~I EIA=EIA-I GO TO ~ kOIL n C~. ~7 DO 70 I : I , ~ DO l~ J = I , N IF 1 - j ) 7 ~ , 7 ~ , 7 ~ 7Q T! I , j ) : I . 0 D 0 GO T0 72 73 IF ] - ~ E I G ) 7 5 , 7 5 , 7 ~ I1 I , J ) = O . O 0 0 75 GO TO 7 2 ?b IF(J-MEIG)77,/7,75 77 II=I-~EIG Tl(l,J)=U(II,J} 72 CONTINUE 7 0 CONIINUF DETER~II~E Q = I Q 1 * Q 2 % * Q J . G IN DUMMY. CALL ~ U L T ( Q , T t , D U ~ Y , N , N , N ) PPINT ~I
INVERSE*A|
AND THF Q TRANSFORMATION.
202
FO~AI(/,bx,'TME QI I N V F ~ S F * ~ M A T P I X IS') DO IU0 I = I , M E I G I00 ~ Q ] T E [ 2 , 7 I ( h l N V A C I , J I , J = I , M F T G ) ~R]TE(@,7~) 78 F O k M A T ( / , S x , ' T M F Q INANSFO~ATION IS') UU 7 9 I = I , N 79 ~ I T E f ~ , 7 ] ( O U ~ M Y ( I , J I , J = I , N ]
C C C
00NE) ],2 I F ( I E N D .LE. ~WITE(2,6~I b~ FO@MAT(IH)) GO T0 lq ENO
0)SrUP
GAUSS ~LI~I~,AIIO~,J
~OUTINF.
CG~PILPW F)UUWLF PP~CIHTLIi~ SUP~flUTI~;E GAUSS(~MATmMMAT,PN~AT,TC,KDI~,ONF~AT,FPS} l)I~,k~Sl*) 'V A ~ A T ( I 0 , 1 U ) , ~ A I f I 0 , 1 0 ] , P M A I ( 1 0 , I 0 I , I R r I 0 ) , I C ( 1 0 ) , T p ( I ~ ) LOGICAL Oh E~AT I~TTTALIZE ~O~ ANI) COLUMN TNDFX ARPAyS. T ~ , I r l a L I Z k NUa TRA~FOW~ATInN MATPTx TO TMF I ~ N T I T Y .
O0 10 I = ] , K O I ~ i~(II=T IC[I]=I IPCI]=T OC ~0 J = l , ~ ) i ~ IF(I-J}@~,~Om~b ~0 P'~AT(I,J)=I.0DO 25 P~A I ( I , J ] =0. OI)U ~0 C ~ T I ,~HE 10 CONTI!~U~ Pk~F(~ OO
~n
(~URS ELI~INATION USING FULL PIVOTING.
l=l,~I
~
I I ~ I F I A L I Z E THE PlvDT
FLFMFNT.
PI~oT:~AT(I~{II,TCfl)) ICOL~=I FIND
T~E PIVOT
ELFmENI.
DO 50 J = I , K D I ~ D0 O0 K=I,KDIM I F ( D A R $ [ S M A T ( I ~ ( K % , I £ ( J ) ) ) . L E . D A B S ( P I V O T I ] G OTO bO PIVOT=B~AT(]~(~),TC{J)}
bO 50
ICOL~=J CONTINUE CONTINUE CMFCK
THF
~AGNITUDE
OF
1F(DA~S[PlVOT].LT.EPS)GO C
TMF
PIVOT
TO 7 0
ELEMENT.
An algorithm for the generalized eigenproblem CHANGF RO~S
IF NEEDED.
IF(IR(IRO~).EQ.IR(1))GO TO 7% JRO~=IR(1) I~(I)=IR(IRO~) IR(IRU^)=JRO~ CHA~GE COLUMNS IF NEEDED.
7S
]F(IC(ICOLM).EQ.IC(I))GO T0 80 JCOLM=IC(II IC(1):IC(ICOLM) IC(ICOLMI=JCOLM ELIMINATION FOLLOWS.
~U
IF(I.bE.~DI,a)~O rn 70 II=I*l i){' ~u , ] = l l , K I ) l ~ OtIg#~lr,*~ A~ EL~M~I OF T~F ~O~ T~ANSFORMATION MATRIX. RtL~M:~,~AI(IW(J)oIC(1))/PIVGI ALOE, R ~OLIIRLE OF IHP PIVOT ~(;~ TO ANOTMER ROt.
~S ¢5
~AI(r~(j],lC(tl)=O.nnO
O0 qb ~ = ] I p g l ) I M
MvATCIR(j),IC(^))=HMAT(I~(J),ICfKI)-RELE~*MMAT(IR(TI,IC(K)) O0 1.)0 K=l*gt) J "i IFt.i~OT.ONE~AIJAMAT(IRfJ)tIC(K))=A~AT(IR(JI,ICCK)) I -RELEM*A~AT(IR(T),IC(K)) tO0 R~AT(I~(J},K) =P~AI(IR(j)'W)'pELEM*PMAT(IR(T)*K) @0 CONTINUE ~0 CO~[INU~ S(JRT TME MATRICES.
70 CALL ~SORT[~MAT*KDIM*IR~IC) CALL uSOWT{P~AI*KOI~ ' | R ' I P ) IF(.NOI.h~EPAI}CALL ~SORT{AMAT,KDIM~IR,IC) RETURN EW ,D shPTT*~'G wOUII~E FO~ ~ATRTCFS. CO~PILf~ hOi,PLE P~t~ISIUN SUhRUU I II~*~ ~'SURI{ S~AToJi)IM, JR, JC) UI N~r,~ I 0,,, S~,AT( I O#IO),DUMMY(IO°I 0 ) ,JR(TO) ,JC(TO) O0 110 I=I,JC)/~ 0{~ ~IO J=loo~)l ~ 110 I ) u ~ Y ( I , J ) = S ~ A I ( J R { I ) ° J C ( J I ) CALL ~:~!~Al@ ( OUt'Y, S',AT, JDI~° JOT~)
C C C'
SUBROUTIN(
FINDS
THE INVERSE OF A M~TRIX.
COMPILER DOUBLE PRECISION SUBROUTINE INVERSE(UMtT,PHAT,ICOL*TINV,LDIM) DIMENSION U H A T ( I 0 , 1 0 I , P H A T I 1 0 , 1 0 I , T I N V ( 1 0 , I O I , I C O L ( 1 0 I * IPCOL(10),TCOL(|0) EXTRACT a COLUMN FROM PM~T TO FINO A COLUMN OF TMF INVERSe.
Sl0
DO 800 J = I , L D I M DO 8 1 0 IffiI,LDIM PCOL(I)=PMAT(I,J) FIND A COLUMN OF THE INVERSF.
CALL BACSUflfUMATePCOL~TCDL,LDIM)
FORM THE INVERSE FROM ITS COLUMN VECTORS. DO 820 K=|~LDIM 820 TINV(ICDL(K)*J)=TCOL(K) ~00 CONTINUE RETURN END eAC~ SUBSTITUTION ROUTINE. COMPILER DOUBLE PRECISION SUBROUTINE ~ACSUB(URMAIoCOLU~TSOL,JLENGI DIMENSIO~ u R M A T t l 0 , I O ) , C O L U M ( 1 0 ) ~ T S O L ( 1 0 ) FINO IME LAST ELEMENT OF THE SOLUTION VECTOR FIRST,
TSOL(JLENG)=COLUM(JLENG)IUPMATIJLENGeJLENG} C
319
320
L . F . GODBOUT,JR. and D. JORDAN FI~O
ThE
OTMEW
ELE~ENIS
nF
TMF
SOLUTION
VFCTO~.
JJ=JL~G'I Ob ~ 0 J : 1 , J J JA:JLF~G-J 5TOT:0.000 Ja]=ja+1 Dn ~ 0 T:JAI,JLENG KW:JAIejLENG-I B~O SILIT:STOT+UPMAI[j~,K~)*TSOL(KK] 830 T $ O L ( J A ) = [ C O L U M [ J A ) - S T O T ) / I I P ~ A T ( J A , J A ) ~TUR~
~OLTIPLIES
920 910 g00
x MM
MATRIX
~Y
AN MR X PP
COMPILEW OUUPLE PRECISION ~U~W~uT]N~ ~ M U L T C X M A T , Y M A T , I M A T , N N , M M , P P ) DI~E~,SION x M A I ( | 0 , 1 0 ) , Y M A T ( 1 0 , 1 0 } , Z M A T ( ~ 0 , 1 0 ) I N T E G E ~ PP OU 900 I : ~ , ' ~ N OO 910 J : I , P P SU~IqT:0.000 O0 qPO K : I , ~ M 5U~'TnT:Su~TOTeXMAT[T,K)*Y~AT(K,j) Z~AI(I,J):SJMTOT CON~]~UF ~ETuw~ E~D
~ULTIFLIE5
2~I 220
AN N N
A '~AIQIX MY A VFCTDR.
COMPILFP DOUBLE PRECISION SUbroutiNE MATvEC{TMAT,INVFC,OUTVFC,DIMET] OIMES;SIO~' I M A I [ 1 0 , 1 0 ) , O U I V F C ( } 0 ) I~TEGEQ OlmEI ~EAL I , , ' v F C ( 1 0 ] O0 2 ~ 0 I : I , O I M E I SUM~Y=0.0D0 ()U ~ I J = [ , i ) I~ET ~UM~=bU~MY+T~AT(I,J}*INVFC(J) O U I V E C ( I J = ~ U ~MY wETuP~ END
~OOS
T~C~ ~ a l " ~ I C ~ S °
COMPILER DOUBLE. P R E C I S I O N S U ~ P U U I INE ~ A O D (MAT I, MA T2, M A T S U ~ , NR, NC) Q E A L ~AT I ( I O, I O } , M A T 2 [ tO° |t)I , M A T S U M { I0° IO) O0 7 I : I , ~ , R DO 7 J : L , N C 7 MATSUM{ I , J ] :~.A I I { I , J ) ÷ M A T 2 ( I , . J } RETUQN END C C C"
EI~L,#TES ONE '.,ATRIX
TO ANOTMFR.
CO,PILaW rJOUBLE PHECISION SU~POU I I:H" ~_UUATE { IN"~A T , (%UT~A T, N 1 , NP) OI ~E {,,S I ~,~ O O T ~ A T ( 1 0 , 1 0 } RLAL I ~ , V A T ( I Q , ] 0 ) i)O 17 I : 1 , ~ ' ~ 1 OU 17 j--],r~_P 17 otJT~a I { I , J ) = I N~'+AI { ] , J J WE f OPt,
MATRIX.
0045-7906/7811201--0311/$02.00/0
AN A L G O R I T H M FOR T H E G E N E R A L I Z E D
EIGENPROBLEMt
Louis F. GODBOUT, JR. Department of Electrical Engineering, University of Hartford, 200 Bloomfield Avenue, West Hartford, CT 06117, U.S.A. and DAVID JORDAN Department of Electrical Engineering & Computer Science, Box U-157, University of Connecticut, Storrs, CT 06268, U.S.A.
(Received 22 May 1978; received[or publication 19 September 1978) Abstract--A new method is presented that determines the eigenvalues and corresponding eigenvectors for
the most general form of the eigenproblem. An example illustrating the technique is included. The appendix contains a listing of the Fortran program developed for use on a Data General Eclipse/S200 computing
system. BACKGROUND
The generalized eigenvalue problem is concerned with the nontrivial solutions A and x of Ax = ABx
(1)
where A and B are n × n constant matrices. The complex values A are called eigenvalues while the complex vectors x are referred to as eigenvectors. When B has an inverse (1) reduces to the standard eigenvalue problem
B-JAx = Ax
(2)
which can be solved using Francis' QR algorithm[l] together with Sridhar and Jordan's method [2] to find the n eigenvalues and eigenvectors respectively. Whenever B does not have an inverse the reduction (2) cannot be accomplished, However, (1) can be rewritten as (AB-A)x = 0.
(3)
[AB-AI = 0
(4)
The A's are given by
where I'l represents the determinant of a matrix. Expansion of (4) yields
(7) (~) (n~-,) Inl A~ + ~ IM,"'IA~-' + ~ IMi(2)IA"-2+ " " + ~ IM~¢"-"IA+1- AI : 0 i=1
i=I
i=1
(5)
where (:) is the notation for a binomial coefficient. M,.°) is the ith possible matrix that can be formed by replacing j rows of B by the corresponding j rows of - A. Notice that AB-A is a regular pencil of matrices[3]. Also note from (5) that there are less than n solutions A to (4) when B is singular. Thus for this case one must also identify the number of solutions, r~ < n, along with the r~ solutions for A and x. tThis work was supported in part by NSF Grants GK-37841 and ENG-75-17962. 311
312
[.. F. GODBOUT, JR. and D, JORDAN
Previous work has been done with regard to the generalized eigenproblem. Moler and Stewart[4] have introduced the QZ algorithm (modified QR algorithm) to solve the problem. Peters and Wilkinson[5] also have a method. The following discussion presents an alternative approach to the problem using the regular pencil properties of hB-A. ALGORITHM Consider the generalized eigenproblem in the form of eqn (3) where B is singular. Since AB-A is a regular pencil of matrices it can be transformed by appropriate, square, nonsingular matrices P and Q to a strictly equivalent quasi-diagonal form where: (i) P(AB-A)Q = diag {AB~- A~, AN-I} (ii) The diagonal block A N - I is (n - fi) x (n - fi) where I is the identity matrix and N is a quasi-diagonal nilpotent matrix of the form 0.
1
"1 "0 0
1 (6)
N= '1 0 0.
1
l 0 (iii) the diagonal block AB~-A~ is ti × fi (r~ being the number of solutions to the problem) where fi = n - £ # ,
(7)
I=1
with the gi's being the v block sizes of N, and (iv) B~ is nonsingular. The proof of these results and an algorithm for generating the canonical pencil along with the P and Q transformations is given in Godbout and Jordan [6]. The application of the canonical form technique to (3) yields P(AB - A)Qz = 0
(8a)
or
{:;J: {:]
['LA where z = Q Ix.
(9)
partitioned vectors zl and z2 of z are ri and (n - ti) dimensional respectively. The partition XBt - At is ti x ri while the A N - I p a r t i t i o n is ( n - ri) X (n - ti). The upper partition of (8b) yields the equation
The
(ABI - A t ) z l = 0.
(10)
313
An algorithm for the generalized eigenproblem
However since B1-t exists one gets (ll)
(AI - Bt-lAl)Zl = 0
which is the standard eigenvalue problem. Francis' OR method[l] can be employed to find the ri eigenvalues A while Sridhar and Jordan's algorithm[2] will generate t h e corresponding ti eigenvectors z~. The lower partition of (8b) gives (12)
(AN - / ) z 2 = 0.
Because of the form of N the [AN- I[ is unity (plus or minus) and no new ,~ values will be found, Also note that ( A N - / ) - I exists forcing z2 = 0.
(13)
With the knowledge of z~ and z2 for each A as found by (11) and (13) the desired eigenvectors x are given by x=Q
[z,] 0 •
(14)
EXAMPLE
Consider the generalized eigenvalue problem (AB - A)x
=0
where
B=
ii0 0, 1 Ii°°°°l i I'212l°°1 0 1 0 0 0 0 -l 0
0 1 0 -1
0 -1 0 2
00 0
and
0 0 -1 0
A=
.
-l
0 0 1 0 0 • 0 0 0 0 0 0
Using Godbout and Jordan's algorithm[6] P and Q transformations can be found to bring ;tB - A to canonical form. They are
I~ 1/2 0 1/2 P =
0 112 0 0
0 0 l 0
0 0 0 - 112 0 0 0 0 -l
and
-1/2
0 1/2 112 0
Q=
o
0 1/2 112 0
0 1 1 0
0 1 0 0 • 0 0 l 0
The result of applying these transformations to AB - A is P ( A B - A) Qz = 0
0 0 01
o o t A I--101/2 -11/2 0 blo 0 0 0 0 olo
o[
VI[2
[. o
ol-
O_l
-1/2
o
0 0 0
zl
o o
1 0
z2
o
0 0
0 1
0 0
=
1
The upper partition of the above yields (AI - B t - I A t ) z l = 0
or
[ -1/2
~ + 1/2'
- 1/2
~+l/2jz~=o.
[:1
314
L . F . (JODBOUT.JR. and D. JORDAN
The eigenvalues a and corresponding eigenvectors z~ are
A=O
with
z~-[{
with
z,=
1
and
a=-I
['1 -I
"
The desired eigenvector x for a = 0 is given by (14) as
X=
]/2 0 I/2 1/2 0
I~
I 0 0 0 1 0 0 I 0 0 0 I 0
-1/2
0 I/2 - 1/2
0
In the same manner the eigenvector x for a =
x =
=
0 ~) .
'1
0
1 is
I
.
o
0
CONCI.USIf)NS
An algorithm has been outlined to find the eigenvalue and eigenvectors of the most general
eigenproblem. The properties of regular matrix pencils were used to transform the problem to a form in which proven methods (QR method, Sridhar and Jordan method) could be applied to obtain the solution. The appeal of the technique lies in its simplicity and use of existing processes. An example has been given to illustrate the computational procedure. REFERENCES 1. J. A. F. Francis, The QR transformation, Part 1 and I1. Comput. J. No. 4. 265-271 and 332-345 (i962L 2. B. Sridhar and D. Jordan, An algorithm for calculation of the Jordan canonical form of a matrix. Comput. Elect. Engng 1,237-254 (1973). 3. F. R. Gantmacher, The Theory o[ Matrices, Vol. 2, Chelsea, New York (1959). 4. C. B. Moler and G. W. Stewart, An algorithm forthe generalized matrix ¢igenvalue problem. Rep. CNA32, Center for Numerical Analysis, Univ. of Texas at Austin. 5. G. Peters and J. H. Wilkinson, Ax = )tBx and the generalized eigenproblem. SIAM ,L Numer Anal. 7, 479-492 (t970). 6. L. F. Godbout and D, Jordan, On state equation descriptions of linear differential systems. Trans. ASME Dynamic Systems Meas. and Contr. 97, 333-344 (1975). APPENDIX The following Fortran program successfully implements the solution to the generalized eigenproblem discussed in this paper. The program has been written to run on a Data General Ecfipse/S200 computer. Only slight modifications are required to have it run on other systems. The program uses as input the A and B matrices. If it is determined that B is nonsingutar the output generated is simply the B-'A product. However, if B is singular the canonical form indicated in the paper is generated. For this case the B, ~A. product together with the Q transformation is the outpuL In either case the user should implement the resultant output to solve the remaining standard eigenproblem. COMPILER DOUBLE PRECISION •DIMENSION A ( | 0 o I 0 ) , B ( I 0 . I 0 % . P t I O , 1 0 I p B I N V A ( I O , 1 0 )
DIMENSION
G(IO,IO)#U(IO,IO),T|f10,IO)
DIMENBION O U M M Y t I 0 , 1 0 ) , S C R A T C H f I O . 1 0 ) DIMENSION X ( 1 0 ) ° X I ( I 0 ) , X 2 ( I ~ } o Y f I 0 ) , Y ! ( 1 0 ) DIMENSION I N C O L t I 0 } , I C L ~ ( I O ] LOGICAL ONEMAT,NOGO INTEGER ~ANK,COLM,ETA,ETANFW C • C
OPENS CARD READER AND LINF PRINTER AS INPUT ANn OUTPUT
An algorithmfor the generalized eigenproblem C C c
DEVICES R E S P E C T I V E L Y . THIS lS A NECESSITY ON A DATA GENERAL E C L I P S F / S ~ 0 0 SYSTE M.
WHEN RUNNING SATCH
CALL O P E N ( | , ' T M P ~ , t e I E R ) CALL OPEN(2,'SYSOUTW,3,1ER) 19 w R I T E ( 2 , t ) ] FORMAT(/,SX,'IHIS ROUTINE DFTERMINES THE EIGENVALUES,LAMRDAo IAND THE EIGENVECTORS,X, THAT SATISFY AtX=LAM~DAtR*X.') C C C C C
READ IN THE OATA. N=SYSTE~ ORDER. EPS=ALLDWABLE ERROR. IF IEND>O ANOTMER SYSTEM DATA SET FOLLOWS T H I S SET. READ(I,2)N,IENO,EPS FOR~AT(2110,D?O.IO) REAO(I,~)((8(I,J),J=I,N),I=I,N] HEAD(I,3]((A(I,J],J=t,N),I=toN) FORMAT(~020.10)
C C C
PPINT THE DATA.
8 q 5
6 7 C C C C
~RITE(2,a)N FORMATf/,~X,'IHE SYSTEM OWnFR TS ° , 1 2 , ' . ' ) ~ITE(2,8] FURMATf/,SX,'THE B ~ATRIX T S ' ) DO G 1 = 1 , ~ , ~RIT~(2,7)(~fI,J),J=I,~) ^wIIE(2,5) FO~HAT(I,~X,'IHF A ~ATRIx I S ' ) DO ~ I = 1 , ~ ~ITE[2,/i(A(I,J),JmloN} FbR~AT(bD25.15) PE~FOW~ G~USS ELIMIwATIhN ON ~. 5Av@ ~ T~ OU~'~Y I~ CASE R INVFRSE EXISTS. I F ( ~ . F : } . I I G O T{I 10 O~F~AI~.FALSE. CALL GAUSS(A,~,P,]HCOL,N,UNFMAT,FPS]
C C C C C C C C
NIiTF:~ A~O A A ~ TRANSF(IR~FD MATRICES OF THF ORIGINALS. Pl AND q l a ~ IRE T R A N S F O P ~ A T I O N S . ~1 I~F(,q~AIIu.~ II~ ]NCOL. Pl TN P. A M A ~ I~ A. H~AR IN ~. F]~,O ] ~
RA~K OF H ANi) OUTPUT
IT.
II=i~-I+1 I F { ~ a ~ S E P t I I , I I ) ) . G I . ~ P S ) G P TO 12 11 CONTINUE 12 WRIIE(2,1]]RANK I ) FORMATt/,SX,'THE RANK OF B TS ' , I 2 , ' . ' ) C C C
LOOK FOR AN ERROR. IF(RANK.NE.O]GO TO tA ~RITE(2,15) 15 FORMAT(/,SX,'THE 8 MATRIX GO TO 5 2
C C C C
IS
NULL.')
IF 8 INVERSE EXISTS THEN DETFRMINF I T . OUTPUT THE 8 INVERSE*A PRODUCT. t~
IF(RANK.NE.N)GO TO 16 wRITE(2,SB) 58 FORaAT(/,$X,'B INVERSE ~ X I S T S . ' ) CALL INVERSE(B,P,INCUL,SCRATCH°N) CALL ~MULI(SCRATCH,DUMMYoMTNVA°N,N,N) ~R]TE(2,59) ~g FOWMATf/,SX~'IHE 8 INVERSE*A MATRIX IS ~) DO bO I = I , N bO ~ R I I E ( 2 , 7 ] ( ~ I N V A ( I , J ) , J = I , N ) GO TO ~2 C C C C E C C C C C
~ OOES Nor HAVE AN INVERSF. T H E R E F O R E GO T H R O U G H THF R E D U C T I O N ~1
PROCESS.
UPPER LEFT PARTITION OF BRAR.
F O ~ M 0 2 ANO P~, 92 IN (~. P2 INVERSE IN P. ETA=MAXIMUM DIAGONAL HLOCK S I I F
COLM=I O,~E~AI=.TwUE. ETA=! O~ ~ d ]=IRAN~,N
OF THE NILPOTFNT m A T R I x .
315
L. F. GODBOUT, JR. and D. JORDAN
316 E T~'~E ,~:O C C C
1;~ TI-~', I : , , f
g
COLu~'r~
IIF
,,2.
]F ( I - j ) 2 ~ , 85,,~,4 2g
~;[J,COL~,~] =1. ohO :{ { J ) : 1 . (h")O G(! 1 0 # ~
#~1
,~, ( j , C OL ~'l ) = 0 . b4)(~ x[J]:O.OnO CI}~ 1 l ~ . u ;
8~
P6
;(# ( J 1 : - ' ~ (
CALL
J,
l )
~ A C S u H ( ~ , X~, x l ,RANK)
DI, ~ 7
J=I,~4~'JK ,.(j, C,3L~.') =X 1 [ J ) X{J)=,l[j]
27 C C
r.f/~,,,~
vF['Ta)R
a
UF
PP" TNVF~SF.
C t3~'
CAt L k A I v [ " C ( A , X , Y,N) r [,(;n=.l- aLSE . I)tJ c~'4 J : l p~i
P(j,COLM]=Y[J) IF ( J - R A N K ] 3 0 , $ O , 3 1 Y]{J}--Y(J)
30
GO
31
TO
29
IF { {OABS{Y ( J ] ) . G T . E P S ) . A N D . [ . N O T . COtx, T INUE
2Q
C C C
C~ECK FOR A SINGULAR PENCIl_.
NOGOI )NOGO=.TRUE.
USE U FOR STOPAGF.
O0 3~ J = I , N DU 3l W : I , N ! F [ K - C O L M ] 38# 3 8 , 3 Q U(j,K]=P{j,K)
38
GO TO 3 7 U(j,K)=O.ODO COP,'T l RuE CONTINUE CALL GAUSS [ DUMMy ~ U, SCRATCH, ]'CLMo N, ONEMAT, EPS) I F ( O A ~ S [ U [ C O L m , C O L N ' I } . G T . F P S ) C , O TO ~0 K.#ITE(2,~I) "l F O R M A T ( I , S X , ' L A ~ B D A * B - A I S A SINGULAR P E N C I L . ' , / , 1SXo'IHERFFO~'E THERE AreE NO SOLUTIONS TO T H I S P P O R L E F . * ) GO TO .~2
39 37 ].~}
C C
SE' FAP TF~E PENCIL fS NOT STNG!JLAP. vECTORS L)E P2 AND 0 ~ .
C
CONTINUE
FINDING
C O L U MN
C COL ~:COL ~", 1 F I ANF ¢,=F f AfgE"~+ 1 IF(,'~OGO)GO TO 21 CALL ~ / ~ C S u H [ R , Y I , X I o RAN~ ) Ou ~ J=I,N
qO
,~5
G( J,CUL~')=O.UUO
~a
x(Jj:o.nno r:-L. [ r , s ' , ~(J,EOL,~')mxI
[d)
X ( J 1::'(1 [,J} CU~' T I ,U~" Gr TO d 'a ]F(EIA,\,FW . ( ; [ .
3~
2.1 C C
]r,,)ILAIF
IHF
ETA)FTA:FTANF~"
NU~"FJlr~ Ok S('ILL,TI[JNS THAT
EXISI.
C
CuL-~:C CI[. '~'" 1 •:'
I IJ
~ ;
C
C C
~=ILL ~1
I),'
(:UT
u~
K:I,>
: i~
,a q
.~:
,J I =hi*
j j:
T~F P 2 ] N V
C~'
I #l,I
} -
)
,I-CI'~L a+J ] I F [ J l -COLM)
~. ~9
~D
~18, ~le~, q 7
IF ( K - I ) ~ 8 , ~ l g , a 8 U(I,J1)~Q(I,,JI) Q(I,JJ}=Q(I,JI} GO TO ~ 5
T~A~;SFFII4~AT~ONS.
An algorithmfor the generalizedeigen~oblem ~8
U{IwJl]wPfIoJl} PtloJJ)sP(l,Jl) GO TO ~5
a7 ~5 Q~
U(IoJI)=0.000 CONTINUE CONTINUE CALL GAUSS(OUMMYeU°Tt°ICLM°N,ONEMAT,EP$)
CALL GAU88(DUMMYwTI.SCRATCH,ICLMeNeONEMAT~EP$) CALL INVERSE{TI,SCRATCM,ICLM,OUMMYeN) DO 50
I=I,N
ICM=COLM+I 00 %1 J = I C M , N JJ= J- COLM
IF(K-I)S2wS~,S2 $3
52 51 50 ~3 C C C C C C
~(I,JJ]=SCRATCH(I.J) GO TO 51 P(I,JJ)=SCRATCMII,J] CONTINUE CONTINUE CONTINUE
OETEwMI~E P~ FROM ITS INVFRSF. T l USED TO STORE P 2 . HwIGGLE IN H. A~IGGLE IN A.
CALL GAUSSCOUM~YeP,SCRATCM, ICLM,N°ONEMAT,EPS) CALL INVERSE(R,SCRATCMoICLM,TI,N) C C C
DETERMINE H~IGGLE AND AWTGGLE. CALL CALL CALL CALL
C C C
MVULTCTI,R,DUMMY,NoN,N) M~ULT(DUMMY,Q°d°NoN°N) ~UULTCTI°A,OUMMY,N,~,N} MMULT(OU~Yo~,A,N°N°N)
FOR~' Q=UI*02
hl
C C C C C C
O0 bl I = l , ~ ICLM(]}=I CALL ~SOwIIO,N,Ir~COL,ICL~)
EXTRACT T~F ~At~Ix PARTITIONS. ~ l ]~] ~. A[ IN A. ~ IN H . - i ~ In T I ° f i l I:~VERSF*AI IN wII~vA. N TN P. TME~. SFI l ~ h LO~E~ LEFT PAWTTTTONS TO ZERO. ~U=N-~EIG DO 71 I = l , , ~ U TI:MFIG+i O0 ~ J=T,MFIG U(I,J)=H(II,J) ~(IIoJ)=O.ODU IT(]°J)=-A(II,J] 6~ A(II,,I)=O.OPn ~U b5 J=t°~U JJ=~EIG+J b5 RtI,JI=H(IIoJJ) 7] CO~TINu~
C C C
P~I~'T THE FINAL TkANSFOR~FD MATRICES, nWIIF(2,b~) 5~ FN~AI(IoSXo'T~E T~ANSFO~F~ ~ ~ATRIX I S ' ) O0 b~ I = I , ' ~ 55 ~ I I ~ ( ~ , / J [ R ( | , J J , j = I ° ~ ) %0 F O ~ V A T ( / , ~ X , ' T ~ F T~A~SFORMF~ A ~ATRIx IS ° ) DO 57 ]=1,~', ~7 ~ | T ~ ( ~ , ? ) C A C I , J } , J : ] ° ~ ]
C C C
FOO~ ~1
I'~VFkSE*A].
CALL GAUS$(DUM~YoH°SCkAtCM°TNC[}L,MEIG#ONEMAToFPS} CALL I'~VFkSE[~,SCWAICM,T.~COI°CiIMMY°MEIG) CALL u~ULT(DIJ~Y,A,MI~A°~FT~°~FTG,~FIG}
C C C C
Ft.~v APHAI:DLIMY~=M~*RINVA-A~. AGHAI l~'J SCRATCH. CALL MMUL T (LJe~[~JVAeDt~Y,NIJ,~FTG°~EIG) CALL ~AOO(OU~MY, TI ° SCRATEN.~II°~F I~)
~UILO Q~ ] ~
~
LN~ER LEFT P A R T I I I N N OF THE TRANSFORMATION Q3. 1 1 . LURER LEFT P A R T I T I O N OF 03 IN U.
CALL E Q u A T E ( S C W A T C H p u , n U , M F I G } I F ( ~ T A . E Q . 1}GI) TO b7
317
L. F. GODBOUT. JR. and D. JORDAN
318
CALL r ~ U L I ( P , S E R A I C M , D u ~ Y , N U , N U , ~ E I G ) CALL ~ U L T I D u ~ Y , ~ I N V A , ~ C R A T ~ , N U , ~ E I G ° ~ E I G ) CALL r~AD0(LI~SC~AIC~,u°NU,~FT~I EIA=EIA-I GO TO ~ kOIL n C~. ~7 DO 70 I : I , ~ DO l~ J = I , N IF 1 - j ) 7 ~ , 7 ~ , 7 ~ 7Q T! I , j ) : I . 0 D 0 GO T0 72 73 IF ] - ~ E I G ) 7 5 , 7 5 , 7 ~ I1 I , J ) = O . O 0 0 75 GO TO 7 2 ?b IF(J-MEIG)77,/7,75 77 II=I-~EIG Tl(l,J)=U(II,J} 72 CONTINUE 7 0 CONIINUF DETER~II~E Q = I Q 1 * Q 2 % * Q J . G IN DUMMY. CALL ~ U L T ( Q , T t , D U ~ Y , N , N , N ) PPINT ~I
INVERSE*A|
AND THF Q TRANSFORMATION.
202
FO~AI(/,bx,'TME QI I N V F ~ S F * ~ M A T P I X IS') DO IU0 I = I , M E I G I00 ~ Q ] T E [ 2 , 7 I ( h l N V A C I , J I , J = I , M F T G ) ~R]TE(@,7~) 78 F O k M A T ( / , S x , ' T M F Q INANSFO~ATION IS') UU 7 9 I = I , N 79 ~ I T E f ~ , 7 ] ( O U ~ M Y ( I , J I , J = I , N ]
C C C
00NE) ],2 I F ( I E N D .LE. ~WITE(2,6~I b~ FO@MAT(IH)) GO T0 lq ENO
0)SrUP
GAUSS ~LI~I~,AIIO~,J
~OUTINF.
CG~PILPW F)UUWLF PP~CIHTLIi~ SUP~flUTI~;E GAUSS(~MATmMMAT,PN~AT,TC,KDI~,ONF~AT,FPS} l)I~,k~Sl*) 'V A ~ A T ( I 0 , 1 U ) , ~ A I f I 0 , 1 0 ] , P M A I ( 1 0 , I 0 I , I R r I 0 ) , I C ( 1 0 ) , T p ( I ~ ) LOGICAL Oh E~AT I~TTTALIZE ~O~ ANI) COLUMN TNDFX ARPAyS. T ~ , I r l a L I Z k NUa TRA~FOW~ATInN MATPTx TO TMF I ~ N T I T Y .
O0 10 I = ] , K O I ~ i~(II=T IC[I]=I IPCI]=T OC ~0 J = l , ~ ) i ~ IF(I-J}@~,~Om~b ~0 P'~AT(I,J)=I.0DO 25 P~A I ( I , J ] =0. OI)U ~0 C ~ T I ,~HE 10 CONTI!~U~ Pk~F(~ OO
~n
(~URS ELI~INATION USING FULL PIVOTING.
l=l,~I
~
I I ~ I F I A L I Z E THE PlvDT
FLFMFNT.
PI~oT:~AT(I~{II,TCfl)) ICOL~=I FIND
T~E PIVOT
ELFmENI.
DO 50 J = I , K D I ~ D0 O0 K=I,KDIM I F ( D A R $ [ S M A T ( I ~ ( K % , I £ ( J ) ) ) . L E . D A B S ( P I V O T I ] G OTO bO PIVOT=B~AT(]~(~),TC{J)}
bO 50
ICOL~=J CONTINUE CONTINUE CMFCK
THF
~AGNITUDE
OF
1F(DA~S[PlVOT].LT.EPS)GO C
TMF
PIVOT
TO 7 0
ELEMENT.
An algorithm for the generalized eigenproblem CHANGF RO~S
IF NEEDED.
IF(IR(IRO~).EQ.IR(1))GO TO 7% JRO~=IR(1) I~(I)=IR(IRO~) IR(IRU^)=JRO~ CHA~GE COLUMNS IF NEEDED.
7S
]F(IC(ICOLM).EQ.IC(I))GO T0 80 JCOLM=IC(II IC(1):IC(ICOLM) IC(ICOLMI=JCOLM ELIMINATION FOLLOWS.
~U
IF(I.bE.~DI,a)~O rn 70 II=I*l i){' ~u , ] = l l , K I ) l ~ OtIg#~lr,*~ A~ EL~M~I OF T~F ~O~ T~ANSFORMATION MATRIX. RtL~M:~,~AI(IW(J)oIC(1))/PIVGI ALOE, R ~OLIIRLE OF IHP PIVOT ~(;~ TO ANOTMER ROt.
~S ¢5
~AI(r~(j],lC(tl)=O.nnO
O0 qb ~ = ] I p g l ) I M
MvATCIR(j),IC(^))=HMAT(I~(J),ICfKI)-RELE~*MMAT(IR(TI,IC(K)) O0 1.)0 K=l*gt) J "i IFt.i~OT.ONE~AIJAMAT(IRfJ)tIC(K))=A~AT(IR(JI,ICCK)) I -RELEM*A~AT(IR(T),IC(K)) tO0 R~AT(I~(J},K) =P~AI(IR(j)'W)'pELEM*PMAT(IR(T)*K) @0 CONTINUE ~0 CO~[INU~ S(JRT TME MATRICES.
70 CALL ~SORT[~MAT*KDIM*IR~IC) CALL uSOWT{P~AI*KOI~ ' | R ' I P ) IF(.NOI.h~EPAI}CALL ~SORT{AMAT,KDIM~IR,IC) RETURN EW ,D shPTT*~'G wOUII~E FO~ ~ATRTCFS. CO~PILf~ hOi,PLE P~t~ISIUN SUhRUU I II~*~ ~'SURI{ S~AToJi)IM, JR, JC) UI N~r,~ I 0,,, S~,AT( I O#IO),DUMMY(IO°I 0 ) ,JR(TO) ,JC(TO) O0 110 I=I,JC)/~ 0{~ ~IO J=loo~)l ~ 110 I ) u ~ Y ( I , J ) = S ~ A I ( J R { I ) ° J C ( J I ) CALL ~:~!~Al@ ( OUt'Y, S',AT, JDI~° JOT~)
C C C'
SUBROUTIN(
FINDS
THE INVERSE OF A M~TRIX.
COMPILER DOUBLE PRECISION SUBROUTINE INVERSE(UMtT,PHAT,ICOL*TINV,LDIM) DIMENSION U H A T ( I 0 , 1 0 I , P H A T I 1 0 , 1 0 I , T I N V ( 1 0 , I O I , I C O L ( 1 0 I * IPCOL(10),TCOL(|0) EXTRACT a COLUMN FROM PM~T TO FINO A COLUMN OF TMF INVERSe.
Sl0
DO 800 J = I , L D I M DO 8 1 0 IffiI,LDIM PCOL(I)=PMAT(I,J) FIND A COLUMN OF THE INVERSF.
CALL BACSUflfUMATePCOL~TCDL,LDIM)
FORM THE INVERSE FROM ITS COLUMN VECTORS. DO 820 K=|~LDIM 820 TINV(ICDL(K)*J)=TCOL(K) ~00 CONTINUE RETURN END eAC~ SUBSTITUTION ROUTINE. COMPILER DOUBLE PRECISION SUBROUTINE ~ACSUB(URMAIoCOLU~TSOL,JLENGI DIMENSIO~ u R M A T t l 0 , I O ) , C O L U M ( 1 0 ) ~ T S O L ( 1 0 ) FINO IME LAST ELEMENT OF THE SOLUTION VECTOR FIRST,
TSOL(JLENG)=COLUM(JLENG)IUPMATIJLENGeJLENG} C
319
320
L . F . GODBOUT,JR. and D. JORDAN FI~O
ThE
OTMEW
ELE~ENIS
nF
TMF
SOLUTION
VFCTO~.
JJ=JL~G'I Ob ~ 0 J : 1 , J J JA:JLF~G-J 5TOT:0.000 Ja]=ja+1 Dn ~ 0 T:JAI,JLENG KW:JAIejLENG-I B~O SILIT:STOT+UPMAI[j~,K~)*TSOL(KK] 830 T $ O L ( J A ) = [ C O L U M [ J A ) - S T O T ) / I I P ~ A T ( J A , J A ) ~TUR~
~OLTIPLIES
920 910 g00
x MM
MATRIX
~Y
AN MR X PP
COMPILEW OUUPLE PRECISION ~U~W~uT]N~ ~ M U L T C X M A T , Y M A T , I M A T , N N , M M , P P ) DI~E~,SION x M A I ( | 0 , 1 0 ) , Y M A T ( 1 0 , 1 0 } , Z M A T ( ~ 0 , 1 0 ) I N T E G E ~ PP OU 900 I : ~ , ' ~ N OO 910 J : I , P P SU~IqT:0.000 O0 qPO K : I , ~ M 5U~'TnT:Su~TOTeXMAT[T,K)*Y~AT(K,j) Z~AI(I,J):SJMTOT CON~]~UF ~ETuw~ E~D
~ULTIFLIE5
2~I 220
AN N N
A '~AIQIX MY A VFCTDR.
COMPILFP DOUBLE PRECISION SUbroutiNE MATvEC{TMAT,INVFC,OUTVFC,DIMET] OIMES;SIO~' I M A I [ 1 0 , 1 0 ) , O U I V F C ( } 0 ) I~TEGEQ OlmEI ~EAL I , , ' v F C ( 1 0 ] O0 2 ~ 0 I : I , O I M E I SUM~Y=0.0D0 ()U ~ I J = [ , i ) I~ET ~UM~=bU~MY+T~AT(I,J}*INVFC(J) O U I V E C ( I J = ~ U ~MY wETuP~ END
~OOS
T~C~ ~ a l " ~ I C ~ S °
COMPILER DOUBLE. P R E C I S I O N S U ~ P U U I INE ~ A O D (MAT I, MA T2, M A T S U ~ , NR, NC) Q E A L ~AT I ( I O, I O } , M A T 2 [ tO° |t)I , M A T S U M { I0° IO) O0 7 I : I , ~ , R DO 7 J : L , N C 7 MATSUM{ I , J ] :~.A I I { I , J ) ÷ M A T 2 ( I , . J } RETUQN END C C C"
EI~L,#TES ONE '.,ATRIX
TO ANOTMFR.
CO,PILaW rJOUBLE PHECISION SU~POU I I:H" ~_UUATE { IN"~A T , (%UT~A T, N 1 , NP) OI ~E {,,S I ~,~ O O T ~ A T ( 1 0 , 1 0 } RLAL I ~ , V A T ( I Q , ] 0 ) i)O 17 I : 1 , ~ ' ~ 1 OU 17 j--],r~_P 17 otJT~a I { I , J ) = I N~'+AI { ] , J J WE f OPt,
MATRIX.