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 SB16BD 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 (sb16bd, SB16BD) 36 (char& DICO, char& JOBD, char& JOBMR, char& JOBCF, 37 char& EQUIL, 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* DC, F77_INT& LDDC, 47 double* HSV, 48 double& TOL1, double& TOL2, 49 F77_INT* IWORK, 50 double* DWORK, F77_INT& LDWORK, 51 F77_INT& IWARN, F77_INT& INFO); 52 } 53 54 // PKG_ADD: autoload ("__sl_sb16bd__", "__control_slicot_functions__.oct"); 55 DEFUN_DLD (__sl_sb16bd__, args, nargout, 56 "-*- texinfo -*-\n\ 57 Slicot SB16BD Release 5.0\n\ 58 No argument checking.\n\ 59 For internal use only.") 60 { 61 octave_idx_type nargin = args.length (); 62 octave_value_list retval; 63 64 if (nargin != 15) 65 { 66 print_usage (); 67 } 68 else 69 { 70 // arguments in 71 char dico; 72 char jobd; 73 char jobmr; 74 char jobcf; 75 char equil; 76 char ordsel; 77 78 Matrix a = args(0).matrix_value (); 79 Matrix b = args(1).matrix_value (); 80 Matrix c = args(2).matrix_value (); 81 Matrix d = args(3).matrix_value (); 82 83 const F77_INT idico = args(4).int_value (); 84 const F77_INT iequil = args(5).int_value (); 85 F77_INT ncr = args(6).int_value (); 86 const F77_INT iordsel = args(7).int_value (); 87 const F77_INT ijobd = args(8).int_value (); 88 const F77_INT ijobmr = args(9).int_value (); 89 90 Matrix f = args(10).matrix_value (); 91 Matrix g = args(11).matrix_value (); 92 93 const F77_INT ijobcf = args(12).int_value (); 94 95 double tol1 = args(13).double_value (); 96 double tol2 = args(14).double_value (); 97 98 if (idico == 0) 99 dico = 'C'; 100 else 101 dico = 'D'; 102 103 if (iequil == 0) 104 equil = 'S'; 105 else 106 equil = 'N'; 107 108 if (iordsel == 0) 109 ordsel = 'F'; 110 else 111 ordsel = 'A'; 112 113 if (ijobd == 0) 114 jobd = 'Z'; 115 else 116 jobd = 'D'; 117 118 if (ijobcf == 0) 119 jobcf = 'L'; 120 else 121 jobcf = 'R'; 122 123 switch (ijobmr) 124 { 125 case 0: 126 jobmr = 'B'; 127 break; 128 case 1: 129 jobmr = 'F'; 130 break; 131 case 2: 132 jobmr = 'S'; 133 break; 134 case 3: 135 jobmr = 'P'; 136 break; 137 default: 138 error ("__sl_sb16bd__: argument jobmr invalid"); 139 } 140 141 142 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 143 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 144 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 145 146 F77_INT lda = max (1, n); 147 F77_INT ldb = max (1, n); 148 F77_INT ldc = max (1, p); 149 F77_INT ldd; 150 151 if (jobd == 'Z') 152 ldd = 1; 153 else 154 ldd = max (1, p); 155 156 F77_INT ldf = max (1, m); 157 F77_INT ldg = max (1, n); 158 F77_INT lddc = max (1, m); 159 160 // arguments out 161 Matrix dc (lddc, p); 162 ColumnVector hsv (n); 163 164 // workspace 165 F77_INT liwork; 166 F77_INT pm; 167 168 F77_INT ldwork; 169 F77_INT lwr = max (1, n*(2*n+max(n,m+p)+5)+n*(n+1)/2); 170 171 switch (jobmr) 172 { 173 case 'B': 174 pm = 0; 175 break; 176 case 'F': 177 pm = n; 178 break; 179 default: // if JOBMR = 'S' or 'P' 180 pm = max (1, 2*n); 181 } 182 183 if (ordsel == 'F' && ncr == n) 184 { 185 liwork = 0; 186 ldwork = p*n; 187 } 188 else if (jobcf == 'L') 189 { 190 liwork = max (pm, m); 191 ldwork = (n+m)*(m+p) + max (lwr, 4*m); 192 } 193 else // if JOBCF = 'R' 194 { 195 liwork = max (pm, p); 196 ldwork = (n+p)*(m+p) + max (lwr, 4*p); 197 } 198 199 200 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 201 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 202 203 // error indicators 204 F77_INT iwarn = 0; 205 F77_INT info = 0; 206 207 208 // SLICOT routine SB16BD 209 F77_XFCN (sb16bd, SB16BD, 210 (dico, jobd, jobmr, jobcf, 211 equil, ordsel, 212 n, m, p, 213 ncr, 214 a.fortran_vec (), lda, 215 b.fortran_vec (), ldb, 216 c.fortran_vec (), ldc, 217 d.fortran_vec (), ldd, 218 f.fortran_vec (), ldf, 219 g.fortran_vec (), ldg, 220 dc.fortran_vec (), lddc, 221 hsv.fortran_vec (), 222 tol1, tol2, 223 iwork, 224 dwork, ldwork, 225 iwarn, info)); 226 227 228 if (f77_exception_encountered) 229 error ("cfconred: exception in SLICOT subroutine SB16BD"); 230 231 static const char* err_msg[] = { 232 "0: OK", 233 "1: the reduction of A-L*C to a real Schur form " 234 "failed", 235 "2: the matrix A-L*C is not stable (if DICO = 'C'), " 236 "or not convergent (if DICO = 'D')", 237 "3: the computation of Hankel singular values failed", 238 "4: the reduction of A-B*F to a real Schur form " 239 "failed", 240 "5: the matrix A-B*F is not stable (if DICO = 'C'), " 241 "or not convergent (if DICO = 'D')"}; 242 243 static const char* warn_msg[] = { 244 "0: OK", 245 "1: with ORDSEL = 'F', the selected order NCR is " 246 "greater than the order of a minimal " 247 "realization of the controller."}; 248 249 error_msg ("cfconred", info, 5, err_msg); 250 warning_msg ("cfconred", iwarn, 1, warn_msg); 251 252 253 // resize 254 a.resize (ncr, ncr); // Ac 255 g.resize (ncr, p); // Bc 256 f.resize (m, ncr); // Cc 257 dc.resize (m, p); // Dc 258 259 // return values 260 retval(0) = a; 261 retval(1) = g; 262 retval(2) = f; 263 retval(3) = dc; 264 retval(4) = octave_value (ncr); 265 retval(5) = hsv; 266 } 267 268 return retval; 269 } 270