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