i HP SBR nSH LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN ULcot GOO-2. ^fi^ZT A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS By Killion Noh January 20, 1972 CAC Document No. 2h DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Digitized by the Internet Archive in 2013 http://archive.org/details/numericalsolutio498nohk CAC Document No. 2k UIUCDCS-R-72-U98 A MJMERICAL 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. DAHCOU 72-C-OOOl. 11 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 ACKNOWLEDGEMENT 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, GLYENIR. 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(t Q )= x Q 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, = |y T (O Fy(tJ 4 f (y T (t) Qy(t) + u T (t )Ru(t ))at tQ (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 l = 2 xT(t f } cTpCx ( t f ) + \ 1 (x T (t)C T QC x(t) + u T (t)Ru(t))dt " t (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 _1 B T K(t)x(t) - (2.10 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 + A T K(t) - K(t)BR _1 B T K(t) + C T QC = (2.5) with the boundary condition K(t f ) = C T PC (2.6) Furthermore, the optimal trajectory is the solution of the system of differential equations x(t) = (A - BR~ 1 B T K(t)) x(t) (2.7) x(t Q ) = x Q The minimum cost is given by V 1 * = |x T (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 = and t = °°, the output regulator problem turns out to be the following: minimize V = | /~(y T (t)Qy(t) + u T (t)Ru(t)) dt (2.9) subject to x(t) = Ax(t) + Bu(t) y(t) = Cx(t) (2.10) x(0) = x Q Then the optimal control is given by u ( t ) = -R _1 B T K 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 + A T K - KBR _1 B T K + C T QC = (2.12) The optimal trajectory is the solution of the system i(t) = Gx(t) -1 T G = A - BR B K (2.13) :(0) ^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.11+) 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) = (2(t) - K(t) X(t)) X _1 (t) (2.16) into (2.5) and setting X(t) = AX(t) - BR _1 B T Z(t) (2.17) we obtain Z(t) = -C T QCX(t) - A T Z(t) (2.18) Therefore (2.5) becomes a system of 2n-dimensional linear homogeneous differential equations "x(t) i r _Z(t) J [. -1 T -BR B T T -C QC -A X(t) Z(t) (2.19) Since we are interested in increasing the time t, let t = t f - t (2.20) Then X(t) Z(t) = ¥ X(t) Z(t) (2.21 where W = -A T C QC -1 T BR B with X(0) = I Z(0) = C X QC (2.22) Therefore, the solution will be given by the following matrix exponential X(t) Z(t) = e Wx T C QC (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(x) are obtained from the similarity transformation which reduces W into the Jordan canonical form. Throughout the remaining discussions, the increasing time variable t will be used instead of x. 3. SOLUTION OF THE MATRIX RICCATI EQUATION We begin with two lemmas for our main result. Lemma 1. The eigenvalues of the matrix W of (2.21) are symmetric with respect to the imaginary axis. Proof. Let I be a 2n x 2n matrix given by x i = -I where I is an n x n identity matrix. Then it is obvious that I - 1 = I T 1 1 T T I WI = -W (3.1) ("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 = A ll A 12 A A L A 21 A 22J (3.M where A , A , A , and A are n x n, m x m, n x m, and m x n matrices, respectively. Assume that A^ is non-singular. Then „-l T T 11 12 T T L 21 22J (3.5) exists iff A - A A A.^ is non-singular, and is given by T__ = A,," 1 - -. . u A „ T. "11 11 1 ^2 T 22 A 21 A ll ' -1 T = -A "" A T 12 A ll A 12 X 22 3.6) T 21 " T 22 A 21 ^1 -1 T = (A - A A 22 VA 22 A 21 11 -1 *u.» -1 where submat rices T^, T^, T ±2 , and T 21 are the same sizes as A A 22' A 12' and A 21' res P ectivel y- 11 -1 -1 Proof. Suppose A g2 - A gl A^ A L2 is non-singular. From TT = I, A ll T ll + A 12 T 21 ' : „ A ll T 12 + ^2 T 22 = ° A 21 T ll + A 22 T 21 " ° (3.7) A 21 T 12 + A 22 T 22 = Im -1 Premultiplying the second equation by A A " and subtracting from the fourth, T 22 (A 22 " A 21 A ll '" A 12 } ' And hence -1 T 12 ^1 h2 T 22 Similarly, from the first and third, T 21 = " T 22 A 21 ^Ll T u ■ A n _1 - A ll * v "I m "I ^2 T 22 A 21 A ll 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 -1 . ? det ^.1 A 12 A A A 21 A 22J = det \i 0l dl m ] F : 12 ° A 22" A 21 \l -1 12 J -1 " det < A ll' det(A 22 " A 21 *U " A 12 ) Therefore, A - A A A is nonsingular. (3.8) Assuming that W has distinct eigenvalues, we now turn to constructing the similarity transformation with which ¥ of (2.2l) is reduced to the 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 W, such that the first n columns are the right eigenvectors of W corresponding to the n eigenvalues of W with positive real parts, then rx„ S _1 ¥S = X. (3.9) 2n 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 = S ll S 12 . S 21 S 22J (3.10) then the solution of the matrix Riccati equation is given by K - S 21 hl^ 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.1l) holds for the case in which W 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 = *! (3.12) m, and in which m. A. 1 1 A. 1 l \ A. l , (m. x m. ) (3.13) If we denote S by its column vectors s, . .... s^ , then 1 2n ¥ Ls l5 s 2 , ..., s 2n ] = [s 1 , s 2 , ..., s 2n ] J (3.1*0 For any Jordan block J , if m. > 1 then, dropping subscripts of s for i convenience. Ws (l) =s (1) A ws = s + A . s 1 (3.15) or Ws (m. ) . (m.-l) (m. ) 1=S1 +A.S1 (¥ - A.I) s (l) = l (¥ - A.I) 2 s (2) = (3.16) / sm. (m. ) (W-A.Ijis i =0 Thus, s , 1 < j < m. , is a principal vector of degree j associated with (l ) (m ) A.. Since S is nonsingular, the m. principal vectors s , ..., s i 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 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) intc J = (3.17) where J and J are diagonal with the eigenvalues of W with positive and negative real parts respectively, and E and E p are matrices with one's on the super-diagonals and zeros elsewhere. To associate with (3.6), let us partition S and S as s S „ 11 12 s = s s„„ 21 22 ^ T T -l 11 12 s = T T 21 22 (3.18) where all submatrices are n x n. Note that, by lemma 2, it is necessary to assume that S and S 22 - S S ~ S are nonsingular. Then Wt SJS -1 t e = e - Se^s" 1 s ll S 12 f S 21 S 22 Jit n e 1 e <= At E 2 t T T 11 12 T T 21 22 S 11 e J l t e E l t T 11 + Sl2e W S \ ^eWl^ + ^e^ J.t E,t J«t E„t J. t E n t r S-_e 1 e 1 T + S./2V2X, S_e u lWT + R„e d 2V2 l T 21 J^t E„t r "11 22 21 21 12 22 "22 (3.19) 11 Since the elements of J have negative real parts, e 2 -*• as t -> °° and thus the terms having e 2 in (3.19) vanish as t -* °° . Therefore, X(t) Z(t) = e Wt T C QC ri T C QC S 1:L eV eV (T 12 + T^ C T QC)' S 21 e* 7 !* eV (T + T C T QC ) (3.20) Therefore the solution of the matrix Riccati equation is given by K = lim Z(t) X _1 (t) = lim [S eV eV (T-.+T C T QC) ] [S eVeV (T +T C T QC) I" 1 ■ S 21 S ll -1 (3-21) provided that the indicated inverse exists. However, we shall show that only the nonsingularity of S and S - S S S ? 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 S ll S 12 S 21 S 22 = J (3.22) Equating the (2,1 ) elements on each side of (3.22), (-T 21 A + T 22 C T QC) S 1X + (T 21 BR _1 B T + T^ A T ) S 21 = 12 or " T 21 AS 11 + T 22 ^ S 21 + T 21 BR_1bT S 21 + T 22 ^^ S ll = ° (3 ' 23) Substituting T = - T 22 S 21 S^" into (3.23), T 22 S 21 S ll" 1+ T 22 ATS 21 " T 22 S 21 S ir lBR " lBTs 21 + T 2 2 cT «° S ll=° (3 ^ Premultiplying by T and postmulti plying by S , S 21 S 11 _1a + aTs 21 S 11 _1 " S 21 S 1i" 1bR " 1bTS 21 S 11 _1 + ^ = ° (3 ' 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 W 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, lM and the Jacobi-like algorithms [15, l6] are already available. 14 Test Result 1, A = 1 1 -1 -1 1 -1 -1 p = B = 1 1 1 Q = 0' 10 10 .0 c = 1 1 1 1 1 R = 10 10 1 The numerical example was tested "by the QR algorithm and the inverse iteration method. The ten distinct eigenvalues of W are ± 1, ± 1.7287601*77185 ± 1.577533767573i, ± 1.35319578U0U6 ± l.l537^989928Ui. 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. i+9^009759 -0.819173651 0.668267901 -O.UU3608958 2.1+9^009759 7.1+351+51161+ -1.8257^1858 1.122910U32 -0.668267901 -0.819173651 -I.8257U1858 1.6383U7303 1.8257^1858 -0.819173651 0.668267901 1.1229101+32 1.8257^1858 7. 1+351+ 511 6k -2.1+91+009759 -0.1+1+ 3608958 -0.668267901 -O.819173651 -2.1+91+009759 I.262782609 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 = P=Q= The mat rix 3 -2 2 -1 1 W = -3 -2 2 1 c = R = 1 1 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.000008*+ 6.000001+3 6.0000031+ -k. 000002T 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 -1+ 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 h 1+ 1 1+ 6 1 k 1+ 1 6 k 1 1+ k 6 B = 1 C=R= 1 1 1 1 P=Q= 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 17 Ik. 99999999930 ^.9999999991+9 -h. 99999999982 0.99999999997 -1^.9999999981+ U. 99999999981 j+. 999999999)48 -0.99999999995 The computed solution K is given "by ■1.999999999811 1.99999999988U 2.000000000276 -2.00000000036U 1.9999999999U2 1.999999999913 -2.000000000015 -2.000000000000 -2.0000000001+07 -2.000000000378 2.000000000509 2.0000000001*80 -2.000000000015 2.000000000087 2.000000000U80 -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^ (1963), pp. 6U3-6U5 [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. SIAM Appl. Math., ik (1966), pp. 1+96-501 [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., ik (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-23^ [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 [ik] 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 t TWO ALGot PROCEDURES' RlCCATU AnO RICCATI2. FOR COMMUTING THE OOOOOlOO? f SOLUTION OF THE MATRIX RICCATI EQUATION ARE LISTED BtLOW* FOR 00000200? ( «RICCATlU. THE RROCtOURE "EIGEN* (BY EBERLEIN AND BUQThNoyD) WAS 00000300? I USEU WHILE A COMBINATION OF "HQR" (BY MARTIN. ET AL) AND 00000400? I MNVITErATION" WAS APPLIED TO "RICCATI2* f LINEAR EqUATIUn SOLVER 00000500? % (BY BOWDlER. E T AL) WAS ALSO USED TO PROVIDE THE INytRSE OF A MATRIX 00000600? « ANO TO CALCULATE THE PRINCIPAL VECTORS FROH ( ( W-LAMBUAM >**K>*X«0. 00000700? 00000800? 00000900? PRDCEUURL RlCCATIl(N.L»M»A»B#C»P»O.R»THX»W»K)l 00001000? VALUE N.L.M.TMXI 00001100? INTEGtR N,L.M.TMX| 00001200? ARRAY A.B.C.P»a«R»W»K[l»ll| 00001300? COMMENT SOLVES ThE MATRIX RlccATI EQUATION . OOOOUOO? KA*(AT)K-KR(RINV)(BT)K*(CT)QC«o, 00001500? A Is N*Nt R IS N*L» C IS M*N» P AND 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 OOOOUOO? Op w ARE (1»1> B *A» (1#2)«B(RINV)". REARRANGE EI GE N VE c TOK C Pr INC I P AL v EcTOR> H A TRlX 00002000? SO THAT ThE rlHST N COLUMNS CORRESPOND TO TO THE EIGENVALUES 00002100? WlT H P0SITI V E REA L pARTS. dENoTINq ThE UPP^R H a Lf Op N COLUMN 00002200 ? VECTORS BY SMTXl ANO THE LOWER HALF BY SMTX2. I HE SOLUTION 00002300? K IS THEN GlvE N BY K»C SMTx2) ( SMTxl INv ) . THE PRUCEDURE 00002a00? "plGEN" IS USED ^OR EIGEnPrUB L EM F Wi 0000250O? 00002600? PEGIN COMMENT plRST DECLARE PROCEOuRESl 00002700? 00002800? PROCEUURE SETUP(N.L»M»A»B»C»Q»RINV.W)| 00002900? VALUE N.L.M, INTEGER N»L»M> 00003000? REAL ARRAY A • B# C » Q . R I N V # W ( 1 » 1] I 00003100? COMMENT THE PROCEDURE SETS UP ThE 2N*2N M ATRIX W. rInV IS INVERSE Op RI00003200? BEGIN OOOO33OO? REAL ARRAY R I NVBTU I L» 1 I N] . QC [ 1 I M. 1 I N ] | 00003*00? REAL SUM) 00003500? INTtGER I.J.KI 00003600? FOR It"i STEP 1 UNTIL N DO 00003700? TOR Jl«l STEP 1 UNTIL N DO 00003800? 8LGIN 00003Q00? W(I.J1I.-A(I,J]| 00004000? W(N*I,N*J]l«AtJ.I]» 00004100? ENUl 00004200? FOR Il«l STEP 1 UNTIL L 00 00004300? FOR J»«l STEP 1 UNTIL N DO 00004*00? BEGIN StlHl.O.Ol 00004500? FOR K i-l STEP 1 UNTIL L DO 00004600? SuM|«SUM*RINV(I»K]*B[J»KJI 00004700? RINvBTd.JJl-SUMl 00004800? ENOl 00004900? TOR I|»l STEP 1 UNTIL N 00 00 ?Sc°2 07 FOR Jl-1 STEP 1 UNTIL N 00 00005100? BttilN SyMl.O.Ol 00005200? FOR Kill STEP 1 UNTIL L 00 22SS! 3 221 SuH|«SUM*B[I,K]*RlNVBTtK»J]| 0000*400? "CI.N*J]l«SUMl 00003500? END| 00005600? FOR If .1 STEP 1 UNTIL M DO 00005700? 22 FOR Jial STEP 1 UNTIL N 00 000058001 BtGlN SljMlaO.Oj 000059001 FOR K t «l STEP 1 UNTIL M 00 000060001 SuMI«SUM*QtI»KJ*C[K»J]l 000061001 OCC I • J] laSUMl 0000*2001 EN0| 000063001 FOR I|«l STEP I UNTIL N DO 00006*001 FOR J is 1 STEP i UNTIL N 00 000065001 BtGlN $(jMl«0,Ot 000066001 FOR K f «l STEP 1 UNTIL N 00 000067001 SijMI"SUM*CCK t lJ*ac[K»JJI 000068001 W[N*I,J)|«SUMI 000069001 ENOi 000070001 ENO OF SETUPj 000071001 000072001 COMMENT INSERT THE PROCEDURES INNPROD. DECOMPOSE. AND ACCSOLvEi 000073001 oooofaoo? PROCEUURE INVERSE(N»A.X)| 000075001 VALUE N| INTEGER N, 000076001 DOUBLE ARRAY A.XU.ljf 0000/7001 COMMENT TmE PROCEDURE COMPUTES THE INVERSE OF AN N*N RtAL 000078001 UNSyMMrTRIc MATRIX A, USINq The PROCeOUReS "DrC0MPU$ E « AND 0000?900l "ABSOLVE", WRITING ON X| 000080001 REGIN 000081001 INTtGER UJ.IT.D2l 000082001 REAL Oil 000083001 ARRAY AA.RB.BCllN.l |N]*PC1 INll 00008*001 FUR Il«l STEP 1 UNTIL N DO 000085001 FuR Jl«l STEP 1 UNTIL N Do 000086001 IF J«I THEN BtI.J]»«1.0 ELSE BCI»J]l«0.0> 000087001 FUR I,«i STEP 1 UNTIL N DO 000088001 FUR J i«l STEP 1 UNTIL N DO 000089001 AAC I.J] laA[I,J]| 000090001 DtC0Mp0sE(N»AA,0l.D2.P)l 000091001 AtCSULVE(N.N.A.AA.P.B.X.BB»IT)| 000092001 ENO INVERsEl 000093001 00009*001 comment Insert the procedure eigeni 000095001 000096001 PROCEUURE MuLT(Nl.N2,NJ»AiB.C)l 000097001 VALUE N1»n2.N3| iNTEfiEH Ni»N2.N3l 000098001 DOUBLt ARRAY A.B.Ctl.UI 000099001 REGIN 000100001 INTtGER I.J.KI 000101001 REAL SUM1 00010200? FOR It«l STEP 1 UNTIL Ni DO 00010300? FOR J i ■ 1 STEP 1 UNTIL N3 00 00010*00? BEGIN SuH|aO.Ol 00010500? FUR k,"1 STEP 1 UNTIL N2 00 00010600? SUMI»SUM*A[I»K)*BCK,JJ| 00010 700? C I I . J ] I.SUMI 00010800? END* 00010900? END MULTl 00011000? 000111007 ARRAY HiN V »SMTXl,SMTX2UlN.llN].wW.TEMP,VECT(llN*N.l!N*NJ| 00011200? INTLGER ArRAy SEARCHlHNll 00011300? iNTtGEH N2»I.J.InDEX» 00011*00? 23 N2I*N*N| OoOlHOO? INVtRS£(L.R.RlNV)l OOOilloo? SETUP"-a. (1.2).b(rinv)(bt)» ( 2» 1 >■< c i >0c» and 00015300? C?»2)»(AT). REARRANGE E1GEnVE C T0R(PrIN C IPAL VECTOR) MATRIX 00015400? SO THAT THE rlKST N COLUMNS CORRESPOND TO TO TM6 EIGENVALUES 00015500? «ITh POSITIVE HEAL PARTS. DENOTING THE UPPER H*LF "F N COLUMN 00015600? VECTORS BY SHTXi AND THE LOWER HALF BY SMTX2» THE SOLUTION 00015700? K Is THEN GIVEN BY K.( SMTX2) ( SMTX1 INV ) . THE PRUCEDURES 00015800? "mShLD"' "hOr"» a N0 "INVITEHATlON" ARE USED TOM EI<»ENPRORlEM 00015900? Op M | 00016000? 00016100? BEGIN COMMENT FIRST DECLARE PROCEDURES! 00016200? 00016300? PROCEDURE SE!TUP(N.L»M.A»B#C»a»RINV.W)| OOOI64OO? VALUE N.L.Mi INTEGER N»L.M» 00016500? RF.AL ARRAY A » B » C » Q . R I N V • W t 1 » 1 J I 00016600? COMMENT T H E PROCEDURE SETS UP ThE 2N*2N MATRIX W. RINV IS INVERSE Op Rl 00016700? PEGIN 00016800? REAL ARRAY R I N VBT t 1 I Li 1 I N J . QC t 1 I H . 1 I N J I 00016900? REAL SUmi 00017000? iNTtGER I.JfKI 00017100? 2k FOR II" for Ji« BtGIN W[I W[N ENOI FoR Il« FOR Ji« BtGIN FOR S RIN ENO| FOR I|» FOR Jl" BLGIN FOR S wc I ENOi FOR Iia FOR J|" BtQlN FOR S act enOi FOR Il» FOR J|« BLGIN FOR S NCN END| END OF SE 1 STCP 1 UNT 1 STEP 1 UNT , J] I«-At I ,J1 ♦I.N*J]I«A[J 1 STEP 1 UNT 1 STEP 1 UNT SUM|«0»0| K|»l STEp 1 (jM|«SUM*RINv VBTC I* J] l«SU 1 STEP 1 UNT 1 STEp 1 UNT SuMl«0,0| K|«l STEP 1 UM,«SUM*BtI. .N + J] |"SUm» 1 STEP 1 UNT 1 STEP 1 UNT Sum«"0«0| K|»l STEP 1 UM|»SUM*Q( I, I# jl laSUMl 1 STEP 1 UNT 1 STEP 1 UNT SUMl«0,Ol K |«1 STEP I UM|«SUM*CtK. ♦I, J] l"SUMl TUPI IL N DO IL N DO I • I Jl I*. I Do II N DO UNTIL L DO f lfK]*B[Jf K]l Ml IL N DO IL N 00 UNTIL L DO Kl*RINVBTCKf JJ| IL M 00 IL N DO UNTIL M DO K]*C[Kf J] I IL N DO IL N DO UNTIL M DO I]*UC[K»JJI COMMENT INSERT THE PROCEDURES INNpROD. DECOMPOSE. AND AtCSOLvEl PROCEDURE INVERSE(N.A.X)! VALUE N| INtEqEH N| DOUBLE ARRAY A.XU.lJI COMMENT ThE PROCEDURE COMPUTES THE INVERSE OF AN N*N RtA|_ uns y mmetri c matrix a. using the procedures "DECQmpUse" "ACCSOL v E". writing ON Xl REGIN INTLGER I.J. IT. 021 REAL ol| ARRAY AA.BB.BlliN.l t NJ.Ptl|Nli FUR Ii-i STEP 1 UNTIL N DO FuR J i "i Step i u n til n d IF J«I THEN B[I.J]I»1.0 else FUR I,«i STEP 1 UNTIL N DO FUR Ji-i STEP 1 UNTIL N DO AAlIf J] l«A[I, J]| DECOMPOSE (Nf AA.0l.D2. P)l ACCSOLVE(N.N.A.AA.P.B.X.BB.IT)i END INVERsEi AND B[ If J] i»0.0J 00017200? 000173001 00017400? 00017500? 00017600? 00017700? 000179007 000179007B 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 comment InSeRt the procedures hshld and hori PROCEDURE HuLT(Nl*N2.Ni*AtBtC)l VALUE N1»n2,N3i INTEQEK Nl.N2.N3f DOUBLE ARRAy A»B*C[1.1)| BEGIN INTtOER I.J»KI REAL SUM! FOR Ii«l STEP 1 UNTIL Nl DO FOR Jl«l STEp 1 UNTIL N3 DO BEGIN SllM|«0,OI FUR Ki«i STEP 1 UNTIL N2 DO SUM|«SUM*A[ I.K)*B[K. J] | C t I . J 1 I -SUM I ENOI ENO MULTl PROCEUURE lNVlTERATION(N»w»X»KK»P)| VALUE N, TNTEGER N.KK| w.xci.n, 'i DOUBLE ARRAy array pchi c mme n t s lveS wx u MXTrI*. wr T H E ElQENy N Is THE NEEoEO for BEGIN INTEGER I»J»D2»R DOUBLE OI.MAXVAL DOUBLE ARRAY MH( LABEL iNVlTiOUTl FOR Ii"i STEP 1 FOR J l - 1 STEP I WW[ i, j]|«rf[ I.J DECUMpOsE(N.WW.D FOR Il«i STEP 1 81 1.1 ] i.l .01 KK I "1 J R | ■ 1 | INVITI ACCSOLVf;(N.R. «» W MAXVAl.l-0 j FOR Ii"l STEP 1 If AB S (xCI.l]) BEGIN MAXVALI"AB MAXjNOEXl* ENDi MAXVA|_|.X[MAXI FOR I i » 1 STEP 1 XI I. 1 ] I.X[ I.l ] FOR It»i ST^P 1 l> ABS(BCI.I)- BtGIN if kk gtr io FOR J|«l STE BtJ.l]fX[ ■0 BY I N VERSE ITERATION* WHERE * IS A N N*N REAL ITING SOLUTION ON X, THUS. THE PROCEUURE COMPUTES ECTOHS CORRESPONDING TO THE GIVEN EltiENv*LUES. RDER UF N AND KK IS THE NUMBER OF iTtRATiONS ITERATING EACH ElGENvEcTORl .LiMAXlNOEXi I |N.l|NJ.B.BBll|N.l |U| UNTIL N 00 UNTIL N DO ]l l.D2»P)» UNTIL N DO W.P.B»X.BB.L)J UNTIL N DO GTR MAXVAL THEN s(xri*n)i ii N0EX»1]I UNTIL N DO /MAXVALI UNTIL N 00 X[I.i]> GTR p-12 THEN then go to outi p 1 until n du JM 1» 00022900? 00023000? 00023100? 00023200? 00023300? 00023400? 00023500? 00023600? 00023700? 00023800? 00023900? 00024000? 00024100? 00024200? 00024300? 00024400? 00024500? 00024600? 00024700? 00024800? 00024900? 00025000? 00025100? 00025200? 00025300? 00025400? 00025500? 00025600? 00025700? 00025600? 00025900? 00026000? 00026100? 00026200? 00026300? 00026400? 00026500? 00026600? 00026700? 00026800? 00026900? 00027000? 00027100? 00027200? 00027300? 00027400? 00027500? 00027*00? 00027700? 00027800? 00027900? 00028000? 00028100? 00028200? 00028300? 00028400? 00028500? 26 KK I «KK* 1 1 GO TO INVITI ENU| OuTl END INVREsElTi INTLGER ARRAY Te EI REAL PER INTtGER N2l"N*N| INVtRSE( S E TUP DEFECT[JJ] THEN N R 1 1*1 STEP 1 UNTIL N2 00 H[I»I]i«W(i,n*pERTURB*MpLIERI LT(N2»N2,N2»W»W,TEMP)I R Ii.l STEP I UNTIL N2 DO R Jlal STEP 1 UNTIL N2 DO *l I»J]««TEMP[ I»J)I DEROGATEtJJI THEN Il«l STEP 1 UNTIL N2 00 I.I]t>MCNl)^PERTURB*JJ| RATION(N2.W.VECTOR.ITCOUNT»INTCH)I ■1 STEP i UNTIL N2 DO ■1.2 DO tl»Jj3l«VEcT0RCl»l), Search n columns of smtx corresponoiNg to n lige n values WITH POSITIVE REAL PARTS. SMTxl AND SMTx2 ARt ITS UPPER and lower halves* respectivelyi I STEP 1 UNTIL N2 DO RtU GTR THEN X,«INDEX+1, ChCINOEXJIbII 00028600? 00028700? 00028800? 00028900? 00029000? 00029100? 00029200? 00029300? 00029400? 00029500? 000296O0? 00029700? 00029ft00? 00029900? 00030000? 00030100? 00030200? OOO3O3OO? 00030400? 00030500? 00030600? 00030700? 00030800? 00030900? 00031000? 0003UOO? 00031200? 00031300? 00031400? 00031500? 00031600? 00031700? 00031000? 00031900? 00032000? 00032100? 00032200? 00032300? 00032400? 00032500? 00032600? 00032700? 00032800? 00032900? 00033000? 00033100? 00033200? 00033300? 00033400? 00033500? 00033600? 00033700? 00033*00? 00033900? 00034000? 00034100? 00034200? 27 FOR J I «i STEP 1 UNTIL N 00 FOR Ii«l STEP 1 UNTIL N DO BEGIN SMTXltl , J) l-SMTXf I.SE, SliTX2(I ,JJ|»SMTX(N*I •SEARCHIJIJJ endi INVtR$E(N .SMTX1, >TEMP)I MULT(N*N» N»SMTX2#TEMP» K), END RlCC A Tl2» 00034300? 00034400? 00034500? 00034600? 00034700? 00034800? 00034900? 00035000? 00035100? 28 APPENDIX II ILLIAC IV PROGRAMS 29 * two iumc iv glypnir programs. »clinsys" ano "invitation" are I LISTED BELOW, "CLlNSYS" IS USEO TO OBTAIN THE INVERS* OF A MATRIX x ano the solution or * complex system of unsymmetric linear eouations. % "INVITERATION" is useo to calculate the principal vectors of degree % K FKOM ( ENDl END) M00EI,TRUE| MOOEI.TRUE AND PEN LSS N+Nl IF L NEO K THEN BEGIN OETRI--OFTRJ OETII.-DFTU 22 a »At K 5 1 AfK] I-aC. H AtLJI»ZZ| INT[L] l«!NTcKj| END| INT[K]|«LI Xl» G RABONE(AtK],pp-l)| Yl«GRABONE(A[K],P-i)| Z|»X*X*Y*Y| Ht«X*DETR-Y*DETI» DETf t»X*DETI*Y*DETRl DETrI b N| IT ABS(DETR) GTR ABS(DETI) TmEn wi»oetr else WI«dETII IF w .0 THEN BEGIN DETE|«0| BEGIN LAREL Lll IF ABS(W) GEQ 1 THEN BEGIN W|.W*0. 06251 DETll"DETI*0t0625l DETEI»DETE*«» GO TO Lll END| ENDl BEGIN LABEL L2l IF ABs(W) LSS 0,0625 BEGIN WI.WM6I DETRl.DETR*16l DETII»DETI*U> DETEl-nETE**l GO TO L?» ENOi ENDl LOOP BEGIN P|»J*J| PP|«p. MnDEUTRUEl MflOEl.TRUE AND PEN lNNpROD( GO TO FAILI ENOl THEN JI«K*l.l»N DO II LSS K*K-2| iNNpR0D(-GRA 8 0NE(AtKJ.PP-l).0.At K ],ATtPPl.H.HHJl I N NPRnD(.H.-HH.R T L(i.. A CK]), AT [p;,v,Hi), NlJpRbD(-6RABONE(ACKl.P-n;oJlTL 1. AtPPn.ATrPPi M uun iNNpROOCH.HH.AtKJ.ATCPJ.H.HH)^ C ] ' HH) ' W|«-Hl MoDEI«TRUEl MnOEl"TRUE AND PEN«PP A[K],»(V*X*W*YJ/Z| MnOEl-TRUEl II 00005800? 00005900? 0000^000? 00006100? 0000620G? 00006300? 00006400? 000065C0? 0000660C? 0000*700? 00006800? 00006900? 0000^000? 00007100? 00007200? 00007300? 00007400? 00007500? 00007600? 00007700? 00007800? 00007900? 00008000? 00008100? 00008200? OOOO83OO? OOOO84OO? 00008500? 00008600? 00008700? 00008800? 00008900? 00009000? 00009100? 00009200? 00009300? 00009400? 00009500? 00009600? 00009700? 0000980O? 00009900? 00010000? 00010100? 00010200? 00010300? 00010400? 00010500? 00010600? 00010700? 00010800? 00010900? 00011000? O0O1H00? 00011200? 00011300? 00011400? 31 MoOtl«TRUE AND PEN»P-1» AfK)l»(K*X-V*Y)/ZI END| ENO| FAILI END) «Or CDECQMPOSE. SUBKOUTINE CACCSOLVECCINT N.CINT R»PCPOINT A.PCPOINT AA.tsiPOlNT PCPOINT B»PCP0INT X#PCPOINT B8»ClNT L) % BEfiiN * SOLVES AX«B» *hErE A IS AN N*2N COMPLEX UNSY"METHI c AN[) * B IS AN N*2R COMPLEX MATRICES* RESPECTIVELY, AA J$ LU * DECOMPOSITION PROOUCEO BY "CDECOMPOSE". THE KESIUUALS % B8"B-Ax ARE ALSO CALCULATED AND AD«B8 IS SOLVED. OVER- * WrITINg ON BB. i IS THE Num»Er Qf ITERATIONS. SUBROUTINE CS0LVE(CINT N.ClNT R»PCPOlNT A.CNPOlNT PCPOINT B)| BLGIN CINT PREA CREA PREA bool Loup If INT , I. J.K.KK.P.PPl L VECTOR BTt 6 «J» L X.Y#Z»ZUZ2»H.HH| L VECTOR Ct6«l| EAN AMODEI I i«l »i*n Dn iNTtll NEQ I THEN LOOP Jt»R*R.*l»l DO begin Xl»GRABONE(BCI],J-l)J MODEI-TRUE ANO PEN"J-1| Bt 13 i«GRABONE(B( INTC11]»J-1)I BClNTtUll'XI endi LOOP K««R*R»-2.2 DO BEGIN KkI.K-1 I LOOP I IM • 1 »N DO BEGIN MODEl-TRUE AND PEN LSS I* 1-2 1 lNNPROD(-GRABONE(B(I)*KK-I)»0*AtI)»BTCKK]»H»HH)| lNNPR0DCH,-HH»RTL(l#tAtl]).BTCK]»X.HH)l INNPROD(-GRA80NE(BCI]»K-1)»0»RTL(1»»A[I1)»BTIKK]'H#HH)I lNNPROO(H«HH*A[I)*BT[K]*HtHH)l Y|«-HJ P|»UII PP|«P-ll Zll«GRABONE( At I)»PP-1)I Z2I«GRAB0NE(ACI)»P-1 )l Zi«Zl*Zl*Z2*Z2l MOOEI«TRUE AND PEN.KKI BCI Jl.(X*Zl*Y*Z2)/Zl MODEI-TRUE AND PEN»KI B[I]I«(V#Z1-X*Z2)/ZI END; LOOP Il"N.-l,l DO BEGIN AMODEI* TRUE AND PEN GEO N| MnDE«"AMODE» lNN»»RGD(-GRABONE(Bt I J.KK-p.OtAt Il»BTtKKJ#H.HH)| 00011500? 00011600? 00011700? 00011800? 00011900? 00012000? 00012100? 00012200? 00012300? 00012400? 00012500? 00012600? 00012700? 00012800? 00012900? 00013000? 00013100? 00013200? OOOI33OO? 00013<»00? 00013500? 00013600? 00013700? 00013800? 00013900? 00014000? 00014100? 000U200? 00014300? 00014400? OO0U50O? 00014600? 00014700? 00014800? 00014900? 00015000? 00015100? 00015200? 00015300? 00015400? 00015500? 00015600? 00015700? 00015800? 00015900? 00016000? 00016100? 00016200? 00016300? 00016400? 00016500? 00016600? 00016700? 00016800? 00016900? 00017000? 00017100? 32 lNNPROO(»H,«HH.RTL(l»»A(!l>fBTtK]»H,HH)f MoDEl'TRUE AND PEN»KK«1I C C I J I »H I MQDCl«AHOOei lNNpROO( a GRABONE(B[I]*K-l)»0*RTL(l»»ACI])(BTlKK]*H'HH)l INNPR00(H.HH»ACI],BTCK] ( H.HH}I MOOEl«TRuE AND PEN«K-ll C[!]I«"HI M0DE|"TRUE AND (PEN»KK-1 OR PEN«K-1)I B[I]|«CCI]J END! ENOj ENDi t OF CSOLVE, CJNT I»J»K#D2»C» CKEAL E.DO.Ol.xM PHEAL VECTOR XTt 9U0LCAN AHOOEl L A B E L L3.1LLI EKSialtOP«lO) AMODE|«pEN GEQ N MUDEUAmOOEI LuQP I »*l • I »N Do BtGlN Xf 1 1 t «0 • I BBC 1 3 •■BC I END| L»«0» Do»*0» MUDEI.TRUEI L3l CJ>OLVE(N»R»AA.P» L««L*1» 02t«0l MnDE I.AMODEI LUOP I I ■ 1 # 1 » N DO X[ IJ«»X[ H*BB[ LUOp jlalflfR 00 BtGjN XMAXl.Ol BBHA C|«J*J| CClaC" LOOp paliltN BEGIN El I.GRABONEC Ei»El*El*E2* Ip E GTR XMA El J-GRaBONeC Et"El*El*E2* IF E GTR BBM MflDE«*TRUEl lNNpROO(-GRA iNNpROOCH.- MODE«»TRUE A M n Dj;l«T R UEI iNNpRQDCGRA lNNpRQD(HiHH M f1DE««THUE A END 1 IE BBHAx/XMAx IF BBMAX GTR ( ENOJ IF Dl GTR 00*0.2 DU , «D 1 I IF o2«l THEN Go CC» AX»BBNAX»EPSfEl«E2«H»HH| 6 4 J, ll ll»Ol Xl«OI II 00 Xtl J#CC-i)l E2|»GRAB0NE(XtI1,C-l)l E2I X THEN XHAXtvEl BBtlJiCOl E2»«GRaB0Ne(BBIi J.OI E2* AX THEN BBHAXI-EI BnNE(B[I],CC).0,AtI).XTtCCl»H f HH)J HH»RTL(lt»AlI]jtXTtC]iH#HH)| NO PEN*CC"1J HbC 1 1 l»HI 80NE(B[I].C).0.RTL(l».AClj)»XT[C]»N»HHJl »A[n»xTCc]'H»HH>) NO PEN"C-H BB[ I ] I --Hi GTR Dl THEN 1 1 ■BBMAX/XMAx ) 2*£PS)*(2*EPS)*XHAX ThEN D2l«ll 5 AND L NEO 1 THEN GO TO ILLI TO 131 00017200? 00017300? 000-17400? 00017500? 00017600? 00017700? 00017800? 00017900? 00018000? 00018100? 00018200? 00018300? OOOI84OO? 00018500? 00018600? 00018700? 00018800? 00018900? 00019000? 00019100? 00019200? 00019300? 00019400? 00019500? 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? 33 endi i r caccsqlvc PREAL VECTOR AA»BB»641l CREAL VECTOR 1NTcH[6*)I CINT I»N»R»OeTE»IT| CREAL OETR.OETIl LOOK u.l.l.N DO aac x } | «ac i ] i cdec mp se(n»aa»detR»oeti»oete»Intch)i CACCS0LvE(NiR.A.AA,INTCH#B»X.B8»IT)| ENDi * OF CLINSYS. SUBRO REGlN PE BEG S* CA END C C P c p L INVIT UTINf % Th X VE % c kCal IN LL SO » I INT I HEAL KEAL HEAL K E AL ABEL MODE LOOP CDEC MODE LOOP BEG! kki. R I ■ 1 I MOOE CACC MODf LOOo BEGi Rn RO EnDi MAXv LOOP IF LOOP If INVITE E SUBRO CTORS B DCCOMPO Su B ROuT RT64( )| F SORT. ■J» I I»D DETR.DE VECTOR VECTOR TR1.TR2 INVIT»N •■PEN L I t • 1 » 1 OMPOSE( I.PEN L I • - 1 • 1 N B2f.Il ll I J.PEN L SOLVECN i.pen l ii.i'i N OT C 1 1 I - OTf I ] *• AlI"ROO I l«2»i mAxval I lll)l ROOTtI RATI0N(CINT N»PCPOlNT WfpCPOINT X.CINT KK)| UTlNE COMPUTES THE EIGENVECTORS ANO PRINCIPAL Y SOLVING ((W-LAM8DA*I)*#K)*X«0. THE SUBROUTINES SE" AND "CACCSOLyE" IN "cLINSyS" ARE USED* INE S 8rT AS rGACPe REA L ARQ A3 R^A)I ETE.R'L* TI, MAXVAL. XXI, XX2, W2-B. tf2.BB»XXC64] , P»ROOT[6«)I I EXT.FINISI SS N*N| • N DO "2CIJI»«U)» N,M2.0ETR*OETI*DETE*P)l SS 21 • N DO i«i. o» Btni.B2r.nj enoi SS N*NI »R, m ,n2»P»B»X»BB»L)I SS 21 • N DO ROWSUM(X[ 1 1*XC 1 1 > | SQRT(ROOTCI])l T[lll • N DO GEO HOOTCH THEN MAXVAL I»mAXVAL ELSE MAXVal I «ROOT[ H I • N DO ]. MAXVAL THEN BEGIN III. II GO TO NEXT* ENUf mEXT I XxI|«GRABONE(Xt in»0)| XX2l«-GRAB0NEl IF -Xx2 NEQ THEN BEGIN LOOp I • ■!• 1#N DO BEGIN 00022900? 00023000? 00023100? 00023200? 00023300? 00023400? 00023300? 00023600? 00023700? 00023800? 00023900? 00024000? 00024100? 00024200? 00024300? 00024400? 00024500? 00024600? 00024700? 00024 8 00? 00024900? 00025000? 00025100? 00025200? 00025300? 00025400? 00025500? 00025600? 00025700? 00025800? 00025900? 00026000? 00026100? 00026200? 00026300? 00026A00? 00026500? 00026600? 00026700? 00026800? 00026900? 00027000? 00027100? 00027200? 00027300? 00027400? 00027500? 00027600? 00027700? 00027800? 00027900? 00028000? 00028100? 00028200? 00028300? 00028400? 00028500? 3^ MoDEl'PEN LSS 21 TRl» ft 6RABnNE(Xtl)»0)*XXl"GRAB NCGRAB0"E(xm»0)*XX2 + GRAB0NE XX(I]I«TR1| MoOEl«PEN»ll XX[n«"TR2| END) MOOEInPEN LSS 21 Xxl»«GRABONE(XX[IIJ»0>l LOUP Il"l»l*N DO XCUI«XXE Il/XXll END ELSE LOOP Ii.l.l'N 00 Xtlll.xm/XXll LOOP IlMl*lfN 00 BEGjN IF ABS(B2[I).X(I}) GTR P. 10 THEN BEGIN IF KK GTR 10 THEN GO TO pINISl Loop ji«i»i.n oo begin b c z 3 t«x c t i s B2CI1I-XCI1I endi KKI"KK*1I Go TO INVITI end else go to finisi ENOi riNISi ENOI UF INVITERATION, 0002860( 0002870' 00028600 00028900 00029000' 00029100' 00029200' 00029300' 00029*00' 00029500' 00029600* 00029700* 00029800' 00029900' 000300001 00030100! 000302001 00030300' 00030400' 00030500! 00030600' 00030700! 00030800! 00030900' 00031000! 00031100' UNCLASSIFIED Security Classification DOCUMENT CONTROL DATA -R&D (S.cuHtr cla.mltlalon at t ill,, body at ab.rmcl and lmHut*4 -mo tatl an mu,l b* antarad whan tha o„, mU tmpot< ,. cla:lfl.d) ORIGINATING ACTIVITY (Corpora f author) Center for Advanced Computation University of Illinois at Urbana-Champaign Urbana, Illinois 6l801 3. REPORT TITLE 2«. REPORT SECURITY C L A ill FIC A TIO* UNCLASSIFIED 2b. GROUP A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS 4. descriptive NOTES (Typ* ol rapott and toeiut/n dstam) Research Report 8 AUTHOR(S) (Firml nam—, mlddl* Initial, laml nam*) Killion Noh «. REPORT DATE January 20, 1972 IB. TOTAL NO- OF PAGES 3L 7b. NO. OP REFS -L£ •a. CONTRACT OR GRANT NO. DAHCOU 72-C-OOOl b. PROJECT NO. ARPA Order 1899 CAC Document No. 2k •b. OTHER REPORT NOI1I (Any othar nimbara thai may b» mmtlgnmd thla report) UIUCDCS-R-72-1+98 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 WORD* HOLE *T HOLE 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 4. Title and Subtitle A NUMERICAL SOLUTION OF THE MATRIX RICCATI EQUATIONS 7. Author(s) Killion Noh 3. Recipient's Accession No. 5- Report Date January 20, 1972 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. DAHC0U 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. 17a. 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 from the address in (9) above. Distribution unlimited. 19. Security Class (This Report) TINT I ASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 38 22. Price NC FORM NTIS-35 ( 10-70) USCOMM-DC 40329-P7I 4^ JUL26 A973 UNIVERSITY OF ILLINOIS-URBANA 510.84 IL6R no. C002 no. 493-498(1972 Tree height reduction tor parallel proce 3 01 12 088400210 ■wSHnhH rM ■ran rami H BBS BB I 9 H KB SDH MHQj ■■■1 ■ : ■ --33 I *^ v 7 .V. Y/