1 /* 2 3 Copyright (C) 2014-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 Matrix exponential and integral for a real matrix. 21 Uses SLICOT MB05ND by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Thomas Vasileiou <thomas-v@wildmail.com> 25 Created: March 2014 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 (mb05nd, MB05ND) 36 (F77_INT& N, double& DELTA, 37 double* A, F77_INT& LDA, 38 double* EX, F77_INT& LDEX, 39 double* EXINT, F77_INT& LDEXINT, 40 double& TOL, 41 F77_INT* IWORK, 42 double* DWORK, F77_INT& LDWORK, 43 F77_INT& INFO); 44 } 45 46 // PKG_ADD: autoload ("__sl_mb05nd__", "__control_slicot_functions__.oct"); 47 DEFUN_DLD (__sl_mb05nd__, args, nargout, 48 "-*- texinfo -*-\n\ 49 Slicot MB05ND 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 double delta = args(1).double_value (); 65 double tol = args(2).double_value (); 66 67 if (a.numel () == 0) 68 { 69 // given matrix is empty, just return the empty matrix, 70 // otherise mb05d would make a [](1x0) matrix from it 71 retval(0) = a; 72 retval(1) = a; 73 return retval; 74 } 75 76 F77_INT n = TO_F77_INT (a.rows ()); 77 F77_INT lda = max (1, n); 78 F77_INT ldex = max (1, n); 79 F77_INT ldexin = max (1, n); 80 81 // arguments out 82 Matrix ex (ldex, n); 83 Matrix exint (ldexin, n); 84 85 // workspace 86 F77_INT ldwork = max (1, 2*n*n); // optimum performance 87 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, n); 88 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 89 90 // error indicators 91 F77_INT info = 0; 92 93 94 // SLICOT routine MB05ND 95 F77_XFCN (mb05nd, MB05ND, 96 (n, delta, 97 a.fortran_vec (), lda, 98 ex.fortran_vec (), ldex, 99 exint.fortran_vec (), ldexin, 100 tol, 101 iwork, 102 dwork, ldwork, 103 info)); 104 105 if (f77_exception_encountered) 106 error ("__sl_mb05nd__: exception in SLICOT subroutine MB05ND"); 107 108 if (info > 0) 109 { 110 if (info == n+1) 111 info = 2; 112 else 113 info = 1; 114 } 115 116 117 static const char* err_msg[] = { 118 "0: OK", 119 "1: some element of the denominator of the Pade " 120 "approximation is zero, so the denominator " 121 "is exactly singular.", 122 "2: DELTA = (delta * frobenius norm of matrix A) is " 123 "probably too large to permit meaningful computation. " 124 "That is, DELTA > SQRT(BIG), where BIG is a " 125 "representable number near the overflow threshold of " 126 "the machine."}; 127 128 error_msg ("__sl_mb05nd__", info, 2, err_msg); 129 130 // return values 131 retval(0) = ex; 132 retval(1) = exint; 133 } 134 135 return retval; 136 } 137