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 Gain of descriptor state-space models. Based on SLICOT TB04BX.f. 21 22 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 23 Created: March 2011 24 Version: 0.2 25 26 */ 27 28 #include <octave/oct.h> 29 #include "common.h" 30 #include <complex> 31 #include <xpow.h> 32 33 extern "C" 34 { 35 int F77_FUNC (tg04bx, TG04BX) 36 (F77_INT& IP, F77_INT& IZ, 37 double* A, F77_INT& LDA, 38 double* E, 39 double* B, 40 double* C, 41 double* D, 42 double* PR, double* PI, 43 double* ZR, double* ZI, 44 double& GAIN, 45 F77_INT* IWORK); 46 } 47 48 // PKG_ADD: autoload ("__sl_tg04bx__", "__control_slicot_functions__.oct"); 49 DEFUN_DLD (__sl_tg04bx__, args, nargout, 50 "-*- texinfo -*-\n\ 51 Slicot TG04BX Release 5.0\n\ 52 No argument checking.\n\ 53 For internal use only.") 54 { 55 octave_idx_type nargin = args.length (); 56 octave_value_list retval; 57 58 if (nargin != 9) 59 { 60 print_usage (); 61 } 62 else 63 { 64 // arguments in 65 Matrix a = args(0).matrix_value (); 66 Matrix e = args(1).matrix_value (); 67 Matrix b = args(2).matrix_value (); 68 Matrix c = args(3).matrix_value (); 69 Matrix d = args(4).matrix_value (); 70 71 ColumnVector pr = args(5).column_vector_value (); 72 ColumnVector pi = args(6).column_vector_value (); 73 74 ColumnVector zr = args(7).column_vector_value (); 75 ColumnVector zi = args(8).column_vector_value (); 76 77 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 78 F77_INT ip = TO_F77_INT (pr.numel ()); // ip: number of finite poles 79 F77_INT iz = TO_F77_INT (zr.numel ()); // iz: number of zeros 80 81 // For ss, IP = n is always true. 82 // However, dss models with poles at infinity 83 // (filtered by pole.m) may have IP <= n 84 85 // Take pr.length == pi.length == ip for granted, 86 // and the same for iz, zr and zi. 87 88 F77_INT lda = max (1, n); 89 90 // arguments out 91 double gain; 92 93 // workspace 94 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, lda); 95 96 97 F77_XFCN (tg04bx, TG04BX, 98 (ip, iz, 99 a.fortran_vec (), lda, 100 e.fortran_vec (), 101 b.fortran_vec (), 102 c.fortran_vec (), 103 d.fortran_vec (), 104 pr.fortran_vec (), pi.fortran_vec (), 105 zr.fortran_vec (), zi.fortran_vec (), 106 gain, 107 iwork)); 108 109 if (f77_exception_encountered) 110 error ("dss: zero: __sl_tg04bx__: exception in TG04BX"); 111 112 // return values 113 retval(0) = octave_value (gain); 114 } 115 116 return retval; 117 } 118