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 discrete-time Sylvester equations. 21 Uses SLICOT SB04QD 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 (sb04qd, SB04QD) 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_sb04qd__", "__control_slicot_functions__.oct"); 47 DEFUN_DLD (__sl_sb04qd__, args, nargout, 48 "-*- texinfo -*-\n\ 49 Slicot SB04QD 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 + 9*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 SB04QD 89 F77_XFCN (sb04qd, SB04QD, 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 ("dlyap: __sl_sb04qd__: exception in SLICOT subroutine SB04QD"); 101 102 if (info != 0) 103 error ("dlyap: __sl_sb04qd__: SB04QD returned info = %d", static_cast<int> (info)); 104 105 // return values 106 retval(0) = c; 107 } 108 109 return retval; 110 } 111