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 generalized Lyapunov equations.
21 Uses SLICOT SG03AD 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 (sg03ad, SG03AD)
36                  (char& DICO, char& JOB,
37                   char& FACT, char& TRANS,
38                   char& UPLO,
39                   F77_INT& N,
40                   double* A, F77_INT& LDA,
41                   double* E, F77_INT& LDE,
42                   double* Q, F77_INT& LDQ,
43                   double* Z, F77_INT& LDZ,
44                   double* X, F77_INT& LDX,
45                   double& SCALE,
46                   double& SEP, double& FERR,
47                   double* ALPHAR, double* ALPHAI,
48                   double* BETA,
49                   F77_INT* IWORK,
50                   double* DWORK, F77_INT& LDWORK,
51                   F77_INT& INFO);
52 }
53 
54 // PKG_ADD: autoload ("__sl_sg03ad__", "__control_slicot_functions__.oct");
55 DEFUN_DLD (__sl_sg03ad__, args, nargout,
56    "-*- texinfo -*-\n\
57 Slicot SG03AD Release 5.0\n\
58 No argument checking.\n\
59 For internal use only.")
60 {
61     octave_idx_type nargin = args.length ();
62     octave_value_list retval;
63 
64     if (nargin != 4)
65     {
66         print_usage ();
67     }
68     else
69     {
70         // arguments in
71         char dico;
72         char job = 'X';
73         char fact = 'N';
74         char trans = 'T';
75         char uplo = 'U';    // ?!?
76 
77         Matrix a = args(0).matrix_value ();
78         Matrix e = args(1).matrix_value ();
79         Matrix x = args(2).matrix_value ();
80         F77_INT discrete = args(3).int_value ();
81 
82         if (discrete == 0)
83           dico = 'C';
84         else
85           dico = 'D';
86 
87         F77_INT n = TO_F77_INT (a.rows ());      // n: number of states
88 
89         F77_INT lda = max (1, n);
90         F77_INT lde = max (1, n);
91         F77_INT ldq = max (1, n);
92         F77_INT ldz = max (1, n);
93         F77_INT ldx = max (1, n);
94 
95         // arguments out
96         double scale;
97         double sep = 0;
98         double ferr = 0;
99 
100         Matrix q (ldq, n);
101         Matrix z (ldz, n);
102         ColumnVector alphar (n);
103         ColumnVector alphai (n);
104         ColumnVector beta (n);
105 
106         // workspace
107         F77_INT* iwork = 0;  // not referenced because job = X
108 
109         F77_INT ldwork = max (1, 8*n+16);
110         OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
111 
112         // error indicator
113         F77_INT info;
114 
115 
116         // SLICOT routine SG03AD
117         F77_XFCN (sg03ad, SG03AD,
118                  (dico, job,
119                   fact, trans,
120                   uplo,
121                   n,
122                   a.fortran_vec (), lda,
123                   e.fortran_vec (), lde,
124                   q.fortran_vec (), ldq,
125                   z.fortran_vec (), ldz,
126                   x.fortran_vec (), ldx,
127                   scale,
128                   sep, ferr,
129                   alphar.fortran_vec (), alphai.fortran_vec (),
130                   beta.fortran_vec (),
131                   iwork,
132                   dwork, ldwork,
133                   info));
134 
135         if (f77_exception_encountered)
136             error ("lyap: __sl_sg03ad__: exception in SLICOT subroutine SG03AD");
137 
138         if (info != 0)
139             error ("lyap: __sl_sg03ad__: SG03AD returned info = %d", static_cast<int> (info));
140 
141         // return values
142         retval(0) = x;
143         retval(1) = octave_value (scale);
144     }
145 
146     return retval;
147 }
148