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 Transfer matrix of a given state-space representation (A,B,C,D), 21 using the pole-zeros method. 22 Uses SLICOT TB04BD by courtesy of NICONET e.V. 23 <http://www.slicot.org> 24 25 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 26 Created: October 2010 27 Version: 0.3 28 29 */ 30 31 #include <octave/oct.h> 32 #include "common.h" 33 34 extern "C" 35 { 36 int F77_FUNC (tb04bd, TB04BD) 37 (char& JOBD, char& ORDER, char& EQUIL, 38 F77_INT& N, F77_INT& M, F77_INT& P, F77_INT& MD, 39 double* A, F77_INT& LDA, 40 double* B, F77_INT& LDB, 41 double* C, F77_INT& LDC, 42 double* D, F77_INT& LDD, 43 F77_INT* IGN, F77_INT& LDIGN, 44 F77_INT* IGD, F77_INT& LDIGD, 45 double* GN, double* GD, 46 double& TOL, 47 F77_INT* IWORK, 48 double* DWORK, F77_INT& LDWORK, 49 F77_INT& INFO); 50 } 51 52 // PKG_ADD: autoload ("__sl_tb04bd__", "__control_slicot_functions__.oct"); 53 DEFUN_DLD (__sl_tb04bd__, args, nargout, 54 "-*- texinfo -*-\n\ 55 Slicot TB04BD Release 5.0\n\ 56 No argument checking.\n\ 57 For internal use only.") 58 { 59 octave_idx_type nargin = args.length (); 60 octave_value_list retval; 61 62 if (nargin != 5) 63 { 64 print_usage (); 65 } 66 else 67 { 68 // arguments in 69 char jobd = 'D'; 70 char order = 'D'; 71 char equil; 72 73 Matrix a = args(0).matrix_value (); 74 Matrix b = args(1).matrix_value (); 75 Matrix c = args(2).matrix_value (); 76 Matrix d = args(3).matrix_value (); 77 const F77_INT scaled = args(4).int_value (); 78 79 if (scaled == 0) 80 equil = 'S'; 81 else 82 equil = 'N'; 83 84 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 85 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 86 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 87 F77_INT md = n + 1; 88 89 F77_INT lda = max (1, n); 90 F77_INT ldb = max (1, n); 91 F77_INT ldc = max (1, p); 92 F77_INT ldd = max (1, p); 93 94 // arguments out 95 F77_INT ldign = max (1, p); 96 F77_INT ldigd = max (1, p); 97 F77_INT lg = p * m * md; 98 99 OCTAVE_LOCAL_BUFFER (F77_INT, ign, ldign*m); 100 OCTAVE_LOCAL_BUFFER (F77_INT, igd, ldigd*m); 101 102 RowVector gn (lg); 103 RowVector gd (lg); 104 105 // tolerance 106 double tol = 0; // use default value 107 108 // workspace 109 F77_INT ldwork = max (1, n*(n + p) + max (n + max (n, p), n*(2*n + 5))); 110 111 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, n); 112 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 113 114 // error indicator 115 F77_INT info; 116 117 118 // SLICOT routine TB04BD 119 F77_XFCN (tb04bd, TB04BD, 120 (jobd, order, equil, 121 n, m, p, md, 122 a.fortran_vec (), lda, 123 b.fortran_vec (), ldb, 124 c.fortran_vec (), ldc, 125 d.fortran_vec (), ldd, 126 ign, ldign, 127 igd, ldigd, 128 gn.fortran_vec (), gd.fortran_vec (), 129 tol, 130 iwork, 131 dwork, ldwork, 132 info)); 133 134 if (f77_exception_encountered) 135 error ("ss2tf: __sl_tb04bd__: exception in SLICOT subroutine TB04BD"); 136 137 if (info != 0) 138 error ("ss2tf: __sl_tb04bd__: TB04BD returned info = %d", static_cast<int> (info)); 139 140 141 // return values 142 Cell num(p, m); 143 Cell den(p, m); 144 F77_INT ik, istr; 145 146 for (ik = 0; ik < p*m; ik++) 147 { 148 istr = ik*md; 149 num.xelem (ik) = gn.extract_n (istr, ign[ik]+1); 150 den.xelem (ik) = gd.extract_n (istr, igd[ik]+1); 151 } 152 153 retval(0) = num; 154 retval(1) = den; 155 } 156 157 return retval; 158 } 159