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