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 TODO 21 Uses SLICOT SB16CD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: November 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 (sb16cd, SB16CD) 36 (char& DICO, char& JOBD, char& JOBMR, char& JOBCF, 37 char& ORDSEL, 38 F77_INT& N, F77_INT& M, F77_INT& P, 39 F77_INT& NCR, 40 double* A, F77_INT& LDA, 41 double* B, F77_INT& LDB, 42 double* C, F77_INT& LDC, 43 double* D, F77_INT& LDD, 44 double* F, F77_INT& LDF, 45 double* G, F77_INT& LDG, 46 double* HSV, 47 double& TOL, 48 F77_INT* IWORK, 49 double* DWORK, F77_INT& LDWORK, 50 F77_INT& IWARN, F77_INT& INFO); 51 } 52 53 // PKG_ADD: autoload ("__sl_sb16cd__", "__control_slicot_functions__.oct"); 54 DEFUN_DLD (__sl_sb16cd__, args, nargout, 55 "-*- texinfo -*-\n\ 56 Slicot SB16CD Release 5.0\n\ 57 No argument checking.\n\ 58 For internal use only.") 59 { 60 octave_idx_type nargin = args.length (); 61 octave_value_list retval; 62 63 if (nargin != 13) 64 { 65 print_usage (); 66 } 67 else 68 { 69 // arguments in 70 char dico; 71 char jobd; 72 char jobmr; 73 char jobcf; 74 char ordsel; 75 76 Matrix a = args(0).matrix_value (); 77 Matrix b = args(1).matrix_value (); 78 Matrix c = args(2).matrix_value (); 79 Matrix d = args(3).matrix_value (); 80 81 const F77_INT idico = args(4).int_value (); 82 F77_INT ncr = args(5).int_value (); 83 const F77_INT iordsel = args(6).int_value (); 84 const F77_INT ijobd = args(7).int_value (); 85 const F77_INT ijobmr = args(8).int_value (); 86 87 Matrix f = args(9).matrix_value (); 88 Matrix g = args(10).matrix_value (); 89 90 const F77_INT ijobcf = args(11).int_value (); 91 double tol = args(12).double_value (); 92 93 if (idico == 0) 94 dico = 'C'; 95 else 96 dico = 'D'; 97 98 if (iordsel == 0) 99 ordsel = 'F'; 100 else 101 ordsel = 'A'; 102 103 if (ijobd == 0) 104 jobd = 'Z'; 105 else 106 jobd = 'D'; 107 108 if (ijobcf == 0) 109 jobcf = 'L'; 110 else 111 jobcf = 'R'; 112 113 switch (ijobmr) 114 { 115 case 0: 116 jobmr = 'B'; 117 break; 118 case 1: 119 jobmr = 'F'; 120 break; 121 default: 122 error ("__sl_sb16cd__: argument jobmr invalid"); 123 } 124 125 126 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 127 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 128 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 129 130 F77_INT lda = max (1, n); 131 F77_INT ldb = max (1, n); 132 F77_INT ldc = max (1, p); 133 F77_INT ldd; 134 135 if (jobd == 'Z') 136 ldd = 1; 137 else 138 ldd = max (1, p); 139 140 F77_INT ldf = max (1, m); 141 F77_INT ldg = max (1, n); 142 143 // arguments out 144 ColumnVector hsv (n); 145 146 // workspace 147 F77_INT liwork; 148 149 if (jobmr == 'B') 150 liwork = 0; 151 else // if JOBMR = 'F' 152 liwork = n; 153 154 F77_INT ldwork; 155 F77_INT mp; 156 157 if (jobcf == 'L') 158 mp = m; 159 else // if JOBCF = 'R' 160 mp = p; 161 162 ldwork = 2*n*n + max (1, 2*n*n + 5*n, n*max(m,p), 163 n*(n + max(n,mp) + min(n,mp) + 6)); 164 165 166 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 167 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 168 169 // error indicators 170 F77_INT iwarn = 0; 171 F77_INT info = 0; 172 173 174 // SLICOT routine SB16CD 175 F77_XFCN (sb16cd, SB16CD, 176 (dico, jobd, jobmr, jobcf, 177 ordsel, 178 n, m, p, 179 ncr, 180 a.fortran_vec (), lda, 181 b.fortran_vec (), ldb, 182 c.fortran_vec (), ldc, 183 d.fortran_vec (), ldd, 184 f.fortran_vec (), ldf, 185 g.fortran_vec (), ldg, 186 hsv.fortran_vec (), 187 tol, 188 iwork, 189 dwork, ldwork, 190 iwarn, info)); 191 192 193 if (f77_exception_encountered) 194 error ("fwcfconred: exception in SLICOT subroutine SB16CD"); 195 196 static const char* err_msg[] = { 197 "0: OK", 198 "1: eigenvalue computation failure", 199 "2: the matrix A-L*C is not stable", 200 "3: the matrix A-B*F is not stable", 201 "4: the Lyapunov equation for computing the " 202 "observability Grammian is (nearly) singular", 203 "5: the Lyapunov equation for computing the " 204 "controllability Grammian is (nearly) singular", 205 "6: the computation of Hankel singular values failed"}; 206 207 static const char* warn_msg[] = { 208 "0: OK", 209 "1: with ORDSEL = 'F', the selected order NCR is " 210 "greater than the order of a minimal realization " 211 "of the controller.", 212 "2: with ORDSEL = 'F', the selected order NCR " 213 "corresponds to repeated singular values, which are " 214 "neither all included nor all excluded from the " 215 "reduced controller. In this case, the resulting NCR " 216 "is set automatically to the largest value such that " 217 "HSV(NCR) > HSV(NCR+1)."}; 218 219 error_msg ("fwcfconred", info, 6, err_msg); 220 warning_msg ("fwcfconred", iwarn, 2, warn_msg); 221 222 223 // resize 224 a.resize (ncr, ncr); // Ac 225 g.resize (ncr, p); // Bc 226 f.resize (m, ncr); // Cc 227 // Dc = 0 228 229 // return values 230 retval(0) = a; 231 retval(1) = g; 232 retval(2) = f; 233 retval(3) = octave_value (ncr); 234 retval(4) = hsv; 235 } 236 237 return retval; 238 } 239