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 Model reduction based on balanced stochastic truncation method. 21 Uses SLICOT AB09HD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: October 2011 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 (ab09hd, AB09HD) 36 (char& DICO, char& JOB, char& EQUIL, char& ORDSEL, 37 F77_INT& N, F77_INT& M, F77_INT& P, 38 F77_INT& NR, 39 double& ALPHA, double& BETA, 40 double* A, F77_INT& LDA, 41 double* B, F77_INT& LDB, 42 double* C, F77_INT& LDC, 43 double* D, F77_INT& LDD, 44 F77_INT& NS, 45 double* HSV, 46 double& TOL1, double& TOL2, 47 F77_INT* IWORK, 48 double* DWORK, F77_INT& LDWORK, 49 F77_LOGICAL* BWORK, 50 F77_INT& IWARN, F77_INT& INFO); 51 } 52 53 // PKG_ADD: autoload ("__sl_ab09hd__", "__control_slicot_functions__.oct"); 54 DEFUN_DLD (__sl_ab09hd__, args, nargout, 55 "-*- texinfo -*-\n\ 56 Slicot AB09HD Release 5.0\n\ 57 No argument checking.\n\ 58 For internal use only.") 59 { 60 octave_idx_type nargin = args.length (); 61 octave_value_list retval; 62 63 if (nargin != 13) 64 { 65 print_usage (); 66 } 67 else 68 { 69 // arguments in 70 char dico; 71 char job; 72 char equil; 73 char ordsel; 74 75 Matrix a = args(0).matrix_value (); 76 Matrix b = args(1).matrix_value (); 77 Matrix c = args(2).matrix_value (); 78 Matrix d = args(3).matrix_value (); 79 80 const F77_INT idico = args(4).int_value (); 81 const F77_INT iequil = args(5).int_value (); 82 const F77_INT ijob = args(6).int_value (); 83 84 F77_INT nr = args(7).int_value (); 85 const F77_INT iordsel = args(8).int_value (); 86 87 double alpha = args(9).double_value (); 88 double beta = args(10).double_value (); 89 double tol1 = args(11).double_value (); 90 double tol2 = args(12).double_value (); 91 92 switch (ijob) 93 { 94 case 0: 95 job = 'B'; 96 break; 97 case 1: 98 job = 'F'; 99 break; 100 case 2: 101 job = 'S'; 102 break; 103 case 3: 104 job = 'P'; 105 break; 106 default: 107 error ("__sl_ab09hd__: argument job invalid"); 108 } 109 110 if (idico == 0) 111 dico = 'C'; 112 else 113 dico = 'D'; 114 115 if (iequil == 0) 116 equil = 'S'; 117 else 118 equil = 'N'; 119 120 if (iordsel == 0) 121 ordsel = 'F'; 122 else 123 ordsel = 'A'; 124 125 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 126 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 127 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 128 129 F77_INT lda = max (1, n); 130 F77_INT ldb = max (1, n); 131 F77_INT ldc = max (1, p); 132 F77_INT ldd = max (1, p); 133 134 // arguments out 135 F77_INT ns; 136 ColumnVector hsv (n); 137 138 // workspace 139 F77_INT liwork = max (1, 2*n); 140 F77_INT mb; 141 142 if (beta == 0) 143 mb = m; 144 else 145 mb = m + p; 146 147 F77_INT ldwork = 2*n*n + mb*(n+p) 148 + max (2, n*(max (n,mb,p)+5), 2*n*p + max (p*(mb+2), 10*n*(n+1))); 149 150 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 151 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 152 OCTAVE_LOCAL_BUFFER (F77_LOGICAL, bwork, 2*n); 153 154 // error indicators 155 F77_INT iwarn = 0; 156 F77_INT info = 0; 157 158 159 // SLICOT routine AB09HD 160 F77_XFCN (ab09hd, AB09HD, 161 (dico, job, equil, ordsel, 162 n, m, p, 163 nr, 164 alpha, beta, 165 a.fortran_vec (), lda, 166 b.fortran_vec (), ldb, 167 c.fortran_vec (), ldc, 168 d.fortran_vec (), ldd, 169 ns, 170 hsv.fortran_vec (), 171 tol1, tol2, 172 iwork, 173 dwork, ldwork, 174 bwork, 175 iwarn, info)); 176 177 if (f77_exception_encountered) 178 error ("bstmodred: exception in SLICOT subroutine AB09HD"); 179 180 181 static const char* err_msg[] = { 182 "0: OK", 183 "1: the computation of the ordered real Schur form of A " 184 "failed", 185 "2: the reduction of the Hamiltonian matrix to real " 186 "Schur form failed", 187 "3: the reordering of the real Schur form of the " 188 "Hamiltonian matrix failed", 189 "4: the Hamiltonian matrix has less than N stable " 190 "eigenvalues", 191 "5: the coefficient matrix U11 in the linear system " 192 "X*U11 = U21 to determine X is singular to working " 193 "precision", 194 "6: BETA = 0 and D has not a maximal row rank", 195 "7: the computation of Hankel singular values failed", 196 "8: the separation of the ALPHA-stable/unstable diagonal " 197 "blocks failed because of very close eigenvalues", 198 "9: the resulting order of reduced stable part is less " 199 "than the number of unstable zeros of the stable " 200 "part"}; 201 202 static const char* warn_msg[] = { 203 "0: OK", 204 "1: with ORDSEL = 'F', the selected order NR is greater " 205 "than NSMIN, the sum of the order of the " 206 "ALPHA-unstable part and the order of a minimal " 207 "realization of the ALPHA-stable part of the given " 208 "system; in this case, the resulting NR is set equal " 209 "to NSMIN.", 210 "2: with ORDSEL = 'F', the selected order NR corresponds " 211 "to repeated singular values for the ALPHA-stable " 212 "part, which are neither all included nor all " 213 "excluded from the reduced model; in this case, the " 214 "resulting NR is automatically decreased to exclude " 215 "all repeated singular values.", 216 "3: with ORDSEL = 'F', the selected order NR is less " 217 "than the order of the ALPHA-unstable part of the " 218 "given system; in this case NR is set equal to the " 219 "order of the ALPHA-unstable part."}; 220 221 error_msg ("bstmodred", info, 9, err_msg); 222 warning_msg ("bstmodred", iwarn, 3, warn_msg); 223 224 // resize 225 a.resize (nr, nr); 226 b.resize (nr, m); 227 c.resize (p, nr); 228 hsv.resize (ns); 229 230 // return values 231 retval(0) = a; 232 retval(1) = b; 233 retval(2) = c; 234 retval(3) = d; 235 retval(4) = octave_value (nr); 236 retval(5) = hsv; 237 retval(6) = octave_value (ns); 238 } 239 240 return retval; 241 } 242