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 Minimal state-space representation (A,B,C,D) for a 21 proper transfer matrix T(s) given as either row or column 22 polynomial vectors over denominator polynomials, possibly with 23 uncancelled common terms. 24 Uses SLICOT TD04AD by courtesy of NICONET e.V. 25 <http://www.slicot.org> 26 27 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 28 Created: August 2011 29 Version: 0.3 30 31 */ 32 33 #include <octave/oct.h> 34 #include "common.h" 35 36 extern "C" 37 { 38 int F77_FUNC (td04ad, TD04AD) 39 (char& ROWCOL, 40 F77_INT& M, F77_INT& P, 41 F77_INT* INDEX, 42 double* DCOEFF, F77_INT& LDDCOE, 43 double* UCOEFF, F77_INT& LDUCO1, F77_INT& LDUCO2, 44 F77_INT& NR, 45 double* A, F77_INT& LDA, 46 double* B, F77_INT& LDB, 47 double* C, F77_INT& LDC, 48 double* D, F77_INT& LDD, 49 double& TOL, 50 F77_INT* IWORK, 51 double* DWORK, F77_INT& LDWORK, 52 F77_INT& INFO); 53 } 54 55 // PKG_ADD: autoload ("__sl_td04ad__", "__control_slicot_functions__.oct"); 56 DEFUN_DLD (__sl_td04ad__, args, nargout, 57 "-*- texinfo -*-\n\ 58 Slicot TD04AD Release 5.0\n\ 59 No argument checking.\n\ 60 For internal use only.") 61 { 62 octave_idx_type nargin = args.length (); 63 octave_value_list retval; 64 65 if (nargin != 4) 66 { 67 print_usage (); 68 } 69 else 70 { 71 // arguments in 72 char rowcol = 'R'; 73 74 NDArray ucoeff = args(0).array_value (); 75 Matrix dcoeff = args(1).matrix_value (); 76 Matrix indexd = args(2).matrix_value (); 77 double tol = args(3).double_value (); 78 79 F77_INT p = TO_F77_INT (ucoeff.rows ()); // p: number of outputs 80 F77_INT m = TO_F77_INT (ucoeff.columns ()); // m: number of inputs 81 82 F77_INT lddcoe = max (1, p); // TODO: handle case ucoeff.rows = 0 83 F77_INT lduco1 = max (1, p); 84 F77_INT lduco2 = max (1, m); 85 86 F77_INT n = 0; 87 OCTAVE_LOCAL_BUFFER (F77_INT, index, p); 88 for (F77_INT i = 0; i < p; i++) 89 { 90 index[i] = TO_F77_INT (indexd.xelem (i)); 91 n += index[i]; 92 } 93 94 // arguments out 95 F77_INT nr = max (1, n); // initialize to prevent crash if info != 0 96 F77_INT lda = max (1, n); 97 F77_INT ldb = max (1, n); 98 F77_INT ldc = max (1, m, p); 99 F77_INT ldd = max (1, p); 100 101 Matrix a (lda, n); 102 Matrix b (ldb, max (m, p)); 103 Matrix c (ldc, n); 104 Matrix d (ldd, m); 105 106 // workspace 107 F77_INT ldwork = max (1, n + max (n, 3*m, 3*p)); 108 109 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, n + max (m, p)); 110 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 111 112 // error indicator 113 F77_INT info; 114 115 116 // SLICOT routine TD04AD 117 F77_XFCN (td04ad, TD04AD, 118 (rowcol, 119 m, p, 120 index, 121 dcoeff.fortran_vec (), lddcoe, 122 ucoeff.fortran_vec (), lduco1, lduco2, 123 nr, 124 a.fortran_vec (), lda, 125 b.fortran_vec (), ldb, 126 c.fortran_vec (), ldc, 127 d.fortran_vec (), ldd, 128 tol, 129 iwork, 130 dwork, ldwork, 131 info)); 132 133 if (f77_exception_encountered) 134 error ("tf2ss: __sl_td04ad__: exception in SLICOT subroutine TD04AD"); 135 136 if (info != 0) 137 error ("tf2ss: __sl_td04ad__: TD04AD returned info = %d", static_cast<int> (info)); 138 139 // resize 140 a.resize (nr, nr); 141 b.resize (nr, m); 142 c.resize (p, nr); 143 d.resize (p, m); 144 145 // return values 146 retval(0) = a; 147 retval(1) = b; 148 retval(2) = c; 149 retval(3) = d; 150 } 151 152 return retval; 153 } 154