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 Minimal realization of state-space models.
21 Uses SLICOT TB01PD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: September 2010
26 Version: 0.5
27 
28 */
29 
30 #include <octave/oct.h>
31 #include "common.h"
32 
33 extern "C"
34 {
35     int F77_FUNC (tb01pd, TB01PD)
36                  (char& JOB, char& EQUIL,
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                   F77_INT& NR,
42                   double& TOL,
43                   F77_INT* IWORK,
44                   double* DWORK, F77_INT& LDWORK,
45                   F77_INT& INFO);
46 }
47 
48 // PKG_ADD: autoload ("__sl_tb01pd__", "__control_slicot_functions__.oct");
49 DEFUN_DLD (__sl_tb01pd__, args, nargout,
50    "-*- texinfo -*-\n\
51 Slicot TB01PD 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 != 5)
59     {
60         print_usage ();
61     }
62     else
63     {
64         // arguments in
65         char job = 'M';
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         double tol = args(3).double_value ();
72         const F77_INT scaled = args(4).int_value ();
73 
74         if (scaled == 0)
75             equil = 'S';
76         else
77             equil = 'N';
78 
79         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
80         F77_INT m = TO_F77_INT (b.columns ());   // m: number of inputs
81         F77_INT p = TO_F77_INT (c.rows ());      // p: number of outputs
82 
83         F77_INT lda = max (1, n);
84         F77_INT ldb = max (1, n);
85         F77_INT ldc;
86 
87         if (n == 0)
88             ldc = 1;
89         else
90             ldc = max (1, m, p);
91 
92         b.resize (ldb, max (m, p));
93         c.resize (ldc, n);
94 
95         // arguments out
96         F77_INT nr = 0;
97 
98         // workspace
99         F77_INT liwork = n + max (m, p);
100         F77_INT ldwork = max (1, n + max (n, 3*m, 3*p));
101 
102         OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork);
103         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
104 
105         // error indicators
106         F77_INT info = 0;
107 
108 
109         // SLICOT routine TB01PD
110         F77_XFCN (tb01pd, TB01PD,
111                  (job, equil,
112                   n, m, p,
113                   a.fortran_vec (), lda,
114                   b.fortran_vec (), ldb,
115                   c.fortran_vec (), ldc,
116                   nr,
117                   tol,
118                   iwork,
119                   dwork, ldwork,
120                   info));
121 
122         if (f77_exception_encountered)
123             error ("ss: minreal: __sl_tb01pd__: exception in SLICOT subroutine TB01PD");
124 
125         if (info != 0)
126             error ("ss: minreal: __sl_tb01pd__: TB01PD returned info = %d", static_cast<int> (info));
127 
128         // resize
129         a.resize (nr, nr);
130         b.resize (nr, m);
131         c.resize (p, nr);
132 
133         // return values
134         retval(0) = a;
135         retval(1) = b;
136         retval(2) = c;
137         retval(3) = octave_value (nr);
138     }
139 
140     return retval;
141 }
142