GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
gepbalance.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2018 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software: you can redistribute it and/or modify it
8 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 Octave is distributed in the hope that it will be useful, but
13 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 Octave; see the file COPYING. If not, see
19 <https://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "CMatrix.h"
28 #include "dMatrix.h"
29 #include "fCMatrix.h"
30 #include "fMatrix.h"
31 #include "gepbalance.h"
32 #include "lo-array-errwarn.h"
33 #include "lo-error.h"
34 #include "lo-lapack-proto.h"
35 #include "oct-locbuf.h"
36 #include "quit.h"
37 
38 namespace octave
39 {
40  namespace math
41  {
42  template <>
43  octave_idx_type
44  gepbalance<Matrix>::init (const Matrix& a, const Matrix& b,
45  const std::string& balance_job)
46  {
47  F77_INT n = to_f77_int (a.cols ());
48 
49  if (a.rows () != n)
50  (*current_liboctave_error_handler)
51  ("GEPBALANCE requires square matrix");
52 
53  if (a.dims () != b.dims ())
54  err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
55 
56  F77_INT info;
57  F77_INT ilo;
58  F77_INT ihi;
59 
60  OCTAVE_LOCAL_BUFFER (double, plscale, n);
61  OCTAVE_LOCAL_BUFFER (double, prscale, n);
62  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
63 
64  balanced_mat = a;
65  double *p_balanced_mat = balanced_mat.fortran_vec ();
66  balanced_mat2 = b;
67  double *p_balanced_mat2 = balanced_mat2.fortran_vec ();
68 
69  char job = balance_job[0];
70 
71  F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
72  n, p_balanced_mat, n, p_balanced_mat2,
73  n, ilo, ihi, plscale, prscale, pwork, info
74  F77_CHAR_ARG_LEN (1)));
75 
76  balancing_mat = Matrix (n, n, 0.0);
77  balancing_mat2 = Matrix (n, n, 0.0);
78  for (F77_INT i = 0; i < n; i++)
79  {
80  octave_quit ();
81  balancing_mat.elem (i ,i) = 1.0;
82  balancing_mat2.elem (i ,i) = 1.0;
83  }
84 
85  double *p_balancing_mat = balancing_mat.fortran_vec ();
86  double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
87 
88  // first left
89  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
90  F77_CONST_CHAR_ARG2 ("L", 1),
91  n, ilo, ihi, plscale, prscale,
92  n, p_balancing_mat, n, info
93  F77_CHAR_ARG_LEN (1)
94  F77_CHAR_ARG_LEN (1)));
95 
96  // then right
97  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
98  F77_CONST_CHAR_ARG2 ("R", 1),
99  n, ilo, ihi, plscale, prscale,
100  n, p_balancing_mat2, n, info
101  F77_CHAR_ARG_LEN (1)
102  F77_CHAR_ARG_LEN (1)));
103 
104  return info;
105  }
106 
107  template <>
108  octave_idx_type
109  gepbalance<FloatMatrix>::init (const FloatMatrix& a, const FloatMatrix& b,
110  const std::string& balance_job)
111  {
112  F77_INT n = to_f77_int (a.cols ());
113 
114  if (a.rows () != n)
115  (*current_liboctave_error_handler)
116  ("FloatGEPBALANCE requires square matrix");
117 
118  if (a.dims () != b.dims ())
119  err_nonconformant ("FloatGEPBALANCE",
120  n, n, b.rows(), b.cols());
121 
122  F77_INT info;
123  F77_INT ilo;
124  F77_INT ihi;
125 
126  OCTAVE_LOCAL_BUFFER (float, plscale, n);
127  OCTAVE_LOCAL_BUFFER (float, prscale, n);
128  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
129 
130  balanced_mat = a;
131  float *p_balanced_mat = balanced_mat.fortran_vec ();
132  balanced_mat2 = b;
133  float *p_balanced_mat2 = balanced_mat2.fortran_vec ();
134 
135  char job = balance_job[0];
136 
137  F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
138  n, p_balanced_mat, n, p_balanced_mat2,
139  n, ilo, ihi, plscale, prscale, pwork, info
140  F77_CHAR_ARG_LEN (1)));
141 
142  balancing_mat = FloatMatrix (n, n, 0.0);
143  balancing_mat2 = FloatMatrix (n, n, 0.0);
144  for (F77_INT i = 0; i < n; i++)
145  {
146  octave_quit ();
147  balancing_mat.elem (i ,i) = 1.0;
148  balancing_mat2.elem (i ,i) = 1.0;
149  }
150 
151  float *p_balancing_mat = balancing_mat.fortran_vec ();
152  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
153 
154  // first left
155  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
156  F77_CONST_CHAR_ARG2 ("L", 1),
157  n, ilo, ihi, plscale, prscale,
158  n, p_balancing_mat, n, info
159  F77_CHAR_ARG_LEN (1)
160  F77_CHAR_ARG_LEN (1)));
161 
162  // then right
163  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
164  F77_CONST_CHAR_ARG2 ("R", 1),
165  n, ilo, ihi, plscale, prscale,
166  n, p_balancing_mat2, n, info
167  F77_CHAR_ARG_LEN (1)
168  F77_CHAR_ARG_LEN (1)));
169 
170  return info;
171  }
172 
173  template <>
174  octave_idx_type
175  gepbalance<ComplexMatrix>::init (const ComplexMatrix& a,
176  const ComplexMatrix& b,
177  const std::string& balance_job)
178  {
179  F77_INT n = to_f77_int (a.cols ());
180 
181  if (a.rows () != n)
182  (*current_liboctave_error_handler)
183  ("ComplexGEPBALANCE requires square matrix");
184 
185  if (a.dims () != b.dims ())
186  err_nonconformant ("ComplexGEPBALANCE",
187  n, n, b.rows(), b.cols());
188 
189  F77_INT info;
190  F77_INT ilo;
191  F77_INT ihi;
192 
193  OCTAVE_LOCAL_BUFFER (double, plscale, n);
194  OCTAVE_LOCAL_BUFFER (double, prscale, n);
195  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
196 
197  balanced_mat = a;
198  Complex *p_balanced_mat = balanced_mat.fortran_vec ();
199  balanced_mat2 = b;
200  Complex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
201 
202  char job = balance_job[0];
203 
204  F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
205  n, F77_DBLE_CMPLX_ARG (p_balanced_mat),
206  n, F77_DBLE_CMPLX_ARG (p_balanced_mat2),
207  n, ilo, ihi, plscale, prscale, pwork, info
208  F77_CHAR_ARG_LEN (1)));
209 
210  balancing_mat = Matrix (n, n, 0.0);
211  balancing_mat2 = Matrix (n, n, 0.0);
212  for (F77_INT i = 0; i < n; i++)
213  {
214  octave_quit ();
215  balancing_mat.elem (i ,i) = 1.0;
216  balancing_mat2.elem (i ,i) = 1.0;
217  }
218 
219  double *p_balancing_mat = balancing_mat.fortran_vec ();
220  double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
221 
222  // first left
223  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
224  F77_CONST_CHAR_ARG2 ("L", 1),
225  n, ilo, ihi, plscale, prscale,
226  n, p_balancing_mat, n, info
227  F77_CHAR_ARG_LEN (1)
228  F77_CHAR_ARG_LEN (1)));
229 
230  // then right
231  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
232  F77_CONST_CHAR_ARG2 ("R", 1),
233  n, ilo, ihi, plscale, prscale,
234  n, p_balancing_mat2, n, info
235  F77_CHAR_ARG_LEN (1)
236  F77_CHAR_ARG_LEN (1)));
237 
238  return info;
239  }
240 
241  template <>
242  octave_idx_type
243  gepbalance<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
244  const FloatComplexMatrix& b,
245  const std::string& balance_job)
246  {
247  F77_INT n = to_f77_int (a.cols ());
248 
249  if (a.rows () != n)
250  {
251  (*current_liboctave_error_handler)
252  ("FloatComplexGEPBALANCE requires square matrix");
253  return -1;
254  }
255 
256  if (a.dims () != b.dims ())
257  err_nonconformant ("FloatComplexGEPBALANCE",
258  n, n, b.rows(), b.cols());
259 
260  F77_INT info;
261  F77_INT ilo;
262  F77_INT ihi;
263 
264  OCTAVE_LOCAL_BUFFER (float, plscale, n);
265  OCTAVE_LOCAL_BUFFER (float, prscale, n);
266  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
267 
268  balanced_mat = a;
269  FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
270  balanced_mat2 = b;
271  FloatComplex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
272 
273  char job = balance_job[0];
274 
275  F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
276  n, F77_CMPLX_ARG (p_balanced_mat),
277  n, F77_CMPLX_ARG (p_balanced_mat2),
278  n, ilo, ihi, plscale, prscale, pwork, info
279  F77_CHAR_ARG_LEN (1)));
280 
281  balancing_mat = FloatMatrix (n, n, 0.0);
282  balancing_mat2 = FloatMatrix (n, n, 0.0);
283  for (F77_INT i = 0; i < n; i++)
284  {
285  octave_quit ();
286  balancing_mat.elem (i ,i) = 1.0;
287  balancing_mat2.elem (i ,i) = 1.0;
288  }
289 
290  float *p_balancing_mat = balancing_mat.fortran_vec ();
291  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
292 
293  // first left
294  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
295  F77_CONST_CHAR_ARG2 ("L", 1),
296  n, ilo, ihi, plscale, prscale,
297  n, p_balancing_mat, n, info
298  F77_CHAR_ARG_LEN (1)
299  F77_CHAR_ARG_LEN (1)));
300 
301  // then right
302  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
303  F77_CONST_CHAR_ARG2 ("R", 1),
304  n, ilo, ihi, plscale, prscale,
305  n, p_balancing_mat2, n, info
306  F77_CHAR_ARG_LEN (1)
307  F77_CHAR_ARG_LEN (1)));
308 
309  return info;
310  }
311 
312  // Instantiations we need.
313 
314  template class gepbalance<Matrix>;
315 
316  template class gepbalance<FloatMatrix>;
317 
318  template class gepbalance<ComplexMatrix>;
319 
320  template class gepbalance<FloatComplexMatrix>;
321  }
322 }
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:41
a
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
octave::err_nonconformant
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-errwarn.cc:61
Matrix
Definition: dMatrix.h:36
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
octave::math::gepbalance::init
octave_idx_type init(const T &a, const T &b, const std::string &job)
b
b
Definition: cellfun.cc:400
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:41
octave_idx_type
i
for i
Definition: data.cc:5264
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
string
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:888