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