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 H-2 norm of a SS model. 21 Uses SLICOT AB13BD 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 "common.h" 32 33 extern "C" 34 { 35 double F77_FUNC (ab13bd, AB13BD) 36 (char& DICO, char& JOBN, 37 F77_INT& N, F77_INT& M, F77_INT& P, 38 double* A, F77_INT& LDA, 39 double* B, F77_INT& LDB, 40 double* C, F77_INT& LDC, 41 double* D, F77_INT& LDD, 42 F77_INT& NQ, 43 double& TOL, 44 double* DWORK, F77_INT& LDWORK, 45 F77_INT& IWARN, 46 F77_INT& INFO); 47 } 48 49 // PKG_ADD: autoload ("__sl_ab13bd__", "__control_slicot_functions__.oct"); 50 DEFUN_DLD (__sl_ab13bd__, args, nargout, 51 "-*- texinfo -*-\n\ 52 Slicot AB13BD Release 5.\n\ 53 No argument checking.\n\ 54 For internal use only.") 55 { 56 octave_idx_type nargin = args.length (); 57 octave_value_list retval; 58 59 if (nargin != 5) 60 { 61 print_usage (); 62 } 63 else 64 { 65 // arguments in 66 char dico; 67 char jobn = 'H'; 68 69 Matrix a = args(0).matrix_value (); 70 Matrix b = args(1).matrix_value (); 71 Matrix c = args(2).matrix_value (); 72 Matrix d = args(3).matrix_value (); 73 F77_INT discrete = args(4).int_value (); 74 75 if (discrete == 0) 76 dico = 'C'; 77 else 78 dico = 'D'; 79 80 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 81 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 82 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 83 84 F77_INT lda = max (1, TO_F77_INT (a.rows ())); 85 F77_INT ldb = max (1, TO_F77_INT (b.rows ())); 86 F77_INT ldc = max (1, TO_F77_INT (c.rows ())); 87 F77_INT ldd = max (1, TO_F77_INT (d.rows ())); 88 89 // arguments out 90 double norm; 91 F77_INT nq; 92 93 // tolerance 94 double tol = 0; 95 96 // workspace 97 F77_INT ldwork = max (1, m*(n+m) + max (n*(n+5), m*(m+2), 4*p ), 98 n*(max (n, p) + 4 ) + min (n, p)); 99 100 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 101 102 // error indicator 103 F77_INT iwarn; 104 F77_INT info; 105 106 107 // SLICOT routine AB13BD 108 norm = F77_FUNC (ab13bd, AB13BD) 109 (dico, jobn, 110 n, m, p, 111 a.fortran_vec (), lda, 112 b.fortran_vec (), ldb, 113 c.fortran_vec (), ldc, 114 d.fortran_vec (), ldd, 115 nq, 116 tol, 117 dwork, ldwork, 118 iwarn, 119 info); 120 121 if (f77_exception_encountered) 122 error ("lti: norm: __sl_ab13bd__: exception in SLICOT subroutine AB13BD"); 123 124 if (info != 0) 125 error ("lti: norm: __sl_ab13bd__: AB13BD returned info = %d", static_cast<int> (info)); 126 127 if (iwarn != 0) 128 warning ("lti: norm: __sl_ab13bd__: AB13BD returned iwarn = %d", static_cast<int> (iwarn)); 129 130 // return value 131 retval(0) = octave_value (norm); 132 } 133 134 return retval; 135 } 136