1 /* 2 3 Copyright (C) 2009-2016 Lukas F. Reichlin 4 5 This file is part of LTI Syncope. 6 7 LTI Syncope is free software: you can redistribute it and/or modify 8 it under the terms of the GNU General Public License as published by 9 the Free Software Foundation, either version 3 of the License, or 10 (at your option) any later version. 11 12 LTI Syncope is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU General Public License for more details. 16 17 You should have received a copy of the GNU General Public License 18 along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. 19 20 Solution of algebraic Riccati equations for descriptor systems. 21 Uses SLICOT SG02AD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: October 2010 26 Version: 0.4 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 #include <complex> 33 34 extern "C" 35 { 36 int F77_FUNC (sg02ad, SG02AD) 37 (char& DICO, char& JOBB, 38 char& FACT, char& UPLO, 39 char& JOBL, char& SCAL, 40 char& SORT, char& ACC, 41 F77_INT& N, F77_INT& M, F77_INT& P, 42 double* A, F77_INT& LDA, 43 double* E, F77_INT& LDE, 44 double* B, F77_INT& LDB, 45 double* Q, F77_INT& LDQ, 46 double* R, F77_INT& LDR, 47 double* L, F77_INT& LDL, 48 double& RCONDU, 49 double* X, F77_INT& LDX, 50 double* ALFAR, double* ALFAI, 51 double* BETA, 52 double* S, F77_INT& LDS, 53 double* T, F77_INT& LDT, 54 double* U, F77_INT& LDU, 55 double& TOL, 56 F77_INT* IWORK, 57 double* DWORK, F77_INT& LDWORK, 58 F77_LOGICAL* BWORK, 59 F77_INT& IWARN, F77_INT& INFO); 60 } 61 62 // PKG_ADD: autoload ("__sl_sg02ad__", "__control_slicot_functions__.oct"); 63 DEFUN_DLD (__sl_sg02ad__, args, nargout, 64 "-*- texinfo -*-\n\ 65 Slicot SG02AD Release 5.0\n\ 66 No argument checking.\n\ 67 For internal use only.") 68 { 69 octave_idx_type nargin = args.length (); 70 octave_value_list retval; 71 72 if (nargin != 8) 73 { 74 print_usage (); 75 } 76 else 77 { 78 // arguments in 79 char dico; 80 char jobb = 'B'; 81 char fact = 'N'; 82 char uplo = 'U'; 83 char jobl; 84 char scal = 'N'; 85 char sort = 'S'; 86 char acc = 'N'; 87 88 Matrix a = args(0).matrix_value (); 89 Matrix e = args(1).matrix_value (); 90 Matrix b = args(2).matrix_value (); 91 Matrix q = args(3).matrix_value (); 92 Matrix r = args(4).matrix_value (); 93 Matrix l = args(5).matrix_value (); 94 F77_INT discrete = args(6).int_value (); 95 F77_INT ijobl = args(7).int_value (); 96 97 if (discrete == 0) 98 dico = 'C'; 99 else 100 dico = 'D'; 101 102 if (ijobl == 0) 103 jobl = 'Z'; 104 else 105 jobl = 'N'; 106 107 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 108 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 109 F77_INT p = 0; // p: number of outputs, not used because FACT = 'N' 110 111 F77_INT lda = max (1, n); 112 F77_INT lde = max (1, n); 113 F77_INT ldb = max (1, n); 114 F77_INT ldq = max (1, n); 115 F77_INT ldr = max (1, m); 116 F77_INT ldl = max (1, n); 117 118 // arguments out 119 double rcondu; 120 121 F77_INT ldx = max (1, n); 122 Matrix x (ldx, n); 123 124 F77_INT nu = 2*n; 125 ColumnVector alfar (nu); 126 ColumnVector alfai (nu); 127 ColumnVector beta (nu); 128 129 // unused output arguments 130 F77_INT lds = max (1, 2*n + m); 131 OCTAVE_LOCAL_BUFFER (double, s, lds * lds); 132 133 F77_INT ldt = max (1, 2*n + m); 134 OCTAVE_LOCAL_BUFFER (double, t, ldt * 2*n); 135 136 F77_INT ldu = max (1, 2*n); 137 OCTAVE_LOCAL_BUFFER (double, u, ldu * 2*n); 138 139 // tolerance 140 double tol = 0; // use default value 141 142 // workspace 143 F77_INT liwork = max (1, m, 2*n); 144 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 145 146 F77_INT ldwork = max (7*(2*n + 1) + 16, 16*n, 2*n + m, 3*m); 147 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 148 149 OCTAVE_LOCAL_BUFFER (F77_LOGICAL, bwork, 2*n); 150 151 // error indicator 152 F77_INT iwarn; 153 F77_INT info; 154 155 156 // SLICOT routine SG02AD 157 F77_XFCN (sg02ad, SG02AD, 158 (dico, jobb, 159 fact, uplo, 160 jobl, scal, 161 sort, acc, 162 n, m, p, 163 a.fortran_vec (), lda, 164 e.fortran_vec (), lde, 165 b.fortran_vec (), ldb, 166 q.fortran_vec (), ldq, 167 r.fortran_vec (), ldr, 168 l.fortran_vec (), ldl, 169 rcondu, 170 x.fortran_vec (), ldx, 171 alfar.fortran_vec (), alfai.fortran_vec (), 172 beta.fortran_vec (), 173 s, lds, 174 t, ldt, 175 u, ldu, 176 tol, 177 iwork, 178 dwork, ldwork, 179 bwork, 180 iwarn, info)); 181 182 if (f77_exception_encountered) 183 error ("are: __sl_sg02ad__: exception in SLICOT subroutine SG02AD"); 184 185 static const char* err_msg[] = { 186 "0: OK", 187 "1: the computed extended matrix pencil is singular, " 188 "possibly due to rounding errors", 189 "2: the QZ algorithm failed", 190 "3: reordering of the generalized eigenvalues failed", 191 "4: after reordering, roundoff changed values of " 192 "some complex eigenvalues so that leading eigenvalues " 193 "in the generalized Schur form no longer satisfy the " 194 "stability condition; this could also be caused due " 195 "to scaling", 196 "5: the computed dimension of the solution does not " 197 "equal N", 198 "6: the spectrum is too close to the boundary of " 199 "the stability domain", 200 "7: a singular matrix was encountered during the " 201 "computation of the solution matrix X"}; 202 203 static const char* warn_msg[] = { 204 "0: OK", 205 "1: solution may be inaccurate due to poor scaling " 206 "or eigenvalues too close to the boundary of the stability domain " 207 "(the imaginary axis, if DICO = 'C', or the unit circle, if DICO = 'D')"}; 208 209 error_msg ("are", info, 7, err_msg); 210 warning_msg ("are", iwarn, 1, warn_msg); 211 212 213 // assemble complex vector - adapted from DEFUN complex in data.cc 214 alfar.resize (n); 215 alfai.resize (n); 216 beta.resize (n); 217 218 ColumnVector poler (n); 219 ColumnVector polei (n); 220 221 poler = quotient (alfar, beta); 222 polei = quotient (alfai, beta); 223 224 ComplexColumnVector pole (n, Complex ()); 225 226 for (F77_INT i = 0; i < n; i++) 227 pole.xelem (i) = Complex (poler(i), polei(i)); 228 229 // return value 230 retval(0) = x; 231 retval(1) = pole; 232 } 233 234 return retval; 235 } 236