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 int F77_FUNC (ib01bd, IB01BD) 52 (char& METH, char& JOB, char& JOBCK, 53 F77_INT& NOBR, F77_INT& N, F77_INT& M, F77_INT& L, 54 F77_INT& NSMPL, 55 double* R, F77_INT& LDR, 56 double* A, F77_INT& LDA, 57 double* C, F77_INT& LDC, 58 double* B, F77_INT& LDB, 59 double* D, F77_INT& LDD, 60 double* Q, F77_INT& LDQ, 61 double* RY, F77_INT& LDRY, 62 double* S, F77_INT& LDS, 63 double* K, F77_INT& LDK, 64 double& TOL, 65 F77_INT* IWORK, 66 double* DWORK, F77_INT& LDWORK, 67 F77_LOGICAL* BWORK, 68 F77_INT& IWARN, F77_INT& INFO); 69 70 int F77_FUNC (ib01cd, IB01CD) 71 (char& JOBX0, char& COMUSE, char& JOB, 72 F77_INT& N, F77_INT& M, F77_INT& L, 73 F77_INT& NSMP, 74 double* A, F77_INT& LDA, 75 double* B, F77_INT& LDB, 76 double* C, F77_INT& LDC, 77 double* D, F77_INT& LDD, 78 double* U, F77_INT& LDU, 79 double* Y, F77_INT& LDY, 80 double* X0, 81 double* V, F77_INT& LDV, 82 double& TOL, 83 F77_INT* IWORK, 84 double* DWORK, F77_INT& LDWORK, 85 F77_INT& IWARN, F77_INT& INFO); 86 } 87 88 // PKG_ADD: autoload ("__sl_ident__", "__control_slicot_functions__.oct"); 89 DEFUN_DLD (__sl_ident__, args, nargout, 90 "-*- texinfo -*-\n\ 91 Slicot IB01AD Release 5.0\n\ 92 No argument checking.\n\ 93 For internal use only.") 94 { 95 octave_idx_type nargin = args.length (); 96 octave_value_list retval; 97 98 if (nargin != 10) 99 { 100 print_usage (); 101 } 102 else 103 { 104 //////////////////////////////////////////////////////////////////////////////////// 105 // SLICOT IB01AD - preprocess the input-output data // 106 //////////////////////////////////////////////////////////////////////////////////// 107 108 // arguments in 109 char meth_a; 110 char meth_b; 111 char alg; 112 char jobd; 113 char batch; 114 char conct; 115 char ctrl; 116 117 const Cell y_cell = args(0).cell_value (); 118 const Cell u_cell = args(1).cell_value (); 119 F77_INT nobr = args(2).int_value (); 120 F77_INT nuser = args(3).int_value (); 121 122 const F77_INT imeth = args(4).int_value (); 123 const F77_INT ialg = args(5).int_value (); 124 const F77_INT iconct = args(6).int_value (); 125 const F77_INT ictrl = args(7).int_value (); 126 127 double rcond = args(8).double_value (); 128 double tol_a = args(9).double_value (); 129 130 double tol_b = rcond; 131 double tol_c = rcond; 132 133 134 switch (imeth) 135 { 136 case 0: 137 meth_a = 'M'; 138 meth_b = 'M'; 139 break; 140 case 1: 141 meth_a = 'N'; 142 meth_b = 'N'; 143 break; 144 case 2: 145 meth_a = 'N'; // no typo here 146 meth_b = 'C'; 147 break; 148 default: 149 error ("__sl_ib01ad__: argument 'meth' invalid"); 150 } 151 152 switch (ialg) 153 { 154 case 0: 155 alg = 'C'; 156 break; 157 case 1: 158 alg = 'F'; 159 break; 160 case 2: 161 alg = 'Q'; 162 break; 163 default: 164 error ("__sl_ib01ad__: argument 'alg' invalid"); 165 } 166 167 if (meth_a == 'M') 168 jobd = 'M'; 169 else // meth_a == 'N' 170 jobd = 'N'; // IB01AD.f says: This parameter is not relevant for METH = 'N' 171 172 if (iconct == 0) 173 conct = 'C'; 174 else 175 conct = 'N'; 176 177 if (ictrl == 0) 178 ctrl = 'C'; 179 else 180 ctrl = 'N'; 181 182 // m and l are equal for all experiments, checked by iddata class 183 F77_INT n_exp = TO_F77_INT (y_cell.numel ()); // number of experiments 184 F77_INT m = TO_F77_INT (u_cell.elem(0).columns ()); // m: number of inputs 185 F77_INT l = TO_F77_INT (y_cell.elem(0).columns ()); // l: number of outputs 186 F77_INT nsmpl = 0; // total number of samples 187 188 // arguments out 189 F77_INT n; 190 F77_INT ldr; 191 192 if (meth_a == 'M' && jobd == 'M') 193 ldr = max (2*(m+l)*nobr, 3*m*nobr); 194 else if (meth_a == 'N' || (meth_a == 'M' && jobd == 'N')) 195 ldr = 2*(m+l)*nobr; 196 else 197 error ("__sl_ib01ad__: could not handle 'ldr' case"); 198 199 Matrix r (ldr, 2*(m+l)*nobr); 200 ColumnVector sv (l*nobr); 201 202 203 // repeat for every experiment in the dataset 204 for (F77_INT i = 0; i < n_exp; i++) 205 { 206 if (n_exp == 1) 207 batch = 'O'; // one block only 208 else if (i == 0) 209 batch = 'F'; // first block 210 else if (i == n_exp-1) 211 batch = 'L'; // last block 212 else 213 batch = 'I'; // intermediate block 214 215 Matrix y = y_cell.elem(i).matrix_value (); 216 Matrix u = u_cell.elem(i).matrix_value (); 217 218 // y.rows == u.rows is checked by iddata class 219 // F77_INT m = TO_F77_INT (u.columns ()); // m: number of inputs 220 // F77_INT l = TO_F77_INT (y.columns ()); // l: number of outputs 221 F77_INT nsmp = TO_F77_INT (y.rows ()); // nsmp: number of samples in the current experiment 222 nsmpl += nsmp; // nsmpl: total number of samples of all experiments 223 224 // minimal nsmp size checked by __slicot_identification__.m 225 if (batch == 'O') 226 { 227 if (nsmp < 2*(m+l+1)*nobr - 1) 228 error ("__sl_ident__: require NSMP >= 2*(M+L+1)*NOBR - 1"); 229 } 230 else 231 { 232 if (nsmp < 2*nobr) 233 error ("__sl_ident__: require NSMP >= 2*NOBR"); 234 } 235 236 F77_INT ldu; 237 238 if (m == 0) 239 ldu = 1; 240 else // m > 0 241 ldu = nsmp; 242 243 F77_INT ldy = nsmp; 244 245 // workspace 246 F77_INT liwork_a; 247 248 if (meth_a == 'N') // if METH = 'N' 249 liwork_a = (m+l)*nobr; 250 else if (alg == 'F') // if METH = 'M' and ALG = 'F' 251 liwork_a = m+l; 252 else // if METH = 'M' and ALG = 'C' or 'Q' 253 liwork_a = 0; 254 255 // TODO: Handle 'k' for DWORK 256 257 F77_INT ldwork_a; 258 F77_INT ns = nsmp - 2*nobr + 1; 259 260 if (alg == 'C') 261 { 262 if (batch == 'F' || batch == 'I') 263 { 264 if (conct == 'C') 265 ldwork_a = (4*nobr-2)*(m+l); 266 else // (conct == 'N') 267 ldwork_a = 1; 268 } 269 else if (meth_a == 'M') // && (batch == 'L' || batch == 'O') 270 { 271 if (conct == 'C' && batch == 'L') 272 ldwork_a = max ((4*nobr-2)*(m+l), 5*l*nobr); 273 else if (jobd == 'M') 274 ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr); 275 else // (jobd == 'N') 276 ldwork_a = 5*l*nobr; 277 } 278 else // meth_b == 'N' && (batch == 'L' || batch == 'O') 279 { 280 ldwork_a = 5*(m+l)*nobr + 1; 281 } 282 } 283 else if (alg == 'F') 284 { 285 if (batch != 'O' && conct == 'C') 286 ldwork_a = (m+l)*2*nobr*(m+l+3); 287 else if (batch == 'F' || batch == 'I') // && conct == 'N' 288 ldwork_a = (m+l)*2*nobr*(m+l+1); 289 else // (batch == 'L' || '0' && conct == 'N') 290 ldwork_a = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr; 291 } 292 else // (alg == 'Q') 293 { 294 // F77_INT ns = nsmp - 2*nobr + 1; 295 296 if (ldr >= ns && batch == 'F') 297 { 298 ldwork_a = 4*(m+l)*nobr; 299 } 300 else if (ldr >= ns && batch == 'O') 301 { 302 if (meth_a == 'M') 303 ldwork_a = max (4*(m+l)*nobr, 5*l*nobr); 304 else // (meth == 'N') 305 ldwork_a = 5*(m+l)*nobr + 1; 306 } 307 else if (conct == 'C' && (batch == 'I' || batch == 'L')) 308 { 309 ldwork_a = 4*(nobr+1)*(m+l)*nobr; 310 } 311 else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N') 312 { 313 ldwork_a = 6*(m+l)*nobr; 314 } 315 } 316 317 /* 318 IB01AD.f Lines 438-445 319 C FURTHER COMMENTS 320 C 321 C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the 322 C calculations could be rather inefficient if only minimal workspace 323 C (see argument LDWORK) is provided. It is advisable to provide as 324 C much workspace as possible. Almost optimal efficiency can be 325 C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the 326 C cache size is large enough to accommodate R, U, Y, and DWORK. 327 */ 328 329 ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr)); 330 331 /* 332 IB01AD.f Lines 291-195: 333 c the workspace used for alg = 'q' is 334 c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr, 335 c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended 336 c value ldrwrk = ns, assuming a large enough cache size. 337 c for good performance, ldwork should be larger. 338 339 somehow ldrwrk and ldwork must have been mixed up here 340 341 */ 342 343 344 OCTAVE_LOCAL_BUFFER (F77_INT, iwork_a, liwork_a); 345 OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a); 346 347 // error indicators 348 F77_INT iwarn_a = 0; 349 F77_INT info_a = 0; 350 351 352 // SLICOT routine IB01AD 353 F77_XFCN (ib01ad, IB01AD, 354 (meth_a, alg, jobd, 355 batch, conct, ctrl, 356 nobr, m, l, 357 nsmp, 358 u.fortran_vec (), ldu, 359 y.fortran_vec (), ldy, 360 n, 361 r.fortran_vec (), ldr, 362 sv.fortran_vec (), 363 rcond, tol_a, 364 iwork_a, 365 dwork_a, ldwork_a, 366 iwarn_a, info_a)); 367 368 369 if (f77_exception_encountered) 370 error ("ident: exception in SLICOT subroutine IB01AD"); 371 372 static const char* err_msg[] = { 373 "0: OK", 374 "1: a fast algorithm was requested (ALG = 'C', or 'F') " 375 "in sequential data processing, but it failed; the " 376 "routine can be repeatedly called again using the " 377 "standard QR algorithm", 378 "2: the singular value decomposition (SVD) algorithm did " 379 "not converge"}; 380 381 static const char* warn_msg[] = { 382 "0: OK", 383 "1: the number of 100 cycles in sequential data " 384 "processing has been exhausted without signaling " 385 "that the last block of data was get; the cycle " 386 "counter was reinitialized", 387 "2: a fast algorithm was requested (ALG = 'C' or 'F'), " 388 "but it failed, and the QR algorithm was then used " 389 "(non-sequential data processing)", 390 "3: all singular values were exactly zero, hence N = 0 " 391 "(both input and output were identically zero)", 392 "4: the least squares problems with coefficient matrix " 393 "U_f, used for computing the weighted oblique " 394 "projection (for METH = 'N'), have a rank-deficient " 395 "coefficient matrix", 396 "5: the least squares problem with coefficient matrix " 397 "r_1 [6], used for computing the weighted oblique " 398 "projection (for METH = 'N'), has a rank-deficient " 399 "coefficient matrix"}; 400 401 402 error_msg ("ident: IB01AD", info_a, 2, err_msg); 403 warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg); 404 } 405 406 407 // resize 408 F77_INT rs = 2*(m+l)*nobr; 409 r.resize (rs, rs); 410 411 if (nuser > 0) 412 { 413 if (nuser < nobr) 414 { 415 n = nuser; 416 // warning ("ident: nuser (%d) < nobr (%d), n = nuser", nuser, nobr); 417 } 418 else 419 error ("ident: 'nuser' invalid"); 420 } 421 422 //////////////////////////////////////////////////////////////////////////////////// 423 // SLICOT IB01BD - estimating system matrices, Kalman gain, and covariances // 424 //////////////////////////////////////////////////////////////////////////////////// 425 426 // arguments in 427 char job = 'A'; 428 char jobck = 'K'; 429 430 //F77_INT nsmpl = nsmp; 431 432 if (nsmpl < 2*(m+l)*nobr) 433 error ("__sl_ident__: nsmpl (%d) < 2*(m+l)*nobr (%d)", 434 static_cast<int> (nsmpl), static_cast<int> (nobr)); 435 436 // arguments out 437 F77_INT lda = max (1, n); 438 F77_INT ldc = max (1, l); 439 F77_INT ldb = max (1, n); 440 F77_INT ldd = max (1, l); 441 F77_INT ldq = n; // if JOBCK = 'C' or 'K' 442 F77_INT ldry = l; // if JOBCK = 'C' or 'K' 443 F77_INT lds = n; // if JOBCK = 'C' or 'K' 444 F77_INT ldk = n; // if JOBCK = 'K' 445 446 Matrix a (lda, n); 447 Matrix c (ldc, n); 448 Matrix b (ldb, m); 449 Matrix d (ldd, m); 450 451 Matrix q (ldq, n); 452 Matrix ry (ldry, l); 453 Matrix s (lds, l); 454 Matrix k (ldk, l); 455 456 // workspace 457 F77_INT liwork_b; 458 F77_INT liw1; 459 F77_INT liw2; 460 461 liw1 = max (n, m*nobr+n, l*nobr, m*(n+l)); 462 liw2 = n*n; // if JOBCK = 'K' 463 liwork_b = max (liw1, liw2); 464 465 F77_INT ldwork_b; 466 F77_INT ldw1; 467 F77_INT ldw2; 468 F77_INT ldw3; 469 470 if (meth_b == 'M') 471 { 472 F77_INT ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); 473 F77_INT ldw1b = max (2*(l*nobr-l)*n+n*n+7*n, 474 (l*nobr-l)*n+n+6*m*nobr, 475 (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m))); 476 ldw1 = max (ldw1a, ldw1b); 477 478 F77_INT aw; 479 480 if (m == 0 || job == 'C') 481 aw = n + n*n; 482 else 483 aw = 0; 484 485 ldw2 = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l ); 486 } 487 else if (meth_b == 'N') 488 { 489 ldw1 = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, 490 2*(l*nobr-l)*n+n*n+8*n, 491 n+4*(m*nobr+n)+1, 492 m*nobr+3*n+l); 493 494 if (m == 0 || job == 'C') 495 ldw2 = 0; 496 else 497 ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); 498 499 } 500 else // (meth_b == 'C') 501 { 502 F77_INT ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n); 503 F77_INT ldw1b = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l, 504 2*(l*nobr-l)*n+n*n+8*n, 505 n+4*(m*nobr+n)+1, 506 m*nobr+3*n+l); 507 508 ldw1 = max (ldw1a, ldw1b); 509 510 ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1); 511 512 } 513 514 ldw3 = max(4*n*n + 2*n*l + l*l + max (3*l, n*l), 14*n*n + 12*n + 5); 515 ldwork_b = max (ldw1, ldw2, ldw3); 516 517 518 OCTAVE_LOCAL_BUFFER (F77_INT, iwork_b, liwork_b); 519 OCTAVE_LOCAL_BUFFER (double, dwork_b, ldwork_b); 520 OCTAVE_LOCAL_BUFFER (F77_LOGICAL, bwork, 2*n); 521 522 523 // error indicators 524 F77_INT iwarn_b = 0; 525 F77_INT info_b = 0; 526 527 528 // SLICOT routine IB01BD 529 F77_XFCN (ib01bd, IB01BD, 530 (meth_b, job, jobck, 531 nobr, n, m, l, 532 nsmpl, 533 r.fortran_vec (), ldr, 534 a.fortran_vec (), lda, 535 c.fortran_vec (), ldc, 536 b.fortran_vec (), ldb, 537 d.fortran_vec (), ldd, 538 q.fortran_vec (), ldq, 539 ry.fortran_vec (), ldry, 540 s.fortran_vec (), lds, 541 k.fortran_vec (), ldk, 542 tol_b, 543 iwork_b, 544 dwork_b, ldwork_b, 545 bwork, 546 iwarn_b, info_b)); 547 548 549 if (f77_exception_encountered) 550 error ("ident: exception in SLICOT subroutine IB01BD"); 551 552 static const char* err_msg_b[] = { 553 "0: OK", 554 "1: error message not specified", 555 "2: the singular value decomposition (SVD) algorithm did " 556 "not converge", 557 "3: a singular upper triangular matrix was found", 558 "4: matrix A is (numerically) singular in discrete-" 559 "time case", 560 "5: the Hamiltonian or symplectic matrix H cannot be " 561 "reduced to real Schur form", 562 "6: the real Schur form of the Hamiltonian or " 563 "symplectic matrix H cannot be appropriately ordered", 564 "7: the Hamiltonian or symplectic matrix H has less " 565 "than N stable eigenvalues", 566 "8: the N-th order system of linear algebraic " 567 "equations, from which the solution matrix X would " 568 "be obtained, is singular to working precision", 569 "9: the QR algorithm failed to complete the reduction " 570 "of the matrix Ac to Schur canonical form, T", 571 "10: the QR algorithm did not converge"}; 572 573 static const char* warn_msg_b[] = { 574 "0: OK", 575 "1: warning message not specified", 576 "2: warning message not specified", 577 "3: warning message not specified", 578 "4: a least squares problem to be solved has a " 579 "rank-deficient coefficient matrix", 580 "5: the computed covariance matrices are too small. " 581 "The problem seems to be a deterministic one; the " 582 "gain matrix is set to zero"}; 583 584 585 error_msg ("ident: IB01BD", info_b, 10, err_msg_b); 586 warning_msg ("ident: IB01BD", iwarn_b, 5, warn_msg_b); 587 588 // resize 589 a.resize (n, n); 590 c.resize (l, n); 591 b.resize (n, m); 592 d.resize (l, m); 593 594 q.resize (n, n); 595 ry.resize (l, l); 596 s.resize (n, l); 597 k.resize (n, l); 598 599 //////////////////////////////////////////////////////////////////////////////////// 600 // SLICOT IB01CD - estimating the initial state // 601 //////////////////////////////////////////////////////////////////////////////////// 602 603 // arguments in 604 char jobx0 = 'X'; 605 char comuse = 'U'; 606 char jobbd = 'D'; 607 608 // arguments out 609 Cell x0_cell (n_exp, 1); // cell of initial state vectors x0 610 611 // repeat for every experiment in the dataset 612 // compute individual initial state vector x0 for every experiment 613 for (F77_INT i = 0; i < n_exp; i++) 614 { 615 Matrix y = y_cell.elem(i).matrix_value (); 616 Matrix u = u_cell.elem(i).matrix_value (); 617 618 F77_INT nsmp = TO_F77_INT (y.rows ()); // nsmp: number of samples 619 F77_INT ldv = max (1, n); 620 621 F77_INT ldu; 622 623 if (m == 0) 624 ldu = 1; 625 else // m > 0 626 ldu = nsmp; 627 628 F77_INT ldy = nsmp; 629 630 // arguments out 631 ColumnVector x0 (n); 632 Matrix v (ldv, n); 633 634 // workspace 635 F77_INT liwork_c = n; // if JOBX0 = 'X' and COMUSE <> 'C' 636 F77_INT ldwork_c; 637 F77_INT t = nsmp; 638 639 F77_INT ldw1_c = 2; 640 F77_INT ldw2_c = t*l*(n + 1) + 2*n + max (2*n*n, 4*n); 641 F77_INT ldw3_c = n*(n + 1) + 2*n + max (n*l*(n + 1) + 2*n*n + l*n, 4*n); 642 643 ldwork_c = ldw1_c + n*( n + m + l ) + max (5*n, ldw1_c, min (ldw2_c, ldw3_c)); 644 645 OCTAVE_LOCAL_BUFFER (F77_INT, iwork_c, liwork_c); 646 OCTAVE_LOCAL_BUFFER (double, dwork_c, ldwork_c); 647 648 // error indicators 649 F77_INT iwarn_c = 0; 650 F77_INT info_c = 0; 651 652 // SLICOT routine IB01CD 653 F77_XFCN (ib01cd, IB01CD, 654 (jobx0, comuse, jobbd, 655 n, m, l, 656 nsmp, 657 a.fortran_vec (), lda, 658 b.fortran_vec (), ldb, 659 c.fortran_vec (), ldc, 660 d.fortran_vec (), ldd, 661 u.fortran_vec (), ldu, 662 y.fortran_vec (), ldy, 663 x0.fortran_vec (), 664 v.fortran_vec (), ldv, 665 tol_c, 666 iwork_c, 667 dwork_c, ldwork_c, 668 iwarn_c, info_c)); 669 670 671 if (f77_exception_encountered) 672 error ("ident: exception in SLICOT subroutine IB01CD"); 673 674 static const char* err_msg_c[] = { 675 "0: OK", 676 "1: the QR algorithm failed to compute all the " 677 "eigenvalues of the matrix A (see LAPACK Library " 678 "routine DGEES); the locations DWORK(i), for " 679 "i = g+1:g+N*N, contain the partially converged " 680 "Schur form", 681 "2: the singular value decomposition (SVD) algorithm did " 682 "not converge"}; 683 684 static const char* warn_msg_c[] = { 685 "0: OK", 686 "1: warning message not specified", 687 "2: warning message not specified", 688 "3: warning message not specified", 689 "4: the least squares problem to be solved has a " 690 "rank-deficient coefficient matrix", 691 "5: warning message not specified", 692 "6: the matrix A is unstable; the estimated x(0) " 693 "and/or B and D could be inaccurate"}; 694 695 696 error_msg ("ident: IB01CD", info_c, 2, err_msg_c); 697 warning_msg ("ident: IB01CD", iwarn_c, 6, warn_msg_c); 698 699 x0_cell.elem(i) = x0; // add x0 from the current experiment to cell of initial state vectors 700 } 701 702 703 // return values 704 retval(0) = a; 705 retval(1) = b; 706 retval(2) = c; 707 retval(3) = d; 708 709 retval(4) = q; 710 retval(5) = ry; 711 retval(6) = s; 712 retval(7) = k; 713 714 retval(8) = x0_cell; 715 } 716 717 return retval; 718 } 719