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 Solution of continuous-time Sylvester equations.
21 Uses SLICOT SB04MD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: January 2010
26 Version: 0.3
27 
28 */
29 
30 #include <octave/oct.h>
31 #include "common.h"
32 
33 extern "C"
34 {
35     int F77_FUNC (sb04md, SB04MD)
36                  (F77_INT& N, F77_INT& M,
37                   double* A, F77_INT& LDA,
38                   double* B, F77_INT& LDB,
39                   double* C, F77_INT& LDC,
40                   double* Z, F77_INT& LDZ,
41                   F77_INT* IWORK,
42                   double* DWORK, F77_INT& LDWORK,
43                   F77_INT& INFO);
44 }
45 
46 // PKG_ADD: autoload ("__sl_sb04md__", "__control_slicot_functions__.oct");
47 DEFUN_DLD (__sl_sb04md__, args, nargout,
48    "-*- texinfo -*-\n\
49 Slicot SB04MD 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 != 3)
57     {
58         print_usage ();
59     }
60     else
61     {
62         // arguments in
63         Matrix a = args(0).matrix_value ();
64         Matrix b = args(1).matrix_value ();
65         Matrix c = args(2).matrix_value ();
66 
67         F77_INT n = TO_F77_INT (a.rows ());
68         F77_INT m = TO_F77_INT (b.rows ());
69 
70         F77_INT lda = max (1, n);
71         F77_INT ldb = max (1, m);
72         F77_INT ldc = max (1, n);
73         F77_INT ldz = max (1, m);
74 
75         // arguments out
76         Matrix z (ldz, m);
77 
78         // workspace
79         F77_INT ldwork = max (1, 2*n*n + 8*n, 5*m, n + m);
80 
81         OCTAVE_LOCAL_BUFFER (F77_INT, iwork, 4*n);
82         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
83 
84         // error indicator
85         F77_INT info;
86 
87 
88         // SLICOT routine SB04MD
89         F77_XFCN (sb04md, SB04MD,
90                  (n, m,
91                   a.fortran_vec (), lda,
92                   b.fortran_vec (), ldb,
93                   c.fortran_vec (), ldc,
94                   z.fortran_vec (), ldz,
95                   iwork,
96                   dwork, ldwork,
97                   info));
98 
99         if (f77_exception_encountered)
100             error ("lyap: __sl_sb04md__: exception in SLICOT subroutine SB04MD");
101 
102         if (info != 0)
103             error ("lyap: __sl_sb04md__: SB04MD returned info = %d", static_cast<int> (info));
104 
105         // return values
106         retval(0) = c;
107     }
108 
109     return retval;
110 }
111