A numerical solution of the matrix Riccati equations

LIBRARY OF THE

UNIVERSITY OF ILLINOIS

AT URBANA-CHAMPAIGN

510.84

T£63c
*o. ai-30

is

AUG. 5 1.976

ihe person charging this material is re-
sponsible for its return to the library from
which it was withdrawn on or before the
Latest Date stamped below.

Theft, mutilation, and underlining of books
are reasons for disciplinary action and may
result in dismissal from the University.

UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN

reftl

f

R

n

JUL

12

DEC * & 8BT&

1980

L161 — O-1096

510.814

I^63c

no.2U

Engin.

1

00

GINGERING LfBRARy
y OF ILLINOIS

MS

BIBLIOGRAPHIC DATA
SHLET

1. K, port No.

[JIUC-CAC-.DN-72-2U

3. Recipient's Accession Ni

•J. 1 it le ami Suht it Ic

A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS

5- Report Date

January] 20, 1972

7. Autlu>r(s)

Killion Noh

8. Performing Organ i /.at ion Kept.

No

CAC-2U

9. Performing Organization Name and Address

Center for Advanced Computation
University of Illinois at Urb ana- Champaign
Urbana, Illinois 61801

10. Project/Task^ Work Unit No.

11. Contract .'Grant No.

DAHC0i+ 72-C-OOOl

12. Sponsoring Organization Name and Address

U.S. Army Research Office-Durham

Duke Station

Durham, North Carolina

13. Type of Report & Period
Covered

Research

14.

15. Supplementary Notes

16. Abstracts

The eigenvector solution of the time-invarient matrix Riccati equation is
discussed.- The coefficient matrix of the canonical equation is allowed to have
multiple eigenvalues, namely, the matrix could be either derogatory or defective.
The solution of matrix Riccati equation is then calculated from a part of similar-
ity transformation which should reduce the coefficient matrix to the Jordan
canonical form.

1

17. Key Words and Document Analysis. 17a. Descriptors

Ordinary and Partial Differential Equations
; Matrix Algebra

Nonlinear and Functional Equations /

17b. ldcntif iers/Open-Ended Terms
17c. COS ATI Field/Group

18. Availability Statement .-. . . , . ,...,,.

No restriction on distribution.
Available from National Technical Information
Service, Springfield, Virginia 22151

19. Security Class (This
Report )

UNri.ASSIFM-n

21. No. of Pages

38

20. Security Class ( 1 his
Page

UNri.ASSIFIF.D

22. Price
- ______

FOI'M NTlS-iS INCV. J-//I

THIS FORM MAY 1)1- KFPRODUCHP

USCOMM-OC I4832-P72

CAC Document No. 2k
UIUCDCS-R-72-U98

A NUMERICAL SOLUTION OF THE
MATRIX RICCATI EQUATIONS

By
Killion Noh

Center for Advanced Computation
University of Illinois at Urb ana-Champaign
Urbana, Illinois 6l801

January 20, 1972

Submitted in partial fulfillment of the requirements for the degree of
Master of Science in Computer Science in the Graduate College of the
University of Illinois at Urbana-Champaign and supported in part by the
Advanced Research Projects Agency of the Department of Defense and "was
monitored by the U.S. Army Research Office-Durham under Contract
No. DAHCO^ 72-C-OOOl.

Digitized by the Internet Archive

in 2012 with funding from

University of Illinois Urbana-Champaign

http://archive.org/details/numericalsolutioOOnohk

&NGINEEJUNG U8RAA*

ii

ABSTRACT

The eigenvector solution of the time-invariant matrix Riccati equation
is discussed. The coefficient matrix of the canonical equation is allowed
to have multiple eigenvalues, namely, the matrix could be either derog-
atory or defective. The solution of matrix Riccati equation is then cal-
culated from a part of similarity transformation which should reduce
the coefficient matrix to the Jordan canonical form.

Ill

ACKNOWLEDGEMEIW

The author would like to express his thanks to his advisor. Professor
Daniel L. Slotnick, and to the ILLIAC IV Project and the Center for Advanced
Computation under whose employ this research was conducted. A note of
appreciation is also extended to the Department of Computer Science for its
cooperation and support during the time of the author's studies at the
University of Illinois.

Deepest gratitude also goes to Professor Ahmed H. Sameh for suggesting
this thesis topic and for his guidance throughout its preparation. Without
his encouragement and patience, this work would have never been done.

IV

TABLE OF CONTENTS

Page

1. INTRODUCTION 1

2. FORMULATION OF THE PROBLEM 2

3. SOLUTION OF THE MATRIX RICCATI EQUATION 6

k. COMPUTATIONAL METHODS AND NUMERICAL RESULTS 13

LIST OF REFERENCES 18

APPENDIX I. ALGOL PROGRAMS 20

APPENDIX II. ILLIAC IV PROGRAMS 28

1 . INTRODUCTION

It is well known [1,2] that the linear regulator problem with quadratic
cost functional is reduced to the problem of the matrix Riccati equation.

Although the solution of the matrix Riccati equation exists and is
unique [2,3], it is not easy to obtain a numerical solution for the
high-order system.

One of the numerical methods is the eigenvector solution which has been
studied by several authors [U , 5 »6 ] . In their studies, however, the coeffi-
cient matrix of the canonical equation is assumed to have distinct eigen-
values, which is a restriction to the method of eigenvector solutions.
Recently, Martensson [7] discussed the case in which the coefficient matrix
is nondiagonable. Similarly, in this paper, the assumption of distinct
eigenvalues will also be removed so that the coefficient matrix is allowed
to have the multiple eigenvalues, namely, the matrix is either derogatory
or defective. To do this, we explicitly construct the similarity transfor-
mation which should reduce the coefficient matrix to the Jordan canonical
form, either diagonal or nondiagonal, and then the solution of the matrix
Riccati equation is calculated from a part of the similarity transformation
matrix.

Since matrix Riccati equations of fairly large sizes appear in many
engineering applications and since the solution algorithm can be easily
constructed such that it becomes suitable for a parallel computer [8], we
provide in Appendix II the basic codes written in ILLIAC IV language,
GLYPNIR.

2. FORMULATION OF THE PROBLEM

In this section, the source of the matrix Riccati equation and its
relation to optimal control problems are summarized in short for our
further discussions [1,2,3].

Let us consider the linear time- invariant dynamical system

x(t) = Ax(t) + Bu(t)

y(t) = Cx(t) (2.1)

x(tQ)=x0
where A, B, and C are n x n, n x r, and m x n constant matrices, respec-
tively, and x(t) is an n-dimensional state vector, u(t) is an r-dimensional
control vector, and y(t) is an m-dimensional output vector, respectively.
Further, we assume that the system (2.1) is completely observable. The
optimal output regulator problem is then to determine the control u(t) which
minimizes the quadratic cost functional,

V = |yT (t ) IV(t ) + \ Jf (yT(t) Qy(t) + uT(t)Ru(t))dt

t0 (2.2)

where P and Q are m x m symmetric positive semidefinite matrices and R is
an r x r symmetric positive definite matrix. Substituting y(t) = Cx(t)
into (2.2), V is rewritten as

V = \ xT(tf ) CTPCx(tf) + | / (xT(t)CTQC x(t) + uT(t)Ru(t))dt

t0 (2.3)

T T
Since the system (2.1) is completely observable, C PC and C QC are symmetric

and positive semidefinite [l].

The optimal feedback control, therefore, is given by

u*(t) = -R_1BTK(t)x(t) - (2.U)

where K(t) is an n x n symmetric positive definite matrix which is the

solution of the matrix Riccati equation

K(t) + K(t) A + ATK(t) - K(t)BR_1BTK(t) + CTQC = 0 (2.5)
with the boundary condition

K(t ) = CTPC ' (2.6)

Furthermore, the optimal trajectory is the solution of the system of
differential equations

