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 L-infinity norm of a SS model. 21 Uses SLICOT AB13DD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: November 2009 26 Version: 0.5 27 28 */ 29 30 #include <octave/oct.h> 31 #include <complex> 32 #include "common.h" 33 34 extern "C" 35 { 36 int F77_FUNC (ab13dd, AB13DD) 37 (char& DICO, char& JOBE, 38 char& EQUIL, char& JOBD, 39 F77_INT& N, F77_INT& M, F77_INT& P, 40 double* FPEAK, 41 double* A, F77_INT& LDA, 42 double* E, F77_INT& LDE, 43 double* B, F77_INT& LDB, 44 double* C, F77_INT& LDC, 45 double* D, F77_INT& LDD, 46 double* GPEAK, 47 double& TOL, 48 F77_INT* IWORK, double* DWORK, F77_INT& LDWORK, 49 Complex* CWORK, F77_INT& LCWORK, 50 F77_INT& INFO); 51 } 52 53 // PKG_ADD: autoload ("__sl_ab13dd__", "__control_slicot_functions__.oct"); 54 DEFUN_DLD (__sl_ab13dd__, args, nargout, 55 "-*- texinfo -*-\n\ 56 Slicot AB13DD 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 != 9) 64 { 65 print_usage (); 66 } 67 else 68 { 69 // arguments in 70 char dico; 71 char jobe; 72 char equil; 73 char jobd = 'D'; 74 75 Matrix a = args(0).matrix_value (); 76 Matrix e = args(1).matrix_value (); 77 Matrix b = args(2).matrix_value (); 78 Matrix c = args(3).matrix_value (); 79 Matrix d = args(4).matrix_value (); 80 F77_INT discrete = args(5).int_value (); 81 F77_INT descriptor = args(6).int_value (); 82 double tol = args(7).double_value (); 83 const F77_INT scaled = args(8).int_value (); 84 85 if (discrete == 0) 86 dico = 'C'; 87 else 88 dico = 'D'; 89 90 if (descriptor == 0) 91 jobe = 'I'; 92 else 93 jobe = 'G'; 94 95 if (scaled == 0) 96 equil = 'S'; 97 else 98 equil = 'N'; 99 100 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 101 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 102 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 103 104 F77_INT lda = max (1, n); 105 F77_INT lde = max (1, n); 106 F77_INT ldb = max (1, n); 107 F77_INT ldc = max (1, p); 108 F77_INT ldd = max (1, p); 109 110 ColumnVector fpeak (2); 111 ColumnVector gpeak (2); 112 113 fpeak(0) = 0; 114 fpeak(1) = 1; 115 116 // workspace 117 F77_INT ldwork = max (1, 15*n*n + p*p + m*m + (6*n+3)*(p+m) + 4*p*m + 118 n*m + 22*n + 7*min(p,m)); 119 F77_INT lcwork = max (1, (n+m)*(n+p) + 2*min(p,m) + max(p,m)); 120 121 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, n); 122 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 123 OCTAVE_LOCAL_BUFFER (Complex, cwork, lcwork); 124 125 // error indicator 126 F77_INT info; 127 128 129 // SLICOT routine AB13DD 130 F77_XFCN (ab13dd, AB13DD, 131 (dico, jobe, 132 equil, jobd, 133 n, m, p, 134 fpeak.fortran_vec (), 135 a.fortran_vec (), lda, 136 e.fortran_vec (), lde, 137 b.fortran_vec (), ldb, 138 c.fortran_vec (), ldc, 139 d.fortran_vec (), ldd, 140 gpeak.fortran_vec (), 141 tol, 142 iwork, dwork, ldwork, 143 cwork, lcwork, 144 info)); 145 146 if (f77_exception_encountered) 147 error ("lti: norm: __sl_ab13dd__: exception in SLICOT subroutine AB13DD"); 148 149 static const char* err_msg[] = { 150 "0: OK", 151 "1: the matrix E is (numerically) singular", 152 "2: the (periodic) QR (or QZ) algorithm for computing " 153 "eigenvalues did not converge", 154 "3: the SVD algorithm for computing singular values did " 155 "not converge", 156 "4: the tolerance is too small and the algorithm did " 157 "not converge"}; 158 159 error_msg ("__sl_ab13dd__", info, 4, err_msg); 160 161 162 // return values 163 retval(0) = fpeak; 164 retval(1) = gpeak; 165 } 166 167 return retval; 168 } 169