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 Balance state-space model.
21 Uses SLICOT TB01ID by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: May 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 (tb01id, TB01ID)
36                  (char& JOB,
37                   F77_INT& N, F77_INT& M, F77_INT& P,
38                   double& MAXRED,
39                   double* A, F77_INT& LDA,
40                   double* B, F77_INT& LDB,
41                   double* C, F77_INT& LDC,
42                   double* SCALE,
43                   F77_INT& INFO);
44 }
45 
46 // PKG_ADD: autoload ("__sl_tb01id__", "__control_slicot_functions__.oct");
47 DEFUN_DLD (__sl_tb01id__, args, nargout,
48    "-*- texinfo -*-\n\
49 Slicot TB01ID Release 5.0\n\
50 No argument checking.\n\
51 For internal use only.")
52 {
53     octave_idx_type nargin = args.length ();
54     octave_value_list retval;
55 
56     if (nargin != 4)
57     {
58         print_usage ();
59     }
60     else
61     {
62         // arguments in
63         char job = 'A';
64 
65         Matrix a = args(0).matrix_value ();
66         Matrix b = args(1).matrix_value ();
67         Matrix c = args(2).matrix_value ();
68         double maxred = args(3).double_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 p = TO_F77_INT (c.rows ());      // p: number of outputs
73 
74         F77_INT lda = max (1, n);
75         F77_INT ldb = max (1, n);
76         F77_INT ldc = max (1, p);
77 
78 
79         // arguments out
80         ColumnVector scale (n);
81 
82 
83         // error indicators
84         F77_INT info = 0;
85 
86 
87         // SLICOT routine TB01ID
88         F77_XFCN (tb01id, TB01ID,
89                  (job,
90                   n, m, p,
91                   maxred,
92                   a.fortran_vec (), lda,
93                   b.fortran_vec (), ldb,
94                   c.fortran_vec (), ldc,
95                   scale.fortran_vec (),
96                   info));
97 
98         if (f77_exception_encountered)
99             error ("ss: prescale: __sl_tb01id__: exception in SLICOT subroutine TB01ID");
100 
101         if (info != 0)
102             error ("ss: prescale: __sl_tb01id__: TB01ID returned info = %d", static_cast<int> (info));
103 
104 
105         // return values
106         retval(0) = a;
107         retval(1) = b;
108         retval(2) = c;
109         retval(3) = octave_value (maxred);
110         retval(4) = scale;
111     }
112 
113     return retval;
114 }
115