1 /* 2 3 Copyright (C) 2013-2015 Thomas Vasileiou 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 Orthogonal reduction of a descriptor system to a SVD-like coordinate form. 21 Uses SLICOT TG01FD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Thomas Vasileiou <thomas-v@wildmail.com> 25 Created: September 2013 26 Version: 0.2 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 33 extern "C" 34 { 35 int F77_FUNC (tg01fd, TG01FD) 36 (char& COMPQ, char& COMPZ, char& JOBA, 37 F77_INT& L, F77_INT& N, F77_INT& M, F77_INT& P, 38 double* A, F77_INT& LDA, 39 double* E, F77_INT& LDE, 40 double* B, F77_INT& LDB, 41 double* C, F77_INT& LDC, 42 double* Q, F77_INT& LDQ, 43 double* Z, F77_INT& LDZ, 44 F77_INT& RANKE, F77_INT& RNKA22, 45 double& TOL, 46 F77_INT* IWORK, 47 double* DWORK, F77_INT& LDWORK, 48 F77_INT& INFO); 49 } 50 51 // PKG_ADD: autoload ("__sl_tg01fd__", "__control_slicot_functions__.oct"); 52 DEFUN_DLD (__sl_tg01fd__, args, nargout, 53 "-*- texinfo -*-\n\ 54 Slicot TG01FD Release 5.0\n\ 55 No argument checking.\n\ 56 For internal use only.") 57 { 58 octave_idx_type nargin = args.length (); 59 octave_value_list retval; 60 61 if (nargin != 6) 62 { 63 print_usage (); 64 } 65 else 66 { 67 // arguments in 68 char compq; 69 char compz; 70 char joba = 'T'; 71 72 Matrix a = args(0).matrix_value (); 73 Matrix e = args(1).matrix_value (); 74 Matrix b = args(2).matrix_value (); 75 Matrix c = args(3).matrix_value (); 76 const F77_INT qz_flag = args(4).int_value (); 77 double tol = args(5).double_value (); 78 79 if (qz_flag == 0) 80 { 81 compq = 'N'; 82 compz = 'N'; 83 } 84 else 85 { 86 compq = 'I'; 87 compz = 'I'; 88 } 89 90 F77_INT l = TO_F77_INT (a.rows ()); 91 F77_INT n = l; 92 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 93 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 94 95 F77_INT lda = max (1, l); 96 F77_INT lde = max (1, l); 97 F77_INT ldb = max (1, l); 98 F77_INT ldc = max (1, p); 99 F77_INT ldq = max (1, l); 100 F77_INT ldz = max (1, n); 101 102 // arguments out 103 Matrix q(l, l, 0.); 104 Matrix z(n, n, 0.); 105 Matrix empty(0, 0); 106 F77_INT ranke, rnka22; 107 108 // workspace 109 F77_INT ldwork = max (1, n+p, min (l,n) + max (3*n-1, m, l)); 110 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, n); 111 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 112 113 // error indicators 114 F77_INT info = 0; 115 116 117 // SLICOT routine TG01FD 118 F77_XFCN (tg01fd, TG01FD, 119 (compq, compz, joba, 120 l, n, m, p, 121 a.fortran_vec (), lda, 122 e.fortran_vec (), lde, 123 b.fortran_vec (), ldb, 124 c.fortran_vec (), ldc, 125 q.fortran_vec (), ldq, 126 z.fortran_vec (), ldz, 127 ranke, rnka22, 128 tol, 129 iwork, 130 dwork, ldwork, 131 info)); 132 133 if (f77_exception_encountered) 134 error ("__sl_tg01fd__: exception in SLICOT subroutine TG01FD"); 135 136 if (info != 0) 137 error ("__sl_tg01fd__: TG01FD returned info = %d", static_cast<int> (info)); 138 139 // return values 140 retval(0) = a; 141 retval(1) = e; 142 retval(2) = b; 143 retval(3) = c; 144 retval(4) = octave_value (ranke); 145 retval(5) = octave_value (rnka22); 146 if (qz_flag == 0) 147 { 148 retval(6) = empty; 149 retval(7) = empty; 150 } 151 else 152 { 153 retval(6) = q; 154 retval(7) = z; 155 } 156 } 157 158 return retval; 159 } 160