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