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 realization of state-space models. 21 Uses SLICOT TB01PD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: September 2010 26 Version: 0.5 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 33 extern "C" 34 { 35 int F77_FUNC (tb01pd, TB01PD) 36 (char& JOB, char& EQUIL, 37 F77_INT& N, F77_INT& M, F77_INT& P, 38 double* A, F77_INT& LDA, 39 double* B, F77_INT& LDB, 40 double* C, F77_INT& LDC, 41 F77_INT& NR, 42 double& TOL, 43 F77_INT* IWORK, 44 double* DWORK, F77_INT& LDWORK, 45 F77_INT& INFO); 46 } 47 48 // PKG_ADD: autoload ("__sl_tb01pd__", "__control_slicot_functions__.oct"); 49 DEFUN_DLD (__sl_tb01pd__, args, nargout, 50 "-*- texinfo -*-\n\ 51 Slicot TB01PD 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 != 5) 59 { 60 print_usage (); 61 } 62 else 63 { 64 // arguments in 65 char job = 'M'; 66 char equil; 67 68 Matrix a = args(0).matrix_value (); 69 Matrix b = args(1).matrix_value (); 70 Matrix c = args(2).matrix_value (); 71 double tol = args(3).double_value (); 72 const F77_INT scaled = args(4).int_value (); 73 74 if (scaled == 0) 75 equil = 'S'; 76 else 77 equil = 'N'; 78 79 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 80 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 81 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 82 83 F77_INT lda = max (1, n); 84 F77_INT ldb = max (1, n); 85 F77_INT ldc; 86 87 if (n == 0) 88 ldc = 1; 89 else 90 ldc = max (1, m, p); 91 92 b.resize (ldb, max (m, p)); 93 c.resize (ldc, n); 94 95 // arguments out 96 F77_INT nr = 0; 97 98 // workspace 99 F77_INT liwork = n + max (m, p); 100 F77_INT ldwork = max (1, n + max (n, 3*m, 3*p)); 101 102 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 103 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 104 105 // error indicators 106 F77_INT info = 0; 107 108 109 // SLICOT routine TB01PD 110 F77_XFCN (tb01pd, TB01PD, 111 (job, equil, 112 n, m, p, 113 a.fortran_vec (), lda, 114 b.fortran_vec (), ldb, 115 c.fortran_vec (), ldc, 116 nr, 117 tol, 118 iwork, 119 dwork, ldwork, 120 info)); 121 122 if (f77_exception_encountered) 123 error ("ss: minreal: __sl_tb01pd__: exception in SLICOT subroutine TB01PD"); 124 125 if (info != 0) 126 error ("ss: minreal: __sl_tb01pd__: TB01PD returned info = %d", static_cast<int> (info)); 127 128 // resize 129 a.resize (nr, nr); 130 b.resize (nr, m); 131 c.resize (p, nr); 132 133 // return values 134 retval(0) = a; 135 retval(1) = b; 136 retval(2) = c; 137 retval(3) = octave_value (nr); 138 } 139 140 return retval; 141 } 142