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 SLICOT system identification 21 Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: March 2012 26 Version: 0.2 27 28 */ 29 30 #include <octave/oct.h> 31 #include <octave/Cell.h> 32 #include "common.h" 33 34 extern "C" 35 { 36 int F77_FUNC (ib01ad, IB01AD) 37 (char& METH, char& ALG, char& JOBD, 38 char& BATCH, char& CONCT, char& CTRL, 39 F77_INT& NOBR, F77_INT& M, F77_INT& L, 40 F77_INT& NSMP, 41 double* U, F77_INT& LDU, 42 double* Y, F77_INT& LDY, 43 F77_INT& N, 44 double* R, F77_INT& LDR, 45 double* SV, 46 double& RCOND, double& TOL, 47 F77_INT* IWORK, 48 double* DWORK, F77_INT& LDWORK, 49 F77_INT& IWARN, F77_INT& INFO); 50 } 51 52 // PKG_ADD: autoload ("__sl_ib01ad__", "__control_slicot_functions__.oct"); 53 DEFUN_DLD (__sl_ib01ad__, args, nargout, 54 "-*- texinfo -*-\n\ 55 Slicot IB01AD Release 5.0\n\ 56 No argument checking.\n\ 57 For internal use only.") 58 { 59 octave_idx_type nargin = args.length (); 60 octave_value_list retval; 61 62 if (nargin != 10) 63 { 64 print_usage (); 65 } 66 else 67 { 68 //////////////////////////////////////////////////////////////////////////////////// 69 // SLICOT IB01AD - preprocess the input-output data // 70 //////////////////////////////////////////////////////////////////////////////////// 71 72 // arguments in 73 char meth_a; 74 char alg; 75 char jobd; 76 char batch; 77 char conct; 78 char ctrl = 'N'; 79 80 const Cell y_cell = args(0).cell_value (); 81 const Cell u_cell = args(1).cell_value (); 82 F77_INT nobr = args(2).int_value (); 83 // F77_INT nuser = args(3).int_value (); 84 85 const F77_INT imeth = args(4).int_value (); 86 const F77_INT ialg = args(5).int_value (); 87 const F77_INT iconct = args(6).int_value (); 88 // const F77_INT ictrl = args(7).int_value (); // ignored 89 90 double rcond = args(8).double_value (); 91 double tol_a = args(9).double_value (); 92 93 // double tol_b = rcond; 94 // doublet ol_c = rcond; 95 96 97 switch (imeth) 98 { 99 case 0: 100 meth_a = 'M'; 101 break; 102 case 1: 103 meth_a = 'N'; 104 break; 105 case 2: 106 meth_a = 'N'; // no typo here 107 break; 108 default: 109 error ("__sl_ib01ad__: argument 'meth' invalid"); 110 } 111 112 switch (ialg) 113 { 114 case 0: 115 alg = 'C'; 116 break; 117 case 1: 118 alg = 'F'; 119 break; 120 case 2: 121 alg = 'Q'; 122 break; 123 default: 124 error ("__sl_ib01ad__: argument 'alg' invalid"); 125 } 126 127 if (meth_a == 'M') 128 jobd = 'M'; 129 else // meth_a == 'N' 130 jobd = 'N'; // IB01AD.f says: This parameter is not relevant for METH = 'N' 131 132 if (iconct == 0) 133 conct = 'C'; 134 else 135 conct = 'N'; 136 /* 137 if (ictrl == 0) 138 ctrl = 'C'; 139 else 140 ctrl = 'N'; 141 */ 142 // m and l are equal for all experiments, checked by iddata class 143 F77_INT n_exp = TO_F77_INT (y_cell.numel ()); // number of experiments 144 F77_INT m = TO_F77_INT (u_cell.elem(0).columns ()); // m: number of inputs 145 F77_INT l = TO_F77_INT (y_cell.elem(0).columns ()); // l: number of outputs 146 F77_INT nsmpl = 0; // total number of samples 147 148 // arguments out 149 F77_INT n; 150 F77_INT ldr; 151 152 if (meth_a == 'M' && jobd == 'M') 153 ldr = max (2*(m+l)*nobr, 3*m*nobr); 154 else if (meth_a == 'N' || (meth_a == 'M' && jobd == 'N')) 155 ldr = 2*(m+l)*nobr; 156 else 157 error ("__sl_ib01ad__: could not handle 'ldr' case"); 158 159 Matrix r (ldr, 2*(m+l)*nobr); 160 ColumnVector sv (l*nobr); 161 162 163 // repeat for every experiment in the dataset 164 for (F77_INT i = 0; i < n_exp; i++) 165 { 166 if (n_exp == 1) 167 batch = 'O'; // one block only 168 else if (i == 0) 169 batch = 'F'; // first block 170 else if (i == n_exp-1) 171 batch = 'L'; // last block 172 else 173 batch = 'I'; // intermediate block 174 175 Matrix y = y_cell.elem(i).matrix_value (); 176 Matrix u = u_cell.elem(i).matrix_value (); 177 178 // y.rows == u.rows is checked by iddata class 179 // F77_INT m = TO_F77_INT (u.columns ()); // m: number of inputs 180 // F77_INT l = TO_F77_INT (y.columns ()); // l: number of outputs 181 F77_INT nsmp = TO_F77_INT (y.rows ()); // nsmp: number of samples in the current experiment 182 nsmpl += nsmp; // nsmpl: total number of samples of all experiments 183 184 // minimal nsmp size checked by __slicot_identification__.m 185 if (batch == 'O') 186 { 187 if (nsmp < 2*(m+l+1)*nobr - 1) 188 error ("__sl_ident__: require NSMP >= 2*(M+L+1)*NOBR - 1"); 189 } 190 else 191 { 192 if (nsmp < 2*nobr) 193 error ("__sl_ident__: require NSMP >= 2*NOBR"); 194 } 195 196 F77_INT ldu; 197 198 if (m == 0) 199 ldu = 1; 200 else // m > 0 201 ldu = nsmp; 202 203 F77_INT ldy = nsmp; 204 205 // workspace 206 F77_INT liwork_a; 207 208 if (meth_a == 'N') // if METH = 'N' 209 liwork_a = (m+l)*nobr; 210 else if (alg == 'F') // if METH = 'M' and ALG = 'F' 211 liwork_a = m+l; 212 else // if METH = 'M' and ALG = 'C' or 'Q' 213 liwork_a = 0; 214 215 // TODO: Handle 'k' for DWORK 216 217 F77_INT ldwork_a; 218 F77_INT ns = nsmp - 2*nobr + 1; 219 220 if (alg == 'C') 221 { 222 if (batch == 'F' || batch == 'I') 223 { 224 if (conct == 'C') 225 ldwork_a = (4*nobr-2)*(m+l); 226 else // (conct == 'N') 227 ldwork_a = 1; 228 } 229 else if (meth_a == 'M') // && (batch == 'L' || batch == 'O') 230 { 231 if (conct == 'C' && batch == 'L') 232 ldwork_a = max ((4*nobr-2)*(m+l), 5*l*nobr); 233 else if (jobd == 'M') 234 ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); 235 else // (jobd == 'N') 236 ldwork_a = 5*l*nobr; 237 } 238 else // meth_a == 'N' && (batch == 'L' || batch == 'O') 239 { 240 ldwork_a = 5*(m+l)*nobr + 1; 241 } 242 } 243 else if (alg == 'F') 244 { 245 if (batch != 'O' && conct == 'C') 246 ldwork_a = (m+l)*2*nobr*(m+l+3); 247 else if (batch == 'F' || batch == 'I') // && conct == 'N' 248 ldwork_a = (m+l)*2*nobr*(m+l+1); 249 else // (batch == 'L' || '0' && conct == 'N') 250 ldwork_a = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; 251 } 252 else // (alg == 'Q') 253 { 254 // F77_INT ns = nsmp - 2*nobr + 1; 255 256 if (ldr >= ns && batch == 'F') 257 { 258 ldwork_a = 4*(m+l)*nobr; 259 } 260 else if (ldr >= ns && batch == 'O') 261 { 262 if (meth_a == 'M') 263 ldwork_a = max (4*(m+l)*nobr, 5*l*nobr); 264 else // (meth == 'N') 265 ldwork_a = 5*(m+l)*nobr + 1; 266 } 267 else if (conct == 'C' && (batch == 'I' || batch == 'L')) 268 { 269 ldwork_a = 4*(nobr+1)*(m+l)*nobr; 270 } 271 else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') 272 { 273 ldwork_a = 6*(m+l)*nobr; 274 } 275 } 276 277 /* 278 IB01AD.f Lines 438-445 279 C FURTHER COMMENTS 280 C 281 C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the 282 C calculations could be rather inefficient if only minimal workspace 283 C (see argument LDWORK) is provided. It is advisable to provide as 284 C much workspace as possible. Almost optimal efficiency can be 285 C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the 286 C cache size is large enough to accommodate R, U, Y, and DWORK. 287 */ 288 289 ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr)); 290 291 /* 292 IB01AD.f Lines 291-195: 293 c the workspace used for alg = 'q' is 294 c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, 295 c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended 296 c value ldrwrk = ns, assuming a large enough cache size. 297 c for good performance, ldwork should be larger. 298 299 somehow ldrwrk and ldwork must have been mixed up here 300 301 */ 302 303 304 OCTAVE_LOCAL_BUFFER (F77_INT, iwork_a, liwork_a); 305 OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a); 306 307 // error indicators 308 F77_INT iwarn_a = 0; 309 F77_INT info_a = 0; 310 311 312 // SLICOT routine IB01AD 313 F77_XFCN (ib01ad, IB01AD, 314 (meth_a, alg, jobd, 315 batch, conct, ctrl, 316 nobr, m, l, 317 nsmp, 318 u.fortran_vec (), ldu, 319 y.fortran_vec (), ldy, 320 n, 321 r.fortran_vec (), ldr, 322 sv.fortran_vec (), 323 rcond, tol_a, 324 iwork_a, 325 dwork_a, ldwork_a, 326 iwarn_a, info_a)); 327 328 329 if (f77_exception_encountered) 330 error ("ident: exception in SLICOT subroutine IB01AD"); 331 332 static const char* err_msg[] = { 333 "0: OK", 334 "1: a fast algorithm was requested (ALG = 'C', or 'F') " 335 "in sequential data processing, but it failed; the " 336 "routine can be repeatedly called again using the " 337 "standard QR algorithm", 338 "2: the singular value decomposition (SVD) algorithm did " 339 "not converge"}; 340 341 static const char* warn_msg[] = { 342 "0: OK", 343 "1: the number of 100 cycles in sequential data " 344 "processing has been exhausted without signaling " 345 "that the last block of data was get; the cycle " 346 "counter was reinitialized", 347 "2: a fast algorithm was requested (ALG = 'C' or 'F'), " 348 "but it failed, and the QR algorithm was then used " 349 "(non-sequential data processing)", 350 "3: all singular values were exactly zero, hence N = 0 " 351 "(both input and output were identically zero)", 352 "4: the least squares problems with coefficient matrix " 353 "U_f, used for computing the weighted oblique " 354 "projection (for METH = 'N'), have a rank-deficient " 355 "coefficient matrix", 356 "5: the least squares problem with coefficient matrix " 357 "r_1 [6], used for computing the weighted oblique " 358 "projection (for METH = 'N'), has a rank-deficient " 359 "coefficient matrix"}; 360 361 362 error_msg ("ident: IB01AD", info_a, 2, err_msg); 363 warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg); 364 } 365 366 367 // resize 368 F77_INT rs = 2*(m+l)*nobr; 369 r.resize (rs, rs); 370 371 372 // return values 373 retval(0) = sv; 374 retval(1) = octave_value (n); 375 } 376 377 return retval; 378 } 379