x(t) = (A - BR"1BTK(t)) x(t) (2.7)

x(tQ) = xQ

The minimum cost is given by

V • = |xT(t) K(t) x(t) (2.8)

In addition to the assumption of complete observability we further
assume that the system is completely controllable. Setting P = 0 and
t = °°, the output regulator problem turns out to be the following:
minimize

V = \ /~(yT(t)Qy(t) + uT(t)Ru(t)) dt (2.9)

subject to

x(t) = Ax(t) + Bu(t)

y(t) = Cx(t) (2.10)

x(0) = xQ
Then the optimal control is given by

u(t) = -R_1BTK x(t) (2.11)

where K is the n x n symmetric positive definite matrix which is the solu-
tion of algebraic Riccati equation

KA + ATK - KBR_1BTK + CTQC = 0 (2.12)

The optimal trajectory is the solution of the system

x(t) = Gx(t)

-1 T
G = A - BR B K

(2.13)

x

(0) =

Note that for an optimal control system the matrix G is diagonable and has'
eigenvalues with negative real parts. Therefore, the symmetric positive
definite solution K of (2.12) must be such' that G has eigenvalues with
negative real parts.

Kalman [2] has shown that if the system (2.1) is completely observable
and completely controllable, then

lim K(t ) = K (2.1U)

t -*»

f

We shall observe this equation in the following.

Let K(t) = Z(t)X-1(t) (2.15)

where X(t) and Z(t) are n x n matrices and X(t) is non-singular. Substituting

K(t) = (Z(t) - K(t) X(t)) X_1(t) (2.16)

into (2.5) and setting

X(t) = AX(t) - BR-1BT Z(t) (2.17)

we obtain

Z(t) = -CTQCX(t) - ATZ(t) (2.18)

Therefore (2.5) becomes a system of 2n-dimensional linear homogeneous
differential equations

x(t)

"i(t) 1= \\

z(t) [-c1

-1 T
-BR B

-C QC -A

Z(t)

(2.19)

Since we are interested in increasing the time t, let

t = tf - t

(2.20)

Then

where

X(r)

= ¥

"*(*)

_Z(t)_

Z(t)

-A

-it"

BR B

¥ =

T
C QC

AT

„

—*

(2.21)

with

X(0) = I

Z(0) = CTQC (2.22)

Therefore, the solution will be given by the following matrix exponential

X(f)"

Z(t)

= e

¥t

T
CTQC

(2.23)

Consequently, the solution of the matrix Riccati equation will be obtained
by (2.15) as t-> °°. In the next section, we shall show that two matrices
X(t) and Z(t) are obtained from the similarity transformation which reduces
¥ into the Jordan canonical form. Throughout . the remaining discussions,
the increasing time variable t will be used instead of t.

3. SOLUTION OF THE MATRIX RICCATI EQUATION

We begin with two lemmas for our main result.

Lemma 1. The eigenvalues of the matrix ¥ of (2.21) are symmetric
with respect to the imaginary axis.

Proof. Let I be a 2n x 2n matrix given by

xi =

-I
0

(3.1)

where I is an n x n identity matrix. Then it is obvious that

-1 T

I = I

1 1

T T

I WI = -¥

(3.2)
(3.3)

Since the eigenvalues of a matrix are unchanged by either similarity trans-
formations or transposing, the eigenvalues of W occur in pairs with opposite

signs.

Lemma 2. Let

T =

All A12

A A
LA21 A22J

(3.10

where A , A , A , and A are n x n, m x m, nxm, and m x n matrices,
respectively. Assume that A.^ is non-singular. Then

-1

T T

11 12

T T

L 21 22J

-1

(3.5)

exists iff A - A A " A is non-singular, and is given by

-1

-1

Tll : All " + All " ^2 T22 A21 All
T12 = "All ' A12 T22

-1

(3.6)

T = -T A A
21 22 21 11

-1

T22 = (A22 " A21 All A12)

7

where submatrices T , T , T , and T are the same sizes as A ,
A22' A12' and A21' resPectively-

Proof. Suppose A - A Ann" ^o is non-sin8'ular' From TT " = I.

All Tll + A12 T21 = Xn

(3.7)

lll T12 + ^2 T22 " °

A21 Tll + A22 T21 " °

A21 T12 + A22 T22 = Im
Premultiplying the second equation by A A and subtracting from the

fourth,

T22 " (A22 " A21 All " A12) '

And hence

T12 = ^1 \2 T22

Similarly, from the first and third,
T21 = "T22 A21 ^1 '

-1

-1

Tll == All " + All " h2 T22 A21 All

-1

Hence T is a right inverse of T. In a similar way we can show that T
is also a left inverse of T. Conversely, suppose T is non-singular. Then

0 ? det

A,

^1 "12

LA21 A22J

= det

A22~A21 hi A12.

-1

-det (A^) det(A22-A21A11- A^)

(3.8)

Therefore, A - A A A is nonsingular.

Assuming that W has distinct eigenvalues, we now turn to constructing
the similarity transformation with which W of (2.2l) is reduced to the

8

diagonal form. By lemma 1, W has n eigenvalues with positive real parts
and n eigenvalues with negative real parts. If we construct the similarity-
transformation S, whose columns are the right eigenvectors of ¥, such that
the first n columns are the right eigenvectors of W corresponding to the
n eigenvalues of W with positive real parts, then

s_1ws =

1

X.

*A

2n

(3.9)

,-1

where, of course, the rows of S " are the left eigenvectors of W.
O'donnell [5] and Potter [6] have shown that if we partition S into four
n x n submatrices,

S =

Sll S12

.321 S22J

(3.10)

then the solution of the matrix Riccati equation is given by

-1

K - S21 S±1
provided that S is non-singular.

(3.11)

We shall now remove the restriction that W has distinct eigenvalues and
present an alternative proof to that of Martensson [7] to show that (3.11 )
holds for the case in which ¥ has nondiagonal Jordan form.

Let us consider the general Jordan form J with the similarity transfor-
mation S, which will be constructed later,
S-1¥S = J

where, J =

"i ,

(3.12)

'J

m

and in which

m.

A. 1

1

A. 1

l

". A.

, (m. x m. )

(3.13)

If we denote S by its column vectors s , . . . , s , then

¥ Ls1, s2,

> S2n] = [S1' S2' •••• S2n] J

(3.1M

For any Jordan block J , if i. > 1 then, dropping subscripts of s for

m. l

l

convenience,

(1) . _(D

Ws

= sv 'A.

l

w«<2> ■ *(l) + > .(2)
ws = s + A . s

(3.15)

or

TT (m. ) _ (m.-l) (m. )

¥si=si +A.si

(W - A.I) s(l) = 0
i

(W - A.I)2 s(2) = 0

(3.16)

/ vin. (m. )

(W-A.IJis l =0

Thus, s , 1 < j <m., is a principal vector of degree j associated with

, 0 . 0 . ... . , n , (l ) (m. )

A.. Since S is nonsmgular, the m. principal vectors s , ..., s l

are independent. If we consider all the Jordan blocks J , 1 < i < £, we

m. -

l

have the 2n principal vectors s , s ,

, s^ which constitute the

2n

similarity transformation S in (3.12). Further, we rearrange the columns
of S such that the first n vectors are the principal vectors corresponding
to the n eigenvalues with positive real parts of W.

10

Let us now split J of (3.12) into

J =

: J2^

• \

(3.17)

where J and J are diagonal with the eigenvalues of W with positive and
negative real parts respectively, and E and E are matrices with one's on
the super-diagonals and zeros elsewhere. To associate with (3.6), let us
partition S and S as

sll

S12

s =

S21

S22

-

_

s-1-

Tll

T12

T21

T22

(3.18)

where all submatrices are n x n. Note that, by lemma 2, it is necessary to
assume that S and S - S S S are nonsingular.

Then

wt

e = e

SJS 1t

= Se^s'1

0m

Sll S12

e x 0

e^ 0

—

T
11

T
12

S21 S22

n J2t
0 e ^

— —

0 ^

T
21

T
22

S11eJlteEltT11 + S12eJ2teE2tT21

S11eJlteEltT12 + S12eJ2teE2tT22

q eJlteEltT + S eJ2teE2tT S eJlteEltT + S eJ2teE2tT
b21 L 111 22e d 6 ^ 121 b21 ± X ^2 b22 d & * X22

(3.19)

11

Since the elements of Jp have negative real parts, e 2 -^ 0 as t + » and
thus the terms having e 2 in (3.19) vanish as t -*■ °° . Therefore,

X(t)
Z(t)

¥t

T
C QC

„ J_ t iii_ tm ^ J. t E., t_

Slie X 6 1 Tll Sll6 X e 2 T12

21 X e X 111 b21 ^ X12

ri

T
C QC

Su eV eV (T1X + T12 CTQC )
S21 eV eV (T^ + T CTQC)

(3.20)

Therefore the solution of the matrix Riccati equation is given by

K = lim Z(t) X_1(t)
t-*»

= lim [S eJlt e¥ (T +T CTQC ) ] [S eJlW (T +T CTQC) I"1

t-K»

■ S21 Sll

-1

(3-21)

provided that the indicated inverse exists. However, we shall show that
only the nonsingularity of S and S - S S S p is needed for the

stationary solution of the matrix Riccati equation. From (3.12),

T T

11 12

T T

21 22

-1 T
-A BR B

T T

C QC A

Sll S12

21 22

= J

(3.22)

Equating the (2,1) elements on each side of (3.22),

(-T21 A - T2£1 CTQC) Su + (T BR^B1 + ^ AT) 8^ = 0

12
or

■ -T21 ASn + T22 AT S21 + T21 BR^ S21+ T22 CTQC S^ = 0 (3.23)

Substituting T21 = - T22 S21 S^"1 into (3.23),

T22S21Sll'1+ T22ATS21 " T22S21Sll"lBR"lBTs21 + ^2°^° Sll=°

(3.21+)

Premultiplying by T22" and post multiplying by S _1 ,

S21Sll"lA + ^^l3!!"1 " ^l5!!"1^"1^^!3!!"1 + <?<* = ° .0.25)

which is to be proved.

13

k. COMPUTATIONAL METHODS AND NUMERICAL RESULTS

The major three steps for computing the solution of the matrix Riccati
equation are as follows:

1) Compute eigenvalues and eigenvectors of W,

2) Compute principal vectors of ¥ if W is defective,

3) Compute S S " from S.

For computing the eigensystem, we used two methods: (i) the QR
algorithm [9] with an inverse iteration routine [10], and (ii) Jacobi-like
method [ll] for finding the eigenvalues and eigenvectors. Their perfor-
mances are satisfactory and two numerical test examples are given in Test
Result 1 and Test Result 3. When the matrix W is nondiagonable, Test Result
2, the first method requires the evaluation of the principal vectors as
well [12]. This poses some numerical difficulties as we are essentially
trying to solve the ill-conditioned system of equation

(W - Al)k x = x., k > 1
where ^ is an approximation of true eigenvalues A. Thus, the Jacobi-like
method is more attractive in this case for providing an approximate eigen-
system that is suitable for engineering purposes.

Appendix I contains the listings of the ALGOL programs used on the B-6500
Appendix II contains the listings of two ILLIAC IV programs, namely, the
inverse iteration routine and a routine for solving a system of complex
linear equations, both written in GLYPNIR. ILLIAC IV programs for the QR
algorithm [13, 1^] and the Jacobi-like algorithms [15, l6] are already
available.

Test Result 1.

ik

A =

1

0

0

0

0

1

0

-1

0

0

0

0

-1

0

0

0

0

1

0

-1

0

0

0

0

-1

p =

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

B =

1

0

0

0

0

0

0

1

0

0

0

0

0

0

1

Q =

0

0

0

0

0

0

10

0

0

0

0

0

0

0

0

0

0

0

10

0

0

0

0

0

0

c =

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

R =

1

0

0~

0

1

0

0

0

1

The numerical example was tested by the QR algorithm and the inverse
iteration method. The ten distinct eigenvalues of W are

± 1.

± 1.728760U77185 ± 1.577533767573i,

± 1.35319578U0U6 ± 1.1537U989928Ui.
We choose the five eigenvectors corresponding to the five eigenvalues
with positive real parts and computing K = S S from them gives the
solution K as

15

1.262782609

2.1+9^009759

-0.819173651

0.668267901

-O.UU3608958

2.1+9^009759
7.^ 3 5^51161+

-1.8257^1858
1.1229101+32

-0.668267901

-0.819173651
-1.82571+l858
1.63831+7303
1. 8257^1858
-0.819173651

0.668267901
1.1229101+32
1.8257^1858
7.1+ 35^51161+
-2.1+91+009759

-O.Ul+3608958
-0.668267901
-0.819173651
-2.1+91+009759
1.262782609

-1 T
All the eigenvalues of G = A - BR B K have negative real parts and hence

the above solution K is acceptable.

Test Result 2.

A =

-3 '2

-2 1

B =

C =

1 0
0 1

P=Q=

0 0
0 0

The mat

rix

—

3

-2

0

0

2

-1

0

1

¥ =

0

0

-3

-2

0

0

2

1

R =

CO

is defective, that is, W has four eigenvalues ±1 (double) and each +1 and
-1 correspond to a quadratic elementary divisor, respectively. This
example was tested both by Jacobi-like method and by a combination of the
QR algorithm and the inverse iteration method.

16

The Jacobi-like method gives K as
-10.0000081* 6.00000U3

6.000003^ -k. 0000027

and the QR and inverse iteration method give K as
-9. 9999917 5.99999^7

5.9999953 -3.9999983
The exact solution is

~Vlo 6
6 -k

For computing the principal vector by the QR and inverse iteration, we
chose the perturbation for the true eigenvalues as much as 10

Test Result 3.

A =

6

k

k

1

h

6

1

k

h

1

6

k

1

h

h

6

B =

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

1

C=R=

1

0

0

0

0

1

0

0

0

0

1

0

0

0

0

1

P=Q=

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

The exact eigenvalues of W are

±15, ±5, ±5, ±1
Since the matrix W is derogatory, +5 and -5 correspond to a two-dimensional
subspace, respectively. The computed eigenvalues by Jacobi-like method are

IT

lU. 99999999930
k. 9999999991+9

-U. 99999999982
0.99999999997

-1^.9999999984
4. 99999999981
_1+. 999999999^8.

-0.99999999995

The computed solution K is given by

-1.999999999811 1.9999999999^2 1.999999999913
1.999999999884 -2.000000000015 -2.000000000000
2.0000000002T6 -2.000000000^07 -2.000000000378

-2.00000000036U 2.000000000509 2.000000000I180

-2.000000000015
2.000000000087
2.000000000^80

-2.000000000568

in -which the absolute error is less than 10

-9

18

LIST OF REFERENCES

[I] M. Athans and P. L. Falb, Optimal Control, McGraw Hill, New York,
1966

[2] R. E. Kalman, Contributions to the Theory of Optimal Control,
Boll. Soci. Mat. Mexicana, 5 (i960), pp. 102-119

[3] R. S. Bucy and P. D. Joseph, Filtering for Stochastic Processes
with Applications to Guidance, Interscience, New York, 1968

[h] A. G. J. MacFarlane, An Eigenvector Solution of the Optimal Linear
Regulator Problem, J. Electronics and Control, 1^4- (1963), pp.
61*3-6^5

[5] J. J. O'donnell, Asymptotic Solution of the Matrix Riccati Equation
of Optimal Control, Proc . Fourth Annual Allerton Conf . , 1966,
pp. 577-586

[6] J. E. Potter, Matrix Quadratic Solutions, J. SIAMAppl. Math., ll+
(1966), pp. U96-5OI

[7] K. Martensson, On the Matrix Riccati Equation, Inf. Sci., 3 (l97l),
pp. 17-^9

[8] D. L. Slotnick, et al, The ILLIAC IV Computer, IEEE Trans, on
Computer. C-17 (1968), pp. 7^6-757

[9] R. S. Martin, G. Peters, and J. H. Wilkinson, The QR Algorithm for
Real Hessenberg Matrices, Num. Math., lk (1970), pp. 219-231

[10] J. M. Varah, The Calculation of the Eigenvectors of a General Complex
Matrix by Inverse Iteration, Math. Comp. , 22 (1968), pp. 785-791

[II] P. J. Eberlein and J. Boothroyd, Solution to the Eigenproblem by a
Norm Reducing Jacobi Type Method, Num. Math., 11 (1968), pp. 1-12

[12] H. J. Bowdler, R. S. Martin, G. Peters, and J. H. Wilkinson,

Solution of Real and Complex Systems of Linear Equations, Num. Math.,
8 (1966), pp. 217-234

[13] M. Ogura, A Parallel Scheme for Reducing a Real Matrix to the Upper
Hessenberg Form, Center for Advanced Computation Document No. 11,
University of Illinois, 1971

[ill] M. Ogura, The QR-Algorithm, Center for Advanced Computation Document
No. 12, Univ. of 111., 1971

19

[15] W. Bernhard, ILLIAC IV Codes for Jacobi and Jacobi-like Algorithms,
Center for Advanced Computation Document No. 19, Univ. of 111., 1971

[l6] A. H. Sameh, On Jacobi and Jacobi-Like Algorithms for a Parallel
Computer, Math. Comp. , vol. 25 (l9Tl), pp. 579-590

20

APPENDIX I

ALGOL PROGRAMS

21

I TWO ALGflL PROCEDURES' RlCCATU AND RICCATI2. FOR COMMUTING THE OOOOOlOO?

I! SOLUTION OF THE MATRIX RICCATI EOUATION ARE LISTEO BtLOWt FOR 00000200?

% "RICCATlU. THE PROCEDURE "EIGENw (BY EflERLEIN AND bUOTHKoyD) WAS 00000300?

% USEU WHILE A COMBINATION OF "HQR" (BY MaRTINi ET AL) AND 00000400?

I "INVITErATION" WAS APPLIED TO «RICCATI2% LINEAR EqUATIUn SOLVER OOOOO5OO?
% (BY BOWDLER* ET AL ) WAS ALSO USED TO PROVIDE THE INytRSE OF A MATRIX 00000600?

% ANO TO CALCULATE THE PRINCIPAL VECTORS fROM ( ( W-LAMflUA*! )**K>*X«0. 00000700?

00000800?

00000900?

PROCEDURE RlCCATIl(N.L»M»A»B#C»P»Q»R«TMX»W»K)l 00001000?

VALUE N»L#M,TMXI 00001100?

INTEGLR N,L.M.TMX| 00001200?

ARRAY A.e(C.P»Q»R»W.K£i»n| 00001300?

COMMENT SOLVES THE MATRIX RlccATI EQUATION . OOOOUOO?

KA*(AT)K-KR(RINV)(BT)K*(CT)QC«o, 00001500?

A Is N*N» R IS N*L» C IS M*N. P AND 0 ARE M*M. R IS L*L OOOOUOO?

MATRICES* RESPECTIVELY* FORM 2N*2n MATRIx W ANU COMpUTE 00001700?

EIGENVALUES ANU EIGENVECTORS OF W* WHERE FOUR N*N BLOCKS 00001800?

Op w ARE C 1 • 1 >«-A» Cl»2)«8(RlNV>(BT)» ( 2» 1 )■( C ' >QC» AND 00001900?

(2»2)"(AT). REARRANGE EI GEN VEcTOH ( Pr INC I PAL VEcT0R> MATRIX 00002000?

SO THAT THE pIHST N COLUMNS CORRESPOND TO TO THE EIGENVALUES 00002100?

WlTH P0SITIVE REAL rArTS. dENoTINq ThE UPP^R HaLf uF N COLUMN 00002200?

VECTORS BY SMTXl AND THE LOWER HALF BY SMTX2* I HE SOLUTION 00002300?

K IS THEN GlvEN BY «■( SMTx2) ( SMTxl INv ) . THE PRUCEOUrE 00002400?

"FI6EN" IS USED For EIGENPrOBlEM qF Wi 00002500?

00002600?

REGIN COMMENT FIRST OECLARE PROCEOuRESl 00002700?

00002800?

PROCEUURE SETUP(NtL»M.A»B»C»Q»RlNV.W)| 00002900?

VALUE N.L.M, INTEGER N.L»M) 00003000?

REAL ARRAY A»B»C»Q.RINV»W[l#l]j 00003100?

COMMENT THE PROCEDURE SETS UH ThE 2N*2N mATRIX W, rInV Is INVERSE OF RI00003200?

«EGIN 00003300?

REAL ARRAY RINVBTC1 lL»l IN].QCC1 IM»1 IN] » 00003400?

REAL SUMj 00003500?

INTtGER I.J.KJ 00003600?

FOR Ij"l STEP 1 UNTIL N DO 00003700?

FOR J l ■! STEP 1 UNTIL N DO 00003800?

BtGIN 00003QOO?

W[I,J]I«-A[ I. J]} 00004000?

W[N*I,N*J]IbA( J.I1I 00004100?

ENO) 00004200?

FOR 1 1 » 1 STEP 1 UNTIL L 00 00004300?

FOR J»«l STEP 1 UNTIL N DO 00004400?

BEGIN StJMl.O.Ol 00004500?

FOR K,«l STEP 1 UNTIL L DO 00004600?

S|jM|«SUM*RINVCl»K]*BC J#KJ1 00004700?

RINvBTC I.J]l«SUM» 00004800?

ENOi 00004900?

FOR I|M STEP 1 UNTIL N 00 00005000?

FOR Jl»l sTEp 1 UNTIL N 00 00005100?

BtttIN SuHlaO.Ol 00005200?

FOR Ki«l STEP 1 UNTIL L DO 00005300?

SuH«"SUH*B[IsK]*RINVBTtK#J]| 00005400?

Wt I,N*J]l«SUMl 00005500?

END| 00005600?

FOR I|»l STEP I UNTIL M DO 00005700?

22

FOR Jl«l STEP 1 UNTIL N DO

BEGIN SuMl.O.Oi

for K|«i step i until m do

SuMl"SUH*OtItK]*CCK.JJI
QCC I # J) I -SUM l
ENO|
TOR I|«l STEP 1 UNTIL N DO
FOR J|«l STEP 1 UNTIL N DO
BtG]N SuNlsO.Ol

