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 controllability form. 21 Uses SLICOT AB01OD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: August 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 (ab01od, AB01OD) 36 (char& STAGES, 37 char& JOBU, char& JOBV, 38 F77_INT& N, F77_INT& M, 39 double* A, F77_INT& LDA, 40 double* B, F77_INT& LDB, 41 double* U, F77_INT& LDU, 42 double* V, F77_INT& LDV, 43 F77_INT& NCONT, F77_INT& INDCON, 44 F77_INT* KSTAIR, 45 double& TOL, 46 F77_INT* IWORK, 47 double* DWORK, F77_INT& LDWORK, 48 F77_INT& INFO); 49 } 50 51 // PKG_ADD: autoload ("__sl_ab01od__", "__control_slicot_functions__.oct"); 52 DEFUN_DLD (__sl_ab01od__, args, nargout, 53 "-*- texinfo -*-\n\ 54 Slicot AB01OD Release 5.0\n\ 55 No argument checking.\n\ 56 For internal use only.") 57 { 58 octave_idx_type nargin = args.length (); 59 octave_value_list retval; 60 61 if (nargin != 3) 62 { 63 print_usage (); 64 } 65 else 66 { 67 // arguments in 68 char stages = 'F'; 69 char jobu = 'I'; 70 char jobv = 'N'; // not referenced because stages = 'F' 71 72 Matrix a = args(0).matrix_value (); 73 Matrix b = args(1).matrix_value (); 74 double tol = args(2).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 79 F77_INT lda = max (1, n); 80 F77_INT ldb = max (1, n); 81 F77_INT ldu = max (1, n); 82 F77_INT ldv = 1; 83 84 // arguments out 85 Matrix u (ldu, n); 86 double* v = 0; // not referenced because stages = 'F' 87 88 F77_INT ncont; 89 F77_INT indcon; 90 91 OCTAVE_LOCAL_BUFFER (F77_INT, kstair, n); 92 93 // workspace 94 F77_INT ldwork = max (1, n + max (n, 3*m)); 95 96 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, m); 97 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 98 99 // error indicators 100 F77_INT info; 101 102 103 // SLICOT routine AB01OD 104 F77_XFCN (ab01od, AB01OD, 105 (stages, 106 jobu, jobv, 107 n, m, 108 a.fortran_vec (), lda, 109 b.fortran_vec (), ldb, 110 u.fortran_vec (), ldu, 111 v, ldv, 112 ncont, indcon, 113 kstair, 114 tol, 115 iwork, 116 dwork, ldwork, 117 info)); 118 119 if (f77_exception_encountered) 120 error ("__sl_ab01od__: exception in SLICOT subroutine AB01OD"); 121 122 if (info != 0) 123 error ("__sl_ab01od__: AB01OD returned info = %d", static_cast<int> (info)); 124 125 // resize 126 a.resize (n, n); 127 b.resize (n, m); 128 u.resize (n, n); 129 130 // return values 131 retval(0) = a; 132 retval(1) = b; 133 retval(2) = u; 134 retval(3) = octave_value (ncont); 135 } 136 137 return retval; 138 } 139