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 Staircase controllability form.
21 Uses SLICOT AB01OD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: August 2010
26 Version: 0.2
27 
28 */
29 
30 #include <octave/oct.h>
31 #include "common.h"
32 
33 extern "C"
34 {
35     int F77_FUNC (ab01od, AB01OD)
36                  (char& STAGES,
37                   char& JOBU, char& JOBV,
38                   F77_INT& N, F77_INT& M,
39                   double* A, F77_INT& LDA,
40                   double* B, F77_INT& LDB,
41                   double* U, F77_INT& LDU,
42                   double* V, F77_INT& LDV,
43                   F77_INT& NCONT, F77_INT& INDCON,
44                   F77_INT* KSTAIR,
45                   double& TOL,
46                   F77_INT* IWORK,
47                   double* DWORK, F77_INT& LDWORK,
48                   F77_INT& INFO);
49 }
50 
51 // PKG_ADD: autoload ("__sl_ab01od__", "__control_slicot_functions__.oct");
52 DEFUN_DLD (__sl_ab01od__, args, nargout,
53    "-*- texinfo -*-\n\
54 Slicot AB01OD Release 5.0\n\
55 No argument checking.\n\
56 For internal use only.")
57 {
58     octave_idx_type nargin = args.length ();
59     octave_value_list retval;
60 
61     if (nargin != 3)
62     {
63         print_usage ();
64     }
65     else
66     {
67         // arguments in
68         char stages = 'F';
69         char jobu = 'I';
70         char jobv = 'N';    // not referenced because stages = 'F'
71 
72         Matrix a = args(0).matrix_value ();
73         Matrix b = args(1).matrix_value ();
74         double tol = args(2).double_value ();
75 
76         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
77         F77_INT m = TO_F77_INT (b.columns ());   // m: number of inputs
78 
79         F77_INT lda = max (1, n);
80         F77_INT ldb = max (1, n);
81         F77_INT ldu = max (1, n);
82         F77_INT ldv = 1;
83 
84         // arguments out
85         Matrix u (ldu, n);
86         double* v = 0;          // not referenced because stages = 'F'
87 
88         F77_INT ncont;
89         F77_INT indcon;
90 
91         OCTAVE_LOCAL_BUFFER (F77_INT, kstair, n);
92 
93         // workspace
94         F77_INT ldwork = max (1, n + max (n, 3*m));
95 
96         OCTAVE_LOCAL_BUFFER (F77_INT, iwork, m);
97         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
98 
99         // error indicators
100         F77_INT info;
101 
102 
103         // SLICOT routine AB01OD
104         F77_XFCN (ab01od, AB01OD,
105                  (stages,
106                   jobu, jobv,
107                   n, m,
108                   a.fortran_vec (), lda,
109                   b.fortran_vec (), ldb,
110                   u.fortran_vec (), ldu,
111                   v, ldv,
112                   ncont, indcon,
113                   kstair,
114                   tol,
115                   iwork,
116                   dwork, ldwork,
117                   info));
118 
119         if (f77_exception_encountered)
120             error ("__sl_ab01od__: exception in SLICOT subroutine AB01OD");
121 
122         if (info != 0)
123             error ("__sl_ab01od__: AB01OD returned info = %d", static_cast<int> (info));
124 
125         // resize
126         a.resize (n, n);
127         b.resize (n, m);
128         u.resize (n, n);
129 
130         // return values
131         retval(0) = a;
132         retval(1) = b;
133         retval(2) = u;
134         retval(3) = octave_value (ncont);
135     }
136 
137     return retval;
138 }
139