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