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 Hankel singular values.
21 Uses SLICOT AB13AD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: January 2010
26 Version: 0.4
27 
28 */
29 
30 #include <octave/oct.h>
31 #include "common.h"
32 
33 extern "C"
34 {
35     int F77_FUNC (ab13ad, AB13AD)
36                  (char& DICO, char& EQUIL,
37                   F77_INT& N, F77_INT& M, F77_INT& P,
38                   double& ALPHA,
39                   double* A, F77_INT& LDA,
40                   double* B, F77_INT& LDB,
41                   double* C, F77_INT& LDC,
42                   F77_INT& NS,
43                   double* HSV,
44                   double* DWORK, F77_INT& LDWORK,
45                   F77_INT& INFO);
46 }
47 
48 // PKG_ADD: autoload ("__sl_ab13ad__", "__control_slicot_functions__.oct");
49 DEFUN_DLD (__sl_ab13ad__, args, nargout,
50    "-*- texinfo -*-\n\
51 Slicot AB13AD Release 5.0\n\
52 No argument checking.\n\
53 For internal use only.")
54 {
55     octave_idx_type nargin = args.length ();
56     octave_value_list retval;
57 
58     if (nargin != 6)
59     {
60         print_usage ();
61     }
62     else
63     {
64         // arguments in
65         char dico;
66         char equil;
67 
68         Matrix a = args(0).matrix_value ();
69         Matrix b = args(1).matrix_value ();
70         Matrix c = args(2).matrix_value ();
71         F77_INT discrete = args(3).int_value ();
72         double alpha = args(4).double_value ();
73         const F77_INT scaled = args(5).int_value ();
74 
75         if (discrete == 0)
76             dico = 'C';
77         else
78             dico = 'D';
79 
80         if (scaled == 0)
81             equil = 'S';
82         else
83             equil = 'N';
84 
85 
86         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
87         F77_INT m = TO_F77_INT (b.columns ());   // m: number of inputs
88         F77_INT p = TO_F77_INT (c.rows ());      // p: number of outputs
89 
90         F77_INT lda = max (1, TO_F77_INT (a.rows ()));
91         F77_INT ldb = max (1, TO_F77_INT (b.rows ()));
92         F77_INT ldc = max (1, TO_F77_INT (c.rows ()));
93 
94         // arguments out
95         F77_INT ns = 0;
96 
97         ColumnVector hsv (n);
98 
99         // workspace
100         F77_INT ldwork = max (1, n*(max (n, m, p) + 5) + n*(n+1)/2);
101 
102         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
103 
104         // error indicators
105         F77_INT info = 0;
106 
107 
108         // SLICOT routine AB13AD
109         F77_XFCN (ab13ad, AB13AD,
110                  (dico, equil,
111                   n, m, p,
112                   alpha,
113                   a.fortran_vec (), lda,
114                   b.fortran_vec (), ldb,
115                   c.fortran_vec (), ldc,
116                   ns,
117                   hsv.fortran_vec (),
118                   dwork, ldwork,
119                   info));
120 
121         if (f77_exception_encountered)
122             error ("hsvd: __sl_ab13ad__: exception in SLICOT subroutine AB13AD");
123 
124         if (info != 0)
125             error ("hsvd: __sl_ab13ad__: AB13AD returned info = %d", static_cast<int> (info));
126 
127         // resize
128         hsv.resize (ns);
129 
130         // return values
131         retval(0) = hsv;
132         retval(1) = octave_value (ns);
133     }
134 
135     return retval;
136 }
137