1      subroutine rilac(n,nn,a,na,c,d,rcond,x,w,nnw,z,eps,iwrk,wrk1,wrk2,
2     &                 ierr)
3C!purpose
4C
5C        to solve the continuous time algebraic equation
6C
7C                trans(a)*x + x*a + c - x*d*x  =  0
8C
9C        where  trans(a)  denotes the transpose of  a .
10C
11C!method
12C
13C        the method used is laub's variant of the hamiltonian -
14C        eigenvector approach (schur method).
15C
16C!reference
17C
18C        a.j. laub
19C        a schur method for solving algebraic riccati equations
20C        ieee trans. automat. control, vol. ac-25, 1980.
21C
22C! auxiliary routines
23C
24C       orthes,ortran,balanc,balbak (eispack)
25C       dgeco,dgesl (linpack)
26C       hqror2,inva,exchgn,qrstep
27C
28C! calling sequence
29C        subroutine rilac(n,nn,a,na,c,d,rcond,x,w,nnw,z,
30C    +                iwrk,wrk1,wrk2,ierr)
31C
32C        integer n,nn,na,nnw,iwrk(nn),ierr
33C        double precision a(na,n),c(na,n),d(na,n)
34C        double precision rcond,x(na,n),w(nnw,nn),z(nnw,nn)
35C        double precision wrk1(nn),wrk2(nn)
36C
37C arguments in
38C
39C       n       integer
40C               -the order of a,c,d and x
41C
42C       na      integer
43C               -the declared first dimension of a,c,d and x
44C
45C       nn      integer
46C               -the order of w and z
47C                    nn = n + n
48C
49C       nnw     integer
50C               -the declared first dimension of w and z
51C
52C
53C       a       double precision(n,n)
54C
55C       c       double precision(n,n)
56C
57C       d       double precision(n,n)
58C
59C arguments out
60C
61C       x       double precision(n,n)
62C               - x  contains the solution matrix
63C
64C       w       double precision(nn,nn)
65C               - w  contains the ordered real upper-triangular
66C               form of the hamiltonian matrix
67C
68C       z       double precision(nn,nn)
69C               - z  contains the transformation matrix which
70C               reduces the hamiltonian matrix to the ordered
71C               real upper-triangular form
72C
73C       rcond   double precision
74C               - rcond  contains an estimate of the reciprocal
75C               condition of the  n-th order system of algebraic
76C               equations from which the solution matrix is obtained
77C
78C       ierr    integer
79C               -error indicator set on exit
80C
81C               ierr  =  0       successful return
82C
83C               ierr  =  1       the real upper triangular form of
84C                                the hamiltonian matrix cannot be
85C                                appropriately ordered
86C
87C               ierr  =  2       the hamiltonian matrix has less than n
88C                                eigenvalues with negative real parts
89C
90C               ierr  =  3       the  n-th order system of linear
91C                                algebraic equations, from which the
92C                                solution matrix would be obtained, is
93C                                singular to working precision
94C
95C               ierr  =  4       the hamiltonian matrix cannot be
96C                                reduced to upper-triangular form
97C
98C working space
99C
100C       iwrk    integer(nn)
101C
102C       wrk1    double precision(nn)
103C
104C       wrk2    double precision(nn)
105C
106C!originator
107C
108C                control systems research group, dept. eecs, kingston
109C                polytechnic, penrhyn rd.,kingston-upon-thames, england.
110C
111C! comments
112C                if there is a shortage of storage space, then the
113C                matrices  c  and  x  can share the same locations,
114C                but this will, of course, result in the loss of  c.
115C
116C*******************************************************************
117C
118      integer n,nn,na,nnw,iwrk(nn),ierr
119      double precision a(na,n),c(na,n),d(na,n)
120      double precision rcond,x(na,n),w(nnw,nn),z(nnw,nn)
121      double precision wrk1(nn),wrk2(nn)
122C
123C        local declarations:
124C
125      integer i,j,low,igh,ni,nj
126      double precision eps,t(1)
127      integer folhp
128      external folhp
129      logical fail
130C
131C
132C         eps is a machine dependent parameter specifying
133C        the relative precision of realing point arithmetic.
134C
135C        initialise the hamiltonian matrix associated with the problem
136C
137      do 10 j = 1,n
138        nj = n + j
139        do 10 i = 1,n
140          ni = n + i
141          w(i,j) = a(i,j)
142          w(ni,j) = -c(i,j)
143          w(i,nj) = -d(i,j)
144          w(ni,nj) = -a(j,i)
145 10   continue
146C
147      call balanc(nnw,nn,w,low,igh,wrk1)
148C
149      call orthes(nn,nn,1,nn,w,wrk2)
150      call ortran(nn,nn,1,nn,w,wrk2,z)
151      call hqror2(nn,nn,1,nn,w,t,t,z,ierr,11)
152      if (ierr .ne. 0) goto 70
153      call inva(nn,nn,w,z,folhp,eps,ndim,fail,iwrk)
154C
155      if (ierr .ne. 0) goto 40
156      if (ndim .ne. n) goto 50
157C
158      call balbak(nnw,nn,low,igh,wrk1,nn,z)
159C
160C
161      call dgeco(z,nnw,n,iwrk,rcond,wrk1)
162      if (rcond .lt. eps) goto 60
163C
164      do 30 j = 1,n
165        nj = n + j
166        do 20 i = 1,n
167          x(i,j) = z(nj,i)
168 20     continue
169        call dgesl(z,nnw,n,iwrk,x(1,j),1)
170 30   continue
171      goto 100
172C
173 40   ierr = 1
174      goto 100
175C
176 50   ierr = 2
177      goto 100
178C
179 60   ierr = 3
180      goto 100
181C
182 70   ierr = 4
183      goto 100
184C
185 100  return
186      end
187
188