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 Square-root solver for generalized Lyapunov equations.
21 Uses SLICOT SG03BD by courtesy of NICONET e.V.
22 <http://www.slicot.org>
23 
24 Author: Lukas Reichlin <lukas.reichlin@gmail.com>
25 Created: September 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 (sg03bd, SG03BD)
36                  (char& DICO, char& FACT, char& TRANS,
37                   F77_INT& N, F77_INT& M,
38                   double* A, F77_INT& LDA,
39                   double* E, F77_INT& LDE,
40                   double* Q, F77_INT& LDQ,
41                   double* Z, F77_INT& LDZ,
42                   double* B, F77_INT& LDB,
43                   double& SCALE,
44                   double* ALPHAR, double* ALPHAI,
45                   double* BETA,
46                   double* DWORK, F77_INT& LDWORK,
47                   F77_INT& INFO);
48 }
49 
50 // PKG_ADD: autoload ("__sl_sg03bd__", "__control_slicot_functions__.oct");
51 DEFUN_DLD (__sl_sg03bd__, args, nargout,
52    "-*- texinfo -*-\n\
53 Slicot SG03BD Release 5.0\n\
54 No argument checking.\n\
55 For internal use only.")
56 {
57     octave_idx_type nargin = args.length ();
58     octave_value_list retval;
59 
60     if (nargin != 4)
61     {
62         print_usage ();
63     }
64     else
65     {
66         // arguments in
67         char dico;
68         char fact = 'N';
69         char trans = 'N';
70 
71         Matrix a = args(0).matrix_value ();
72         Matrix e = args(1).matrix_value ();
73         Matrix b = args(2).matrix_value ();
74         F77_INT discrete = args(3).int_value ();
75 
76         if (discrete == 0)
77           dico = 'C';
78         else
79           dico = 'D';
80 
81         F77_INT n = TO_F77_INT (a.rows ());
82         F77_INT m = TO_F77_INT (b.rows ());
83 
84         F77_INT lda = max (1, n);
85         F77_INT lde = max (1, n);
86         F77_INT ldq = max (1, n);
87         F77_INT ldz = max (1, n);
88         F77_INT ldb = max (1, m, n);
89 
90         F77_INT n1 = max (m, n);
91         b.resize (ldb, n1);
92 
93         // arguments out
94         double scale;
95 
96         Matrix q (ldq, n);
97         Matrix z (ldz, n);
98         ColumnVector alphar (n);
99         ColumnVector alphai (n);
100         ColumnVector beta (n);
101 
102         // workspace
103         F77_INT ldwork = 8*n + 16;
104 
105         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
106 
107         // error indicator
108         F77_INT info;
109 
110 
111         // SLICOT routine SG03BD
112         F77_XFCN (sg03bd, SG03BD,
113                  (dico, fact, trans,
114                   n, m,
115                   a.fortran_vec (), lda,
116                   e.fortran_vec (), lde,
117                   q.fortran_vec (), ldq,
118                   z.fortran_vec (), ldz,
119                   b.fortran_vec (), ldb,
120                   scale,
121                   alphar.fortran_vec (), alphai.fortran_vec (),
122                   beta.fortran_vec (),
123                   dwork, ldwork,
124                   info));
125 
126         if (f77_exception_encountered)
127             error ("lyap: __sl_sg03bd__: exception in SLICOT subroutine SG03BD");
128 
129         if (info != 0)
130             error ("lyap: __sl_sg03bd__: SG03BD returned info = %d", static_cast<int> (info));
131 
132         // resize
133         b.resize (n, n);
134 
135         // return values
136         retval(0) = b;  // b has been overwritten by cholesky factor u
137         retval(1) = octave_value (scale);
138     }
139 
140     return retval;
141 }
142