FOR K|»l STEP 1 UNTIL N DO
StjMI»SUM*CtK. i J*«C t K. J3 1
W[N*I.J]|"SUMI
ENDI
ENO OF SETUP|

COMMENT INSERT THE PROCEDURES INNPRODtDECOMPOSE. AND ACCSOLyEi

PROCEUURE lNVERSE(N»AtX)l

VALUE N| INTEGER Nj

DOUBLE ARRAY A»XC 1 • 1 1 I

COMMENT THE PROCEDURE COMPUTES THE INVERSE OF AN N*N rUl

UNSyMMcTRIc MATRIX A. USINq ThE PROCEOURgS "DECQMpUsE* AND
"ABSOLVE". WRITING ON X|
BEGIN

INTtGER I.J.IT.D2I
REAL Oil

ARRAY AA'BB.BU I N* 1 iNJ.P[llNll
FUR I • * 1 STEP 1 UNTIL N DO
1 UNTIL N Oo

BCl»J]l«ltO ELSE B(I»J]l«O.OI
1 UNTIL N DO
1 UNTIL N DO

FuR Jl«l

IF J«I

FUR I,«l

FUR J t -1

STEP
THEN
STEP
STEP
AACT.Jlt-ACI.J1l
DtCOMpOSE(N»AA.Dl.02»P)l
AtCSULVE(N.N.A,AA.P»B.X»BB»IT)|
ENO INVERsEl

