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 Convert descriptor state-space system into regular state-space form. 21 Uses SLICOT SB10JD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: September 2011 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 (sb10jd, SB10JD) 36 (F77_INT& N, F77_INT& M, F77_INT& NP, 37 double* A, F77_INT& LDA, 38 double* B, F77_INT& LDB, 39 double* C, F77_INT& LDC, 40 double* D, F77_INT& LDD, 41 double* E, F77_INT& LDE, 42 F77_INT& NSYS, 43 double* DWORK, F77_INT& LDWORK, 44 F77_INT& INFO); 45 } 46 47 // PKG_ADD: autoload ("__sl_sb10jd__", "__control_slicot_functions__.oct"); 48 DEFUN_DLD (__sl_sb10jd__, args, nargout, 49 "-*- texinfo -*-\n\ 50 Slicot SB10JD Release 5.0\n\ 51 No argument checking.\n\ 52 For internal use only.") 53 { 54 octave_idx_type nargin = args.length (); 55 octave_value_list retval; 56 57 if (nargin != 5) 58 { 59 print_usage (); 60 } 61 else 62 { 63 // arguments in 64 Matrix a = args(0).matrix_value (); 65 Matrix b = args(1).matrix_value (); 66 Matrix c = args(2).matrix_value (); 67 Matrix d = args(3).matrix_value (); 68 Matrix e = args(4).matrix_value (); 69 70 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 71 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 72 F77_INT np = TO_F77_INT (c.rows ()); // np: number of outputs 73 74 F77_INT lda = max (1, n); 75 F77_INT ldb = max (1, n); 76 F77_INT ldc = max (1, np); 77 F77_INT ldd = max (1, np); 78 F77_INT lde = max (1, n); 79 80 // arguments out 81 F77_INT nsys; 82 83 // workspace 84 F77_INT ldwork = max (1, 2*n*n + 2*n + n*max (5, n + m + np)); 85 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 86 87 // error indicator 88 F77_INT info; 89 90 91 // SLICOT routine SB10JD 92 F77_XFCN (sb10jd, SB10JD, 93 (n, m, np, 94 a.fortran_vec (), lda, 95 b.fortran_vec (), ldb, 96 c.fortran_vec (), ldc, 97 d.fortran_vec (), ldd, 98 e.fortran_vec (), lde, 99 nsys, 100 dwork, ldwork, 101 info)); 102 103 if (f77_exception_encountered) 104 error ("__sl_sb10jd__: exception in SLICOT subroutine SB10JD"); 105 106 if (info != 0) 107 error ("__sl_sb10jd__: SB10JD returned info = %d", static_cast<int> (info)); 108 109 // resize 110 a.resize (nsys, nsys); 111 b.resize (nsys, m); 112 c.resize (np, nsys); 113 114 // return values 115 retval(0) = a; 116 retval(1) = b; 117 retval(2) = c; 118 retval(3) = d; 119 } 120 121 return retval; 122 } 123