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 Square-root solver for generalized Lyapunov equations. 21 Uses SLICOT SG03BD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: September 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 (sg03bd, SG03BD) 36 (char& DICO, char& FACT, char& TRANS, 37 F77_INT& N, F77_INT& M, 38 double* A, F77_INT& LDA, 39 double* E, F77_INT& LDE, 40 double* Q, F77_INT& LDQ, 41 double* Z, F77_INT& LDZ, 42 double* B, F77_INT& LDB, 43 double& SCALE, 44 double* ALPHAR, double* ALPHAI, 45 double* BETA, 46 double* DWORK, F77_INT& LDWORK, 47 F77_INT& INFO); 48 } 49 50 // PKG_ADD: autoload ("__sl_sg03bd__", "__control_slicot_functions__.oct"); 51 DEFUN_DLD (__sl_sg03bd__, args, nargout, 52 "-*- texinfo -*-\n\ 53 Slicot SG03BD Release 5.0\n\ 54 No argument checking.\n\ 55 For internal use only.") 56 { 57 octave_idx_type nargin = args.length (); 58 octave_value_list retval; 59 60 if (nargin != 4) 61 { 62 print_usage (); 63 } 64 else 65 { 66 // arguments in 67 char dico; 68 char fact = 'N'; 69 char trans = 'N'; 70 71 Matrix a = args(0).matrix_value (); 72 Matrix e = args(1).matrix_value (); 73 Matrix b = args(2).matrix_value (); 74 F77_INT discrete = args(3).int_value (); 75 76 if (discrete == 0) 77 dico = 'C'; 78 else 79 dico = 'D'; 80 81 F77_INT n = TO_F77_INT (a.rows ()); 82 F77_INT m = TO_F77_INT (b.rows ()); 83 84 F77_INT lda = max (1, n); 85 F77_INT lde = max (1, n); 86 F77_INT ldq = max (1, n); 87 F77_INT ldz = max (1, n); 88 F77_INT ldb = max (1, m, n); 89 90 F77_INT n1 = max (m, n); 91 b.resize (ldb, n1); 92 93 // arguments out 94 double scale; 95 96 Matrix q (ldq, n); 97 Matrix z (ldz, n); 98 ColumnVector alphar (n); 99 ColumnVector alphai (n); 100 ColumnVector beta (n); 101 102 // workspace 103 F77_INT ldwork = 8*n + 16; 104 105 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 106 107 // error indicator 108 F77_INT info; 109 110 111 // SLICOT routine SG03BD 112 F77_XFCN (sg03bd, SG03BD, 113 (dico, fact, trans, 114 n, m, 115 a.fortran_vec (), lda, 116 e.fortran_vec (), lde, 117 q.fortran_vec (), ldq, 118 z.fortran_vec (), ldz, 119 b.fortran_vec (), ldb, 120 scale, 121 alphar.fortran_vec (), alphai.fortran_vec (), 122 beta.fortran_vec (), 123 dwork, ldwork, 124 info)); 125 126 if (f77_exception_encountered) 127 error ("lyap: __sl_sg03bd__: exception in SLICOT subroutine SG03BD"); 128 129 if (info != 0) 130 error ("lyap: __sl_sg03bd__: SG03BD returned info = %d", static_cast<int> (info)); 131 132 // resize 133 b.resize (n, n); 134 135 // return values 136 retval(0) = b; // b has been overwritten by cholesky factor u 137 retval(1) = octave_value (scale); 138 } 139 140 return retval; 141 } 142