COMMENT iNjSfRT THE PRnCEDURE EIGENi

PROCEUURE MuLT(Nl.N2iNJ»A.B»C)l
VALUE N1»n2.N3I INTEgEH N1*N2»n3>
DOUBLt ARRAY A.B.Ctl»lJI
REGIN

INTtGER I.J.KI

REAL SUMI

FOR I l » 1 STEP i

FOR J t • 1 STEP 1

BEGIN SUM|«0.0|

FUR K,"t STEP

SUMI»SUM*ACI»K]*BCK,J]|
CI I .J] l«SUMI
ENO»
FNO MULTl

ARRAY HiNv»SMTXl,SMTX2lllN.llN].WW.TEMP,VECT[llN*N»llN*NJ|
INTLGEH ArRAy search(1»nii
INTtGER N?» I. J- INDEX'

UNTIL
UNTIL

Nt
N3

DO
OO

1 UNTIL N2 OO

00005800?
00005900?
00006000?
00006100?
0000*200?
00006300?
OOOO64OO?
00006500?
00006600?
00006700?
00006800?
00006900?
00007000?
00007100?
00007200?
00007300?
00007400?
00007500?
00007600?
00007700?
00007800?

00007900?
00008000?
00008100?
00008200?
00008300?
00008400?
00008500?
00008600?
00008700?
00008800?
00008900?
00009000?
00009100?
00009200?
00009300?
00009400?
00009500?

00009600?
00009700?
00009800?
00009900?

00010000?
00010100?
00010200?

00010300?

00010400?
00010500?

00010600?
00010700?
00010800?
00010900?

00011000?
oooiuoo?

00011200?
00011300?
00011400?

23

N2I«N+N| 00011500?

INVtRSE(L.R»RlNV)l 0001U00?

SETUP <N,L#M»A»B»C»Q»RINV»W)J 00011700?

TOR It*l STEP 1 UNTIL N2 DO 00011800?

FOR Jl "1 STEP 1 UNTIL N2 DO 00011900?

TEMPI I. JJ t«M[ I.J] I 00012000?

EIQtN(N2»TMX,TEMP,VECT)l 00012100?

COMMENT SEARCH N COLUMNS OF VECT CORRESPONDING TO N EIGENVALUES 00012200?

WITH POSITIVE REAL PARTS. SMTXl AND SMTX2 ARi IT* UPPER 00012300?

AND LOWER HALVES* RESPECTI VELY| 00012400?

INDLX|.0l 00012500?

FOR I » » 1 sTEp 1 UNTIL N2 DO 00012*00?

I* TEMP[I»I1 GTR 0 THEN 00012700?

BLGIN 00012800?

INDEX|«INDEX*1 , 00012900?

SEArChI INDEX] tail 00013000?

ENO| 00013100?

FOR J|»1 STEP 1 UNTIL N DO 00013200?

FOR I, al STEP 1 UNTIL N DO 00013300?

BEGIN 00013400?

SMTxltliJ) |«VECTC I » SEARCH! J]]j 00013500?

SMTX2tI.J]l«VEcTtN*I»SEARcHCJ}]l 00013600?

ENDI 00013^00?

INVLRSE(N, SMTXl. TEMP)I 00013800?

MULHN»N»N»SMTX2.TEMP»K)| 00013900?

END RiCCATH» 00014000?

00014100?

00014200?

