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 Convert descriptor state-space system into regular state-space form.
21 Uses SLICOT SB10JD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: September 2011
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 (sb10jd, SB10JD)
36                  (F77_INT& N, F77_INT& M, F77_INT& NP,
37                   double* A, F77_INT& LDA,
38                   double* B, F77_INT& LDB,
39                   double* C, F77_INT& LDC,
40                   double* D, F77_INT& LDD,
41                   double* E, F77_INT& LDE,
42                   F77_INT& NSYS,
43                   double* DWORK, F77_INT& LDWORK,
44                   F77_INT& INFO);
45 }
46 
47 // PKG_ADD: autoload ("__sl_sb10jd__", "__control_slicot_functions__.oct");
48 DEFUN_DLD (__sl_sb10jd__, args, nargout,
49    "-*- texinfo -*-\n\
50 Slicot SB10JD Release 5.0\n\
51 No argument checking.\n\
52 For internal use only.")
53 {
54     octave_idx_type nargin = args.length ();
55     octave_value_list retval;
56 
57     if (nargin != 5)
58     {
59         print_usage ();
60     }
61     else
62     {
63         // arguments in
64         Matrix a = args(0).matrix_value ();
65         Matrix b = args(1).matrix_value ();
66         Matrix c = args(2).matrix_value ();
67         Matrix d = args(3).matrix_value ();
68         Matrix e = args(4).matrix_value ();
69 
70         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
71         F77_INT m = TO_F77_INT (b.columns ());   // m: number of inputs
72         F77_INT np = TO_F77_INT (c.rows ());     // np: number of outputs
73 
74         F77_INT lda = max (1, n);
75         F77_INT ldb = max (1, n);
76         F77_INT ldc = max (1, np);
77         F77_INT ldd = max (1, np);
78         F77_INT lde = max (1, n);
79 
80         // arguments out
81         F77_INT nsys;
82 
83         // workspace
84         F77_INT ldwork = max (1, 2*n*n + 2*n + n*max (5, n + m + np));
85         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
86 
87         // error indicator
88         F77_INT info;
89 
90 
91         // SLICOT routine SB10JD
92         F77_XFCN (sb10jd, SB10JD,
93                  (n, m, np,
94                   a.fortran_vec (), lda,
95                   b.fortran_vec (), ldb,
96                   c.fortran_vec (), ldc,
97                   d.fortran_vec (), ldd,
98                   e.fortran_vec (), lde,
99                   nsys,
100                   dwork, ldwork,
101                   info));
102 
103         if (f77_exception_encountered)
104             error ("__sl_sb10jd__: exception in SLICOT subroutine SB10JD");
105 
106         if (info != 0)
107             error ("__sl_sb10jd__: SB10JD returned info = %d", static_cast<int> (info));
108 
109         // resize
110         a.resize (nsys, nsys);
111         b.resize (nsys, m);
112         c.resize (np, nsys);
113 
114         // return values
115         retval(0) = a;
116         retval(1) = b;
117         retval(2) = c;
118         retval(3) = d;
119     }
120 
121     return retval;
122 }
123