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 Gain of descriptor state-space models.  Based on SLICOT TB04BX.f.
21 
22 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
23 Created: March 2011
24 Version: 0.2
25 
26 */
27 
28 #include <octave/oct.h>
29 #include "common.h"
30 #include <complex>
31 #include <xpow.h>
32 
33 extern "C"
34 {
35     int F77_FUNC (tg04bx, TG04BX)
36                  (F77_INT& IP, F77_INT& IZ,
37                   double* A, F77_INT& LDA,
38                   double* E,
39                   double* B,
40                   double* C,
41                   double* D,
42                   double* PR, double* PI,
43                   double* ZR, double* ZI,
44                   double& GAIN,
45                   F77_INT* IWORK);
46 }
47 
48 // PKG_ADD: autoload ("__sl_tg04bx__", "__control_slicot_functions__.oct");
49 DEFUN_DLD (__sl_tg04bx__, args, nargout,
50    "-*- texinfo -*-\n\
51 Slicot TG04BX 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 != 9)
59     {
60         print_usage ();
61     }
62     else
63     {
64         // arguments in
65         Matrix a = args(0).matrix_value ();
66         Matrix e = args(1).matrix_value ();
67         Matrix b = args(2).matrix_value ();
68         Matrix c = args(3).matrix_value ();
69         Matrix d = args(4).matrix_value ();
70 
71         ColumnVector pr = args(5).column_vector_value ();
72         ColumnVector pi = args(6).column_vector_value ();
73 
74         ColumnVector zr = args(7).column_vector_value ();
75         ColumnVector zi = args(8).column_vector_value ();
76 
77         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
78         F77_INT ip = TO_F77_INT (pr.numel ());  // ip: number of finite poles
79         F77_INT iz = TO_F77_INT (zr.numel ());  // iz: number of zeros
80 
81         // For ss, IP = n is always true.
82         // However, dss models with poles at infinity
83         // (filtered by pole.m) may have IP <= n
84 
85         // Take pr.length == pi.length == ip for granted,
86         // and the same for iz, zr and zi.
87 
88         F77_INT lda = max (1, n);
89 
90         // arguments out
91         double gain;
92 
93         // workspace
94         OCTAVE_LOCAL_BUFFER (F77_INT, iwork, lda);
95 
96 
97         F77_XFCN (tg04bx, TG04BX,
98                  (ip, iz,
99                   a.fortran_vec (), lda,
100                   e.fortran_vec (),
101                   b.fortran_vec (),
102                   c.fortran_vec (),
103                   d.fortran_vec (),
104                   pr.fortran_vec (), pi.fortran_vec (),
105                   zr.fortran_vec (), zi.fortran_vec (),
106                   gain,
107                   iwork));
108 
109         if (f77_exception_encountered)
110             error ("dss: zero: __sl_tg04bx__: exception in TG04BX");
111 
112         // return values
113         retval(0) = octave_value (gain);
114     }
115 
116     return retval;
117 }
118