PR0CEUURE RICCatI2(n.L»M»a»B»C»P»0»R»W»K»0EFECt»DER0(»AIE)I 00014 300?

VALUE N.L.Mf 00014400?

INTEGtR N.L.Ml 00014500?

ARRAY A»3.C.PfQ»R*H*Ktl'l)l 00014600?

INTEGtH ARRAY DEFECT, DEKUGATEt 1 J J 00014700?

COMMENT SnLvES THE MATRIX RICCATI EQUATION 00014800?

KA*(AT)K-KB(HINV)(BT)K*(CT)QC«0. 00014900?

A Is N*N. b IS N*L» C IS M*N. P AND Q ARE M*M, R IS L*L 00015000?

matrices* respectively. form 2n*2n matrix w and compute 00015100?

eigenvalues anu eigenvectors of w, where four n*n blocks 00015200?

OF W ARE (l.l)«-A. (1.2)«B(RINV)(BT)* (2*1)«(CU0C AND 00015300?

(2»2)«(AT). REARRANGE El GEnVEcTOR( PRlNC I PAL VECTOR* MATRIX 00015400?

SO THAT THE FIRST N COLUMNS CORRESPOND TO TO THE EIGENVALUES 00015500?
WITH POSITIVE HEAL PARTS. DENOTING THE UPPER H*LF Up N COLUMN 00015600?

VECTORS BY SMTXl AND THE LOWER HALF BY SMTX2, THE SOLUTION 00015700?

K IS THEN GIVEN BY K«( SMTX2) < SMTXl INV ) . THE PROCEDURES 00015800?

"hShLD"* "hOR** AnD "INVITEKATION" ARE USEq FOK EH»ENPRORlEM 00015900?

OF Wl 00016000?

00016100?

REGIN COMMENT FIRST DECLARE PROCEDURES! 00016200?

00016300?

PROCEUURE SETUP(N»L*MtA*B,C»Q*RINV*W)| OOOI64OO?

VALUE N.L.Mf INTEGER N.LtM* 00016500?

REAL ARRAY A » B . C » Q . RIN V • W [ 1 i'l ) I 00016600?

COMMENT THE PROCEDURE SETS UP ThE 2N*2N MATRIX W. RINV IS INVERSE Of R! 00016700?

QEGIN 00016800?

REAL ARRAY R I NVBT[ 1 I L« 1 I N J . QC [ 1 I M» 1 1 N ] I 00016900?

REAL SUM1 00017000?

INTLGER I.J.KI 00017100?

2k

FOR I t -1 STEP 1 UNTIL N DO
TOR J|»1 STEP 1 UNTIL N DO
BtGIN

M[I.J]I«-ACI»J]|

WtN+I.N*J]l«ACJ#ni

ENDI

Fo« i»«i step i until l Do

FOR Jial STEP 1 UNTIL N DO
BtGIN SuM|aO»0|

FOR K|«l STEp 1 UNTIL L DO
SuM|«SUM*RJNVt I»KJ*B[ J»K] J

RINvBTtI»J]l«SUMI

END|

FOR I|»l STEP 1 UNTIL N 00
FOR J t"l STEp 1 UNTIL N DO
BtGIN SuMl-0,01

FOH Kt«l STEP 1 UNTIL L 00

SuM,«SUM*BU.K)*RlNVBT[K.JJ|
W[I.N*JJ|«SUmI
ENO|
FOR Iial STEP l UNTIL M 00
FOR J|"l STEP 1 UNTIL N DO
BLqIN SUMI.O,0|

FOR Kf.l STEP 1 UNTIL M DO

SUM|.SUM*Q[ I.K]*C[K»Jj I
«CC If j] I»SUMJ
ENDl
FOR Il»l STEP 1 UNTIL N DO
FOR J t "1 STEP 1 UNTIL N DO
BtGIN SUMl«O.0l

FOR K|b1 STEP I UNTIL M DO
SuM|«SUM*C(K.n*UC[K»Jj;
W[N*I,j) l«SUMI
END>

end of setupi

comment Insert the procedures innprod. decompose* and accsolvei

proceuure inverse(n»a#x)|

VALUE N| INtEgEH N|

nOUBLt ARRAY A.XU.lJI

COMMENT ThE PROCEOURE COMPUTES THE INVERSE OF AN N*N RtAL

unsymmetric hatrix a. using The procedures "DECQmpu$e"
waCcsolvE". writing on XI

BEGIN

INTtGER I.J.IT.D2I
REAL nl |

ARRAY AA'BB.BUtN.l fNJ»Ptl tNJ|
FUR I i - 1 STEP 1 UNTIL N DO

DO

B[I»J]nO«OI

AND

FuR J I -l

IF J-I

FUR I,-!

FUR J i ■ 1

END

step

then
step

STEP

AAt !ij]|iA[I,jii
ntCnMp0sE(N.AA,0l.02#P)|
ACCSULVE(N#N.A.AA.P»B#X»BB»IT)|

i until n

fl[I. J] i-l .0 ELSE
1 UNTIL N DO
1 UNTIL N

00

00017200?

00017300?

00017400?

00017500?

00017600?

00017700?

00017*00?

00017900?

00018000?

00018100?

00018200?

00018300?

00018400?

00018500?

00018600?

00018700?

00018800?

00018900?

00019000?

000i9i00?

00019200?

00019300?

00019400?

00019«00?

00019600?

00019700?

00019800?

00019900?

00020000?

00020100?

00020200?

00020300?

00020400?

00020500?

00020600?

00020700?

00020800?

00020900?

00021000?

00021100?

00021200?

00021300?

00021400?

00021500?

00021600?

00021700?

00021800?

00021900?

00022000?

00022100?

00022200?

00022300?

00022400?

00022500?

00022600?

00022700?

00022800?

25

00022900?

Comment InseRt the procedures hshld and hori ooo23ooo?

00023100?

PROCEDURE MuLT(Nl»N2*Ni'A»B«C)l 00023200?

VALUE N1»N2*N3| INTEGER N1.N2.N3j 00023300?

DOUBLE ARRAy A.B.CU.Ul 00023400?

REGIN 00023500?

!NTtOER I»J»K» 00023600?

REAL SUMS 00023700?

FOR Il«l STEP 1 UNTIL Nl DO 00023800?

FOR Jl«l STEp 1 UNTIL N3 00 00023900?

BEGIN SUM|«0,0I 00024000?

FUR K|«l STEP 1 UNTIL N2 DO 00024100?

SUM i * SUM* A [ I»K]*B[K«J)| 00024200?

CI I • Jl l.SUMI 00024300?

ENOI 00024400?

END MULTl 00024500?

00024600?

PROCEUURE lNVlTERATION<N»w»X»KK»P>| 00024700?

VALUE N, TNTEGER N.KKj 00024800?

DOUBLE ARRAy W.Xtl.lli 00024900?

ARRAy P[l]l 00025000?

CnMMENT SnLVES WX.O BY INVERSE ITERATION* WHERE W IS AN N*N REAL 00025100?

MXTRlX. WRITING SOLUTION ON X. THUS. THE PROCEDURE COMPUTES 00025200?

ThE EIqENvECTOKS CORRESPONDING TO THE GIVEN EIGENVALUES. 00025300?

N Is THE ORDER OF M AND KK IS THE NUMBER OF iTtRATiONS 00025400?

NEEqEO FOR ITERATING EACH EIGENVECTOR! 00025500?

BEGIN 00025600?

INTfcGER I.J»D2»R.L»MAXIN0EX| 00025700?

DOUBLE DUNAXVALi 00025800?

DOUBLE ARRAY MN( 1 1 N. i | NJ »B.B8[ i | N» 1 | U | 00025900?

LABEL INVIT.OUTI 00026000?

FOR Il"l STEP 1 UNTIL N DO 00026100?

FOR Jl ■ l STEP 1 UNTIL N DO 00026200?

ww [ I » J ] t»ri t I » J 3 1 00026300?

DECUMpOsE(N»WW»Dl»02»P)l 00026400?

FOR Il«l STEP 1 UNTIL N DO 00026500?

BII,1]»81,0» 00026600?

KK I "1 1 R| ■ I j 00026700?

INVITI 00026800?

ACCSOLVe(N»R»W»WW»P»B»X.BB*L)I 00026900?

MAXVALUOj 00027000?

FOR Ii"l STEP 1 UNTIL N DO 00027100?

If ABs(X(I*l]> GTR MAXVAL THEN 00027200?

BEGIN 00027300?

MAXvALl«ABS(Xr I#i])l 00027400?

MaXiNDEXI«P 00027500?

ENDi 00027600?

MAXvA|_l«X[MAXlNDEX»ni 00027700?

FOR I|«l STEP 1 UNTIL N DO 00027800?

XlI»l]l«X[Iil]/MAXVALI 00027900?

FOR III} STEP 1 UNTIL N 00 00028Q00?

I* ABS(BCI»1]-Xtl»l]) GTR P"12 THEN 00028100?

BLGIN 00028200?

IF KK GTR 10 THEN GO TO OUTj 00028300?

FOR J|-l STEP 1 UNTIL N DU 00028400?

BtJ,l]l"X[j.i ]l 00028500?

26

KK»»KK+U
GO TO INVITI
END)

nuTi

ENO INVREsElTi

INTtGER I.J»JJ»CNT.N2»INDEX#ITC0UNT|

ARRAY TEMpfORGNtSHTX*veCTORtt|N«N«l|N4N]*RlNViSHTXUSHTX2tl|N«l|N]»

EIGR»EIGI»!NTCH(1IN4N)I
REAL PERTURB. MPUiERl
INTtGER ARRAY SEARCHllIN)!
N2l"N*N|

INVk.RSE(L.R.RlNV)l
SETUP(N,L,M,A»B,C.Q.HINV.H)|
TOR I|«l STEP 1 UNTIL N2 00
FOR J|"l STEP 1 UNTIL N2 00

HSHLO(i,N,,TEMP)|
HQR(N2»TEMP.I

N2
N2

DO
00

N2 00

H

EIQR.EIGI* CNT)|
FOR I | »1 STEP 1 UNTIL N2 00
FOR J I - 1 sTEp 1 UNTIL N2 00

ONGHCli J]l«NCItJ]l
COMMENT COMPUTE EIGENVECTORS OR PRINCIPAL VECTORSi
FOR Jjl-1 STEP 1 UNTIL N2 00
BEGIN

FUR I t»l STEP l UNTIL

FUR Ji-i STEP 1 UNTIL
*[I»J] ««ORGW[I,J]l

FUR I » - 1 STEP 1 UNTIL

•»tIiI]l"WU.i]-ElGR[ JJ]*PERTURB|

JJ«DEFECT[JJ] then
BEGIN

FOR Ii-1 STEP 1 UNTIL N2 00

W[I.I J»«Wti,IJ*pERTURB*MpLIERl

MuLTCN2.N2»N2»W»W,TEMP)l

Ff)R Ital STEP 1 UNTIL N2 00
FOR Jlal STEP 1 UNTIL N2 DO
*t I«J]i«TEMP[ I#J)I
END |
IF JJaOfROGATEUJ] THEN

FOH Ii-l STEP 1 UNTIL N2 00
Mf I • 1 1 |»WC 1. 1 1 ♦PERTURB* J Jl

INVITFRATI0N(N2#W.VECT0R»ITC0UNT»INTCH)»
FUH I t a.j STEP i UNTIL N2 DO

FUR Jiaj.2 DO

SMTx(l»JJ]|avEcTORtI#l Ji
END JJ>

COMMENT SEARCH N COLUMNS OF SMTx CORReSpONOiNG TO
WITH POSITIVE REAL PARTS, SMTxl AND SMTx2
AND LOWER HALVES* RESPECTI VELYl
INDtXiaQi

FOR ll«l STEP 1 UNTIL N2 00
l> ElGRfl] GTR 0 THEN
BLGIN

iNOFXtalNUEX^l |
SEArcmC INDEX] I.I >
ENDl

N tlGENVALUES
ARt ITS UPPER

00028600?
00028700?
00028800?
00028900?
00029000?
00029100?
00029200?
00029300?
00029400?
00029500?
000296O0?
00029700?
00029800?

00029900?
00030000?
00030100?
00030200?
OOO3O1OO?
00030400?
00030500?
00030600?
00030/00?
00030800?
00030900?
00031000?
OOO3UOO?
00031200?
00031300?
00031400?
00031500?
00031600?
00031700?
00031800?
00031900?
00032000?
00032100?
00032200?
00032300?
00032400?
00032500?
00032.00?
00032700?
00032800?
00032900?
00033000?
00033100?
00033200?
00033300?
00033400?
00033500?
00033600?
00033700?
00033^00?
00033900?
00034000?
00034100?
00034200?

27

FOR Jl«i STEP 1 UNTIL N DQ 00034300?

FOR Ii"l STEP 1 UNTIL N DO 00034400?

BEGIN 00034500?

SMTXltI.J]l«SMTX[I»SEARCH[Jlll 00034600?

SMTX2(I,J)I«SMTX(N+I'SEARCH(J])J 00034700?

CNDI 00034800?

INVtRSE(N,SMTXl»TEMP)l 00034900?

MULT(N»N»N»SMTX2.TEMP»K)| 00035000?

END RlCCATl2' 00035100?

28

APPENDIX II

ILLIAC IV PROGRAMS

29

* TWO ILLIAC IV GLYPNIR PROGRAMS. "ClINSYs" AND "INVITE-RATION" ARE

f LISTED BELOW, "CLlNSTS" IS USED TO OBTAIN THE INVERS* OF A MATRIX

% AND ThE SOLUTION Or A COMPLEX SYSTEM OF UNSYMMETRIC LINEAR EQUATIONS.

% MNVITERATION" IS USED TO CALCULATE THE PRINCIPAL VECTORS OF DEGREE

* K FKOM ((WLAMBDAM)**K)*XnO.

SUBROUTINE cLINSyS(CINT N.ClNT R#PCPOINT A'PCPOINT x»PCPOINt B)J
BEGIN x THE SUBROUTINE SOLVES A COMPLEX SySTEM OF UNSyMMETKlC LINEAR

x Equations ax»r. a is n*2n and b is n*2R matrices* respectively.

* ThEir qEnErAi EiEMENTS ArE At I »2J-H*I* At I »2J] AND BCp2j-ll*

* I*B[I,2J], RESPECTIVELY, R IS THE NUMBER OF RIGHT-HAND SIDES,

SUBROUTINE INNpROOf CREAL CWCREAl C2»PREAL AK»PREAL BK»C«EAL Dl#

CREAL D2)»

begin * accumulatEs the sum of products and aods it 10 the initial
% value (c1.c2j.

CHEAL Sl2|

S12,»CUC2|

Si2l»Sl2*ROwSUM(AK*BK)i

01 I .Si 2 t

D2I.S12-D1I

END)

SUBKgUTTNE CDEC0MP0SE(CiNT N»PCPnlNT A'CREAL DErR'CRtAL i>etI»

CI
BEGAN % ThE COMPLEX UNSy

where l is lower

MATRICES. AND Gv
AND ITS GENERAL
KEEpS THE RECOSO
A, THE OETERMIMA

C1NT I» j.K.L.P.PPj

PREAL VICTOR ATI6AJI

PKEAL Z7I

CKEAL X.Y.Z.V.W.H.HHI

LABEL FAlLl
MUDEI«TRUEl

MUDE»«TRUE AND pen lss
LUUP II»1.1»N DO

lNNpRoD(0»O.AtI]»AtI
DLTRI.ii DETII.Oj DETE
LUOP Kl-l.l.N 00
BLGIN

Lf«Kl Pl«K*K| PPI«P-

LOOP II«K»1.N DO
BEGIN

MoDEI«TRUEl

MoDE»«TRUE AND PEN

Innprodc-graboneca
InnpRod(-h,-hh.rtl
Innprodcgraroneu

lNNpROD(H.HH»ACIl.
Y|«-HJ

modei«true»
m0de1.true and pen
mode»«truej

NT DETE»CNP0INT INT) I

MMETRIC MATRIX A IS DECOMPOStD INJO LU.

triangular and u is unit upper triangular
erwritten on a, a is stored in N*2N array

ELEMENTS ARE At I*2j«l 1«I*A( I *2 J] t INTIM
Of ANy INTERCHANGES mAdE To THE ROWS Op
NT (0ETR*I#0ETI)*2**DETE IS ALSO COMPUTFD,

N + NJ

]ilNT(I].H)l
1*0 1

1) Zl»OJ

LSS N*N|
tI].PP"i)*0*A[I]»ATtPP].H*HH^
(l»*A(l]).ATtP].X,HH)|
tIJ.P-t),O.RTL(l»»ACI]).ATtPPj.H»HH)l

AT(Pl.H»HH)l

EQL PP"ll AC 1 3 f «X|

00000100?
00000200?
00000300?
OOOOOaOO?
00000500?
00000600?
00000700?
OOOOOBOO?
00000900?
00001000?
00001100?
00001200?
00001300?

oooouoo?

00001500?

oooouoo?
00001700?

00001800?
00001900?
00002000?
00002100?
00002200?
00002300?
00002400?
00002500?
00002600?
00002700?
00002800?
00002900?
0000300^?
OC003100?
00003200?
00003300?
00003400?
00003500?
00003f00?
00003700?
00003800?
00003900?
0000*000?
00004100?
00004200?
00004300?
00004400?
00004500?

00004600?
00004700?
00004S00?
00004900?
00005000?
00005100?
00005200?
00005300?
00005400?
00005500?
00005600?
00005700?

30

00006800?
00006900?
0000?OOP?

So°oDI!:TTRRUULAN° PEN E°L P'U AI»-TI 00005800,

HODE„TRUE AND PEN.II SSSSfSSSJ

l' * GTR Z ™EN 00006200?

ENB0riN Z,"X' L"" ENDl oZZlll

lFDL,SKUK,T!SSCUTm AN° ^ LSS "*"* OOOoJjSSJ

BEGIN DETHI-OFTRI XXSXJJaXI

22 j »ACK3 1
AfKll-Afi 1»

iNTtLi'^TNTrKii

E *WT|Uj lalNiCKjl 00007200?

lNTfKli«il 00007300?

x.-JrJboUcackj.pp-u, SSJSfJSS?

li:«™*«t*un» Z°o7rl°o°ol

OETj,;JjDETI*V*OETRl SSSS^SSS?

IF ABS(DETR) GTR ABS(OETI) SSSSl?So?

Then whdetr else hi-detii §§§§si§S?

JErG?N"L°ARTEHLNL?EGlN DETE'"01 G° T° FAIU EN°' 0S0O83OO?

Lli IF ABSCW, GEO 1 THEN 22SS2S222-

BEGIN W!.W*0.0625, 000860

°lVATrTi:°a\06'5t SSSS55SSJ

DETE.-DETEMI 00008800?

EN0| G° T° L1» 00008900?

EN0 00009000?

BEGIN LABEL L2i 00009100?

L2, IF IbjCO LSS 0.0625 THEN SSSSSSSS5

BEGIN S^ETR-H, = °0'

DETKI-0ETR*l6l 00009500?

0ETII.0ET1M6I 00009*00?

GO TO l\] 00009700?

end,

MnOEl.TRUE ANO PEN LSS K*K-2| 00010*00*

lNNPROO(-GRABONE(AtK).PP-n.O.AtKl.ATtPP],H»HH)l 00010600?

INNP«nO(-H.-HH.RTL<J..A[K]).AT[Pl,V.MH)l 00010700?

|NlJpR0D(-6RAHnNE(ACK].P-n,0.RTL(l,.A|RPl).AT[PPl,H.MH)l OOOIOSOO?

INNPR00(H.HHPA[K3.AT[P3.H.HH)| OOOloSooJ

SaScVimucf SSSInSS?

MROEI-TRUE ANO PEN-PP-U OOOllJooJ

Mnon'THun oooiuoo?

0000980O?
00009900?

31

modemtrue ANO »EN«P-1I 00011500?

AfK]|»<W*X-V*Y>/Zl 00011600?

END| 00011700?

ENO| 00011800?

FAIL! 00011900?

ENOl «0r CDECOMPOSE. 00012000?

00012100?
SUBHOUTINE CACCS0LVE(CINT N.CINT R»PCP0INT A.PCPOINT AA.CNPOINT P» 00012200?

PCPOINT B.PCPOINT X.PCPOINT BB.CiNT I) | 00012300?

BEqIN * SOLVES AX«B. WhERE A IS AN N*2N cOHPlEX UNSY"NETKIC ANq 00012400?

» b is an n*?r complex matrices* respectively. aa 1$ lu 00012500?

* decomposition proouceo by "cdecompose". the hesiuuals 00012600?
% bb«b-ax are also calculated and ad*88 is solved. over- 00012700?

* Writing d on bb. l is the number of iterations. 00012800?

00012900?

SUBROUTINE CS0LVE(CINT N.CINT R.PCPOINT A.CNPOINT iNT. 00013000?

PCPOINT B)| 00013100?

BtGIN 00013200?

CINT I,J,K,KK,P,PPl OOOI33OO?

PREAL VECTOR BT[6«]I 00013400?

CREaL X.Y»Z.21.Z2»H.HH» 00013500?

PREAL VECTOR Ct6«l| 00013600?

BOOLEAN AMOOEI 00013700?

LouP I««1»1.N On 00013800?

IF INTCII NEQ I THEN 00013900?

LOOP Jt«R*R."l»l 00 00014000?

begin oooiaioo?

Xi«QRABONE(B£IJ»J-l)l 000U200?

MODEI-TRUE AND PEN"J-ll 00014300?

B[Ilt»GRABONE<BClNT[U].J-l)l 0001 4400?

Bt INTC 111 l"XI 00014500?

EnOI 00014600?

LOOP K»«R*R»-2.2 DO 00014700?

BEGIN 00014800?

KkIbK-11 00014900?

LOOP Il'l.l.N 00 00015000?

BEGIN 00015100?

MODEl'TRUE ANO PEN LSS 1*1-2 J 00015200?

lNNPROO(-GRABONE(B(I]*KK-l).0»AtI]»BT(KK)»H.HH)| 00015300?

lNNPR0DCH,-HH.RTL(l..AHl).BTtK].X.HH)l 00015400?

INNPR0D(-GRAB0NECBII1»K-1)»0#RTL(1»»A[I]).BTIKK]»H»HH)J 00015500?

INNPR00(H.HH»A[I J.BTtKl.H.HH)) 00015600?

Y|«-HI 00015700?

P|«I«II PP|«P-i| 00015800?

Zi!»GRABONE( ACIJ.PP-DI 00015900?

Z2»»GRABoNF(AU J.p-1 >1 00016000?

ZiaZl*Zl«Z?*42l 00016100?

MoOEl«TRUE AND PEN.KKJ 00016200?

Btl Jl.(X#Zl*Y*Z2)/Zl 00016300?

M00EI"TRUE AND PENnKI 00016400?

Bf I]lB(V*Zl-X*Z2)/Zf 00016500?

END; 00016600?

LOOp Ii«N.-l,l DO 00016700?

BEGIN 00016800?

AMODE«« TRUE AND PEN GEO NJ 00016900?

M0DE«"AM0DEJ 00017000?

lNNpRGD<-GRABoNE<Btl)»KK-l).0»Atn»BT[KK]»H,HH)i 000171 00?

32

lNNPROD(-H,-HH.RTL(l»»A(l])»BTtK)»H,HH)l 00017200?

MoOEl»TRUE ANO PEN»KK-il C C I } I «HI 00017300?

MoOEl«AMOO£| 00017*00?

lNNpR00(*6RAB0NE(BCn»K"l)»0«RTi.(l**ACI])*BTUK]»H«HH)l 00017500?

INNPROO(H,HH»ACI]»BT[KJ,H.HH)I 000l7ft00?

MOOEl«TRuE AND PEN«K"ll C[I]H-HI 00017700?

M0DE|«TRUE AND <PEN«KK-1 OR PEN»K-1)I B[I)|nC[I]l 00017800?

ENOl 00017900?

EnDj 00018000?

ENDi % OF CSOLVE, 00018100?

00018200?

C1NT I»J»K#D2.C»CC» 00018300?

CKEAL E*00*01iXHAX'BBHAX*EPS*El*E2»H»HH| OOOUaOO?

PHEAL VECTOR XTC64 J j 00018500?

BUOtEAN AMODEl 00018600?

LaB£L L3»lLLl 00018700?

EKS|sl(0P«10l 00018800?

AMODE|«pEN GEO N| 00018900?

MUOEUAmOOEI 00019000?

CUO* I»»1»1«N Dq 00019100?

BtLGlN XtI]l>O.OI 00019200?

BB< I 1 l"BC Hi 00019300?

ENDj 00019400?

L»«0» Doinfll 00019500?

MUDEI-TRUEl 00019600?

L3I CJ»OLVE(N#R.AAiP.BB)l 00019700?

L««L*1» D2t«0l D1I»0I 00019800?

MoDEUAMODE» 00019900?

LUQP H«1#1»N 00 00020000?

XtllUXtU + BBCllI 00020100?

LUOP jl.l.l.R 00 00020200?

BLGiN 00020300?

XMAXI.Ol BBMAXl«0l 00020400?

C»«J*J» CCl-C-ll 00020500?

LOOp Il«l#l.N 00 00020600?

BEGIN 00020700?

El I.GRAB0NE(X[IJ#CC-1)I E2 | ■GRABONEC X 1 1 ] # C-l ) I 00020800?

Ei"El*El*E2*E2l 00020900?

Ip E GTR XMAX THEN XMAX|»E| 00021000?

El «»GRftB0N£(BBtl J»CC)I t2 I "GR aBONE ( BB 1 1 ) » C ) 1 00021100?

Ei«El*El*E2*E2* 00021200?

IF E GTR BBMAX THEN BBHAXlaEI 00021300?

0002U00?

00021500?

00021600?

00021700?

MnOEl.TRUEi 00021800?

lNNpR0D(-GRAB0NE(B[I]»C)»0.RTL(l».ACl J)»XT[C]»M»HHJ| 00021900?

lNNpRODCH.HH.ACll.xT(cl»H»MH)l nn07?nnn«>

MnoE««TwuE and pen»c-u bbi u i«-hi

in .

IF BBHAx/XMAx C,T« 01 THEN 01 | .BBHAX/XM Ax I
IF BBMAX GTR (2*EPS)*C2*EPS)*XHAX THEN D2l»H
101

Ip E GTR XMAX THEN XMAX|»E| 00021000?

El »»GRABONE(BBtl J»CC)I fc2 » -GR aBONE ( BB t I ) » C ) 1 00021100?

E»»El*El*E2*E2* 00021200?

IF E GTR BBMAX THEN BBHAXlaEJ

MnDE'"TRUEl

lNNpR0O(-GRABnNE(Btn.CC)»O,A[n.xTtCCl.H,HH)|

lNNpROOCH»-HH»RTL(l»»AC I 1 ) • XTt C 1 # H. HH) I

MODEI"TRUE AND PEN«CC"1» bBHll^HI
MnDri«TRU£i

00022000?

00022100?

00022200?

00022300?

00022000?
E^OJ 00022500?

It 01 GTR 00*0.25 AND L NEO 1 THEN GO TO ILLI 00022600?

OUi-Op 00022700?

IF qZ«1 THEN Go TO l3| 00022800?

33

ILLI 00022900?

ENDI X 0F CaCCSqLVc. 00023000?

00023100?

PREAL VECTOR AAiBB[6<»}» 00023200?

creal vector iNTcHt6*)i 00023300?

CINT I»N»R»OeTE«ITi 00023400?

CREAL DETR.OETII 00023500?

LOOP Ilal.l.N 00 00023600?

AAC I ] l»AC III 00023700?

COECoMP0Se(n»AA»OEtR»dEtI»OETE»InTCH)I 00023800?

CACCSOLvE(NtR»A*AA,lNTCH*B»X«BB*lT)i 00023900?

END| * OF CLlNSYS. 00024000?

00024100?

00024200?

SUBROUTINE INVITERATIONCCINT N»PCP0INT W.pCPOINT X.CINT KK)| 00024300?

BEGIN X ThE SUBROUTINE COMPUTES THE EIGENVECTORS AND PKlNClPAL 00024400?

X VECTORS BY SOLVING ( ( W-LAMBDA*! )**K )*X»0, THE SUBROUTINES 00024500?

I "cDECOMPOSE" AND "CACCSQLyE" IN "cLINSyS" ARE USED* 00024600?

PE hEaL Su8ROuTlNE sQRT As r6A(Pe rEAL ARQ AS r6A)I 00024700?

BEGIN 00024a00?

S* CALL SqRt64()| 00024900?

ENDI % QF SORT. 00025000?

CINT I»J»II»DETE.R»L» 00025100?

CKEAL DETRtDETI.MAXVAL.XXl,XX2| 00025200?

pkeal vector w2» b» b2» bb» xxt641 , 00025300?

ckeal vector p.rootcmji 00025400?

PHeAL Tri»TR2» 00025500?

LABEL INVIT»NEXT.FINISI 00025600?

MODEi«PEN LSS N*N| 00025700?

LOOP I>»1'1»N 00 W2C X 1 1 «•*£ III 00025800?

CDEC0MP0SE{N,W2.DETR#DETI#DETE.P)I 00025900?

M0DE1.PEN LSS 21 00026000?

LOOP ll«l*l#N DO 00026100?

BEGIN B2(I]lal,0l BC 1 1 l«B2C III ENDI 00026200?

KKl.ll 00026300?

Rl»ll 00026^00?

INVITI 00026500?

MOOE»»PEN LSS N*NI 00026600?

CACcS0LVE(N»R»m»m2»P»B*X*BB»L)I 00026700?

HODe««PeN LSS 21 00026800?

LOOP !I»1»1*N DO 00026900?

BEGIN 00027000?

RnOT[I)l"ROHSUH(XIIJ*X[Il)| 00027100?

ROOT t I ]I«SQRT(ROOTII ])l 00027200?

ENOi 00027300?

MAXvALI«ROOT[lll 00027400?

LOOP lt«2»l*N DO 00027500?

IT m*XVAL GEO HOOTtI] THEN MAXVAL I »MAXVAL ELSE MAXVal l»ROOT[ Ij I 00027600?

LOOP lialiUN DO 00027700?

IF ROOTtlj.MAXVAL THEN BEGIN IllaH GO TO NEXTl ENU| 00027800?

MEXTf 00027900?

XXl|»GRABONE(XUn»0)l 00028000?

xX2l«-6RAB0NE(XtII]»l)l 00028100?

IF -XX2 NEQ 0 THEN 00028200?

BEGIN 00028300?

LQOp I i * 1 • 1 • N DO 00028400?

BEGIN 00028500?

3^

MODEl«PEN LSS 21 00028600?

TRll»GRABnNE(Xtl].0)*XXl-GRABoNC(Xtl].l)*XX2| 00028700?

TR2l«GRAB0NE(X[n,0)*XX2*QRAB0NE(XCI].l)*XAil 00028800?

M00El«PEN«0l XXCIll.TRli 00028900?

M0DEl«PEN«ll XX(I}|»TR2| 00029000?

END| 00029100?

MQDEI.PEN LSS 21 00029200?

Xxl|«GRABONE(XXtII).0>| 00029300?

LOOP IIP1.1.N 00 XUJlOXCll/XXU 00029400?

END 00029*00?

ELSE LOOP Il.l.l»N 00 X[ IJ I.XC H/XXl | 00029600?

LOOP II-l.l.N DO 00029700?

ErlN -. «, . 00029800?

IP ABSCB2[I].X£I]> GTR P.lO THEN 00029900?

WE9IN 00030000?

IP KK GTR 10 THEN GO TO FlNISl 00030100?

loop ji.i.i.n oo begin Btnt.xtn» 00030200?

B2tI1l"X[I]l 00030300?

Kkl.KK*it END' 00030400?

*K»"KK*1I 00030500?

G0 TO INVITI 00030600?

If?r rn Tn r»m<;. 00030700?

LLSE GO TO FINISI 00030800?

„„,., H°l 00030900?

*** 0O031OOD?

CNOJ UF INVITERATION. 00031100?

UNCLASSIFIED

Security Classification

DOCUMENT CONTROL DATA • R & D

(SacuHty claaallleatlor, at till: body ol abahmet and humming annotation muat ba antarad whan tha ormr.lt raport t. claaalllad)

NAT.NG ACTIVITY fCo^Or.,. author) »*• ""OUT .ECU R. TV C L ASS. FIC A TION

originating activity (Corporata author)

Center for Advanced Computation
University of Illinois at Urbana-Champaign
Urbana, Illinois 6l801

UNCLASSIFIED

2b. GROUP

3. REPORT TITLE

A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS

4. DESCRIPTIVE NOTES (Type ot report and htelualra datam)

Research Report

8. AUTHOR(S) (Flrat mm, middla Initial, la a ( nam*)

Killion Noh

8. REPORT DATE

January 20, 1972

•a. CONTRACT OR GRANT NO.

DAHCOU 72-C-OOOl

b. PROJECT NO.

ARPA Order 1899

fm. TOTAL NO. OF PAGES

.22.

17b. NO. Or REFS

•a. ORIGINATOR'S REPORT NUMBERIS)

CAC Document No. 2k

Sb. OTHER REPORT NO<3» (Any other number* that may ba aemlgnad
thla raport)

UIUCDCS-R-72-U98

10. DISTRIBUTION STATEMENT

Copies may "be obtained from the address given in (l) above.
Distribution unlimited; approved for public release.

II. SUPPLEMENTARY NOTES

None

12. SPONSORING MILITARY ACTIVITY

U.S. Army Research Office-Durham
Durham, North Carolina

13. ABSTRACT

The eigenvector solution of the time -invariant matrix Riccati equation
is discussed. The coefficient matrix of the canonical equation is allowed
to have multiple eigenvalues, namely, the matrix could be either derogatory or
defective. The solution of matrix Riccati equation is then calculated from a
part of similarity transformation which should reduce the coefficient matrix
to the Jordan canonical form.

DD

FORM

,1473

UNCLASSIFIED

Security Classification

UNCLASSIFIED

Security Classification

KEY WORDS

LINK A

ROLE WT ROLE WT

ROLE WT

Ordinary and Partial Differential Equations

Matrix Algebra

Nonlinear and Functional Equations

UNCLASSIFIED

Security Classification

BIBLIOGRAPHIC DATA
SHEET

1. Report No.

UIUCDCS-R-72-1+98

3. Recipient's Accession No.

4. Title and Subtitle

A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS

5. Report Date

January 20, 1972

6.

7. Author(s)

Killion Noh

8. Performing Organization Rept.
No.

9. Performing Organization Name and Address

Center for Advanced Computation
University of Illinois at Urbana-Champaign
Urbana, Illinois 6l801

10. Project/Task/Work Unit No.

11. Contract /Grant No.

DAHCOU 72-C-OOOl

12. Sponsoring Organization Name and Address

U.S. Army Research Office-Durham

Duke Station

Durham, North Carolina

13. Type of Report & Period
Covered

Research

14.

15. Supplementary Notes

None

16. Abstracts

The eigenvector solution of the time -invariant matrix Riccati equation
is discussed. The coefficient matrix of the canonical equation is allowed to
have multiple eigenvalues, namely, the matrix could be either derogatory or
defective. The solution of matrix Riccati equation is then calculated from a part
of similarity transformation which should reduce the coefficient matrix to the
Jordan canonical form.

17. Key Words and Document Analysis. 17o. Descriptors

Ordinary and Partial Differential Equations

Matrix Algebra

Nonlinear and Functional Equations

17b. Identifiers/Open-Ended Terms

17c. COSATI Field/Group

18. Availability Statement Copies may be obtained f rom the

address in (9) above.
Distribution unlimited.

19. Security Class (This
Report)

UNCLASSIFIED

20. Security Class (This

Page
UNCLASSIFIED

21. No. of Pages

38

22. Price

NC

FORM NTIS-35 (10-70)

USCOMM-DC 40329-P71