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 Staircase observability form for descriptor models. 21 Uses SLICOT TG01ID 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.2 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 33 extern "C" 34 { 35 int F77_FUNC (tg01id, TG01ID) 36 (char& JOBOBS, 37 char& COMPQ, char& COMPZ, 38 F77_INT& N, F77_INT& M, F77_INT& P, 39 double* A, F77_INT& LDA, 40 double* E, F77_INT& LDE, 41 double* B, F77_INT& LDB, 42 double* C, F77_INT& LDC, 43 double* Q, F77_INT& LDQ, 44 double* Z, F77_INT& LDZ, 45 F77_INT& NOBSV, F77_INT& NIUOBS, 46 F77_INT& NLBLCK, 47 F77_INT* CTAU, 48 double& TOL, 49 F77_INT* IWORK, double* DWORK, 50 F77_INT& INFO); 51 } 52 53 // PKG_ADD: autoload ("__sl_tg01id__", "__control_slicot_functions__.oct"); 54 DEFUN_DLD (__sl_tg01id__, args, nargout, "Slicot TG01ID Release 5.0") 55 { 56 octave_idx_type nargin = args.length (); 57 octave_value_list retval; 58 59 if (nargin != 5) 60 { 61 print_usage (); 62 } 63 else 64 { 65 // arguments in 66 char jobobs = 'O'; 67 char compq = 'I'; 68 char compz = 'I'; 69 70 Matrix a = args(0).matrix_value (); 71 Matrix e = args(1).matrix_value (); 72 Matrix b = args(2).matrix_value (); 73 Matrix c = args(3).matrix_value (); 74 double tol = args(4).double_value (); 75 76 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 77 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 78 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 79 80 F77_INT lda = max (1, n); 81 F77_INT lde = max (1, n); 82 F77_INT ldb = max (1, n); 83 F77_INT ldc = max (1, m, p); 84 F77_INT ldq = max (1, n); 85 F77_INT ldz = max (1, n); 86 87 b.resize (ldb, max (m, p)); 88 89 // arguments out 90 Matrix q (ldq, n); 91 Matrix z (ldz, n); 92 93 F77_INT nobsv; 94 F77_INT niuobs; 95 F77_INT nlblck; 96 97 OCTAVE_LOCAL_BUFFER (F77_INT, ctau, n); 98 99 // workspace 100 F77_INT ldwork = max (n, 2*p); 101 102 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, p); 103 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 104 105 // error indicators 106 F77_INT info; 107 108 109 // SLICOT routine TG01ID 110 F77_XFCN (tg01id, TG01ID, 111 (jobobs, 112 compq, compz, 113 n, m, p, 114 a.fortran_vec (), lda, 115 e.fortran_vec (), lde, 116 b.fortran_vec (), ldb, 117 c.fortran_vec (), ldc, 118 q.fortran_vec (), ldq, 119 z.fortran_vec (), ldz, 120 nobsv, niuobs, 121 nlblck, 122 ctau, 123 tol, 124 iwork, dwork, 125 info)); 126 127 if (f77_exception_encountered) 128 error ("__sl_tg01id__: exception in SLICOT subroutine TG01ID"); 129 130 if (info != 0) 131 error ("__sl_tg01id__: TG01ID returned info = %d", static_cast<int> (info)); 132 133 // resize 134 a.resize (n, n); 135 e.resize (n, n); 136 b.resize (n, m); 137 c.resize (p, n); 138 q.resize (n, n); 139 z.resize (n, n); 140 141 // return values 142 retval(0) = a; 143 retval(1) = e; 144 retval(2) = b; 145 retval(3) = c; 146 retval(4) = q; 147 retval(5) = z; 148 retval(6) = octave_value (nobsv); 149 } 150 151 return retval; 152 } 153