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 Model reduction based on Hankel-norm approximation method. 21 Uses SLICOT AB09JD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: July 2011 26 Version: 0.2 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 33 extern "C" 34 { 35 int F77_FUNC (ab09jd, AB09JD) 36 (char& JOBV, char& JOBW, char& JOBINV, 37 char& DICO, char& EQUIL, char& ORDSEL, 38 F77_INT& N, F77_INT& NV, F77_INT& NW, F77_INT& M, F77_INT& P, 39 F77_INT& NR, 40 double& ALPHA, 41 double* A, F77_INT& LDA, 42 double* B, F77_INT& LDB, 43 double* C, F77_INT& LDC, 44 double* D, F77_INT& LDD, 45 double* AV, F77_INT& LDAV, 46 double* BV, F77_INT& LDBV, 47 double* CV, F77_INT& LDCV, 48 double* DV, F77_INT& LDDV, 49 double* AW, F77_INT& LDAW, 50 double* BW, F77_INT& LDBW, 51 double* CW, F77_INT& LDCW, 52 double* DW, F77_INT& LDDW, 53 F77_INT& NS, 54 double* HSV, 55 double& TOL1, double& TOL2, 56 F77_INT* IWORK, 57 double* DWORK, F77_INT& LDWORK, 58 F77_INT& IWARN, F77_INT& INFO); 59 } 60 61 // PKG_ADD: autoload ("__sl_ab09jd__", "__control_slicot_functions__.oct"); 62 DEFUN_DLD (__sl_ab09jd__, args, nargout, 63 "-*- texinfo -*-\n\ 64 Slicot AB09JD Release 5.0\n\ 65 No argument checking.\n\ 66 For internal use only.") 67 { 68 octave_idx_type nargin = args.length (); 69 octave_value_list retval; 70 71 if (nargin != 22) 72 { 73 print_usage (); 74 } 75 else 76 { 77 // arguments in 78 char jobv; 79 char jobw; 80 char jobinv; 81 char dico; 82 char equil; 83 char ordsel; 84 85 Matrix a = args(0).matrix_value (); 86 Matrix b = args(1).matrix_value (); 87 Matrix c = args(2).matrix_value (); 88 Matrix d = args(3).matrix_value (); 89 90 const F77_INT idico = args(4).int_value (); 91 const F77_INT iequil = args(5).int_value (); 92 F77_INT nr = args(6).int_value (); 93 const F77_INT iordsel = args(7).int_value (); 94 double alpha = args(8).double_value (); 95 96 const F77_INT ijobv = args(9).int_value (); 97 Matrix av = args(10).matrix_value (); 98 Matrix bv = args(11).matrix_value (); 99 Matrix cv = args(12).matrix_value (); 100 Matrix dv = args(13).matrix_value (); 101 102 const F77_INT ijobw = args(14).int_value (); 103 Matrix aw = args(15).matrix_value (); 104 Matrix bw = args(16).matrix_value (); 105 Matrix cw = args(17).matrix_value (); 106 Matrix dw = args(18).matrix_value (); 107 108 const F77_INT ijobinv = args(19).int_value (); 109 double tol1 = args(20).double_value (); 110 double tol2 = args(21).double_value (); 111 112 switch (ijobv) 113 { 114 case 0: 115 jobv = 'N'; 116 break; 117 case 1: 118 jobv = 'V'; 119 break; 120 case 2: 121 jobv = 'I'; 122 break; 123 case 3: 124 jobv = 'C'; 125 break; 126 case 4: 127 jobv = 'R'; 128 break; 129 default: 130 error ("__sl_ab09jd__: argument jobv invalid"); 131 } 132 133 switch (ijobw) 134 { 135 case 0: 136 jobw = 'N'; 137 break; 138 case 1: 139 jobw = 'W'; 140 break; 141 case 2: 142 jobw = 'I'; 143 break; 144 case 3: 145 jobw = 'C'; 146 break; 147 case 4: 148 jobw = 'R'; 149 break; 150 default: 151 error ("__sl_ab09jd__: argument jobw invalid"); 152 } 153 154 switch (ijobinv) 155 { 156 case 0: 157 jobinv = 'N'; 158 break; 159 case 1: 160 jobinv = 'I'; 161 break; 162 case 2: 163 jobinv = 'A'; 164 break; 165 default: 166 error ("__sl_ab09jd__: argument jobinv invalid"); 167 } 168 169 if (idico == 0) 170 dico = 'C'; 171 else 172 dico = 'D'; 173 174 if (iequil == 0) 175 equil = 'S'; 176 else 177 equil = 'N'; 178 179 if (iordsel == 0) 180 ordsel = 'F'; 181 else 182 ordsel = 'A'; 183 184 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 185 F77_INT nv = TO_F77_INT (av.rows ()); 186 F77_INT nw = TO_F77_INT (aw.rows ()); 187 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 188 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 189 190 F77_INT lda = max (1, n); 191 F77_INT ldb = max (1, n); 192 F77_INT ldc = max (1, p); 193 F77_INT ldd = max (1, p); 194 195 F77_INT ldav = max (1, nv); 196 F77_INT ldbv = max (1, nv); 197 F77_INT ldcv = max (1, p); 198 F77_INT lddv = max (1, p); 199 200 F77_INT ldaw = max (1, nw); 201 F77_INT ldbw = max (1, nw); 202 F77_INT ldcw = max (1, m); 203 F77_INT lddw = max (1, m); 204 205 // arguments out 206 F77_INT ns; 207 ColumnVector hsv (n); 208 209 // workspace 210 F77_INT liwork; 211 F77_INT tmpc; 212 F77_INT tmpd; 213 214 if (jobv == 'N') 215 tmpc = 0; 216 else 217 tmpc = max (2*p, nv+p+n+6, 2*nv+p+2); 218 219 if (jobw == 'N') 220 tmpd = 0; 221 else 222 tmpd = max (2*m, nw+m+n+6, 2*nw+m+2); 223 224 if (dico == 'C') 225 liwork = max (1, m, tmpc, tmpd); 226 else 227 liwork = max (1, n, m, tmpc, tmpd); 228 229 F77_INT ldwork; 230 F77_INT nvp = nv + p; 231 F77_INT nwm = nw + m; 232 F77_INT ldw1; 233 F77_INT ldw2; 234 F77_INT ldw3 = n*(2*n + max (n, m, p) + 5) + n*(n+1)/2; 235 F77_INT ldw4 = n*(m+p+2) + 2*m*p + min (n, m) + max (3*m+1, min (n, m) + p); 236 237 if (jobv == 'N') 238 { 239 ldw1 = 0; 240 } 241 else 242 { 243 ldw1 = 2*nvp*(nvp+p) + p*p + max (2*nvp*nvp + max (11*nvp+16, p*nvp), 244 nvp*n + max (nvp*n+n*n, p*n, p*m)); 245 } 246 247 if (jobw == 'N') 248 { 249 ldw2 = 0; 250 } 251 else 252 { 253 ldw2 = 2*nwm*(nwm+m) + m*m + max (2*nwm*nwm + max (11*nwm+16, m*nwm), 254 nwm*n + max (nwm*n+n*n, m*n, p*m)); 255 } 256 257 ldwork = max (ldw1, ldw2, ldw3, ldw4); 258 259 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 260 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 261 262 // error indicators 263 F77_INT iwarn = 0; 264 F77_INT info = 0; 265 266 267 // SLICOT routine AB09JD 268 F77_XFCN (ab09jd, AB09JD, 269 (jobv, jobw, jobinv, 270 dico, equil, ordsel, 271 n, nv, nw, m, p, 272 nr, 273 alpha, 274 a.fortran_vec (), lda, 275 b.fortran_vec (), ldb, 276 c.fortran_vec (), ldc, 277 d.fortran_vec (), ldd, 278 av.fortran_vec (), ldav, 279 bv.fortran_vec (), ldbv, 280 cv.fortran_vec (), ldcv, 281 dv.fortran_vec (), lddv, 282 aw.fortran_vec (), ldaw, 283 bw.fortran_vec (), ldbw, 284 cw.fortran_vec (), ldcw, 285 dw.fortran_vec (), lddw, 286 ns, 287 hsv.fortran_vec (), 288 tol1, tol2, 289 iwork, 290 dwork, ldwork, 291 iwarn, info)); 292 293 if (f77_exception_encountered) 294 error ("hnamodred: exception in SLICOT subroutine AB09JD"); 295 296 297 static const char* err_msg[] = { 298 "0: OK", 299 "1: the computation of the ordered real Schur form of A " 300 "failed", 301 302 "2: the separation of the ALPHA-stable/unstable " 303 "diagonal blocks failed because of very close eigenvalues", 304 305 "3: the reduction of AV to a real Schur form failed", 306 307 "4: the reduction of AW to a real Schur form failed", 308 309 "5: the reduction to generalized Schur form of the " 310 "descriptor pair corresponding to the inverse of V " 311 "failed", 312 313 "6: the reduction to generalized Schur form of the " 314 "descriptor pair corresponding to the inverse of W " 315 "failed", 316 317 "7: the computation of Hankel singular values failed", 318 319 "8: the computation of stable projection in the " 320 "Hankel-norm approximation algorithm failed", 321 322 "9: the order of computed stable projection in the " 323 "Hankel-norm approximation algorithm differs " 324 "from the order of Hankel-norm approximation", 325 326 "10: the reduction of AV-BV*inv(DV)*CV to a " 327 "real Schur form failed", 328 329 "11: the reduction of AW-BW*inv(DW)*CW to a " 330 "real Schur form failed", 331 332 "12: the solution of the Sylvester equation failed " 333 "because the poles of V (if JOBV = 'V') or of " 334 "conj(V) (if JOBV = 'C') are not distinct from " 335 "the poles of G1 (see METHOD)", 336 337 "13: the solution of the Sylvester equation failed " 338 "because the poles of W (if JOBW = 'W') or of " 339 "conj(W) (if JOBW = 'C') are not distinct from " 340 "the poles of G1 (see METHOD)", 341 342 "14: the solution of the Sylvester equation failed " 343 "because the zeros of V (if JOBV = 'I') or of " 344 "conj(V) (if JOBV = 'R') are not distinct from " 345 "the poles of G1sr (see METHOD)", 346 347 "15: the solution of the Sylvester equation failed " 348 "because the zeros of W (if JOBW = 'I') or of " 349 "conj(W) (if JOBW = 'R') are not distinct from " 350 "the poles of G1sr (see METHOD)", 351 352 "16: the solution of the generalized Sylvester system " 353 "failed because the zeros of V (if JOBV = 'I') or " 354 "of conj(V) (if JOBV = 'R') are not distinct from " 355 "the poles of G1sr (see METHOD)", 356 357 "17: the solution of the generalized Sylvester system " 358 "failed because the zeros of W (if JOBW = 'I') or " 359 "of conj(W) (if JOBW = 'R') are not distinct from " 360 "the poles of G1sr (see METHOD)", 361 362 "18: op(V) is not antistable", 363 364 "19: op(W) is not antistable", 365 366 "20: V is not invertible", 367 368 "21: W is not invertible"}; 369 370 static const char* warn_msg[] = { 371 "0: OK", 372 "1: with ORDSEL = 'F', the selected order NR is greater " 373 "than NSMIN, the sum of the order of the " 374 "ALPHA-unstable part and the order of a minimal " 375 "realization of the ALPHA-stable part of the given " 376 "system. In this case, the resulting NR is set equal " 377 "to NSMIN.", 378 "2: with ORDSEL = 'F', the selected order NR is less " 379 "than the order of the ALPHA-unstable part of the " 380 "given system. In this case NR is set equal to the " 381 "order of the ALPHA-unstable part."}; 382 383 error_msg ("hnamodred", info, 21, err_msg); 384 warning_msg ("hnamodred", iwarn, 2, warn_msg); 385 386 // resize 387 a.resize (nr, nr); 388 b.resize (nr, m); 389 c.resize (p, nr); 390 hsv.resize (ns); 391 392 // return values 393 retval(0) = a; 394 retval(1) = b; 395 retval(2) = c; 396 retval(3) = d; 397 retval(4) = octave_value (nr); 398 retval(5) = hsv; 399 retval(6) = octave_value (ns); 400 } 401 402 return retval; 403 } 404