GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
gepbalance.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2017 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 the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 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 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include <string>
28 #include <vector>
29 
30 #include "Array-util.h"
31 #include "CMatrix.h"
32 #include "dMatrix.h"
33 #include "fCMatrix.h"
34 #include "fMatrix.h"
35 #include "gepbalance.h"
36 #include "lo-lapack-proto.h"
37 #include "oct-locbuf.h"
38 
39 namespace octave
40 {
41  namespace math
42  {
43  template <>
44  octave_idx_type
45  gepbalance<Matrix>::init (const Matrix& a, const Matrix& b,
46  const std::string& balance_job)
47  {
48  octave_idx_type n = a.cols ();
49 
50  if (a.rows () != n)
51  (*current_liboctave_error_handler) ("GEPBALANCE requires square matrix");
52 
53  if (a.dims () != b.dims ())
54  octave::err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
55 
56  octave_idx_type info;
57  octave_idx_type ilo;
58  octave_idx_type 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 (octave_idx_type 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  octave_idx_type n = 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  octave::err_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols());
120 
121  octave_idx_type info;
122  octave_idx_type ilo;
123  octave_idx_type ihi;
124 
125  OCTAVE_LOCAL_BUFFER (float, plscale, n);
126  OCTAVE_LOCAL_BUFFER (float, prscale, n);
127  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
128 
129  balanced_mat = a;
130  float *p_balanced_mat = balanced_mat.fortran_vec ();
131  balanced_mat2 = b;
132  float *p_balanced_mat2 = balanced_mat2.fortran_vec ();
133 
134  char job = balance_job[0];
135 
136  F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
137  n, p_balanced_mat, n, p_balanced_mat2,
138  n, ilo, ihi, plscale, prscale, pwork, info
139  F77_CHAR_ARG_LEN (1)));
140 
141  balancing_mat = FloatMatrix (n, n, 0.0);
142  balancing_mat2 = FloatMatrix (n, n, 0.0);
143  for (octave_idx_type i = 0; i < n; i++)
144  {
145  octave_quit ();
146  balancing_mat.elem (i ,i) = 1.0;
147  balancing_mat2.elem (i ,i) = 1.0;
148  }
149 
150  float *p_balancing_mat = balancing_mat.fortran_vec ();
151  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
152 
153  // first left
154  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
155  F77_CONST_CHAR_ARG2 ("L", 1),
156  n, ilo, ihi, plscale, prscale,
157  n, p_balancing_mat, n, info
158  F77_CHAR_ARG_LEN (1)
159  F77_CHAR_ARG_LEN (1)));
160 
161  // then right
162  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
163  F77_CONST_CHAR_ARG2 ("R", 1),
164  n, ilo, ihi, plscale, prscale,
165  n, p_balancing_mat2, n, info
166  F77_CHAR_ARG_LEN (1)
167  F77_CHAR_ARG_LEN (1)));
168 
169  return info;
170  }
171 
172  template <>
173  octave_idx_type
174  gepbalance<ComplexMatrix>::init (const ComplexMatrix& a,
175  const ComplexMatrix& b,
176  const std::string& balance_job)
177  {
178  octave_idx_type n = a.cols ();
179 
180  if (a.rows () != n)
181  (*current_liboctave_error_handler)
182  ("ComplexGEPBALANCE requires square matrix");
183 
184  if (a.dims () != b.dims ())
185  octave::err_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols());
186 
187  octave_idx_type info;
188  octave_idx_type ilo;
189  octave_idx_type ihi;
190 
191  OCTAVE_LOCAL_BUFFER (double, plscale, n);
192  OCTAVE_LOCAL_BUFFER (double, prscale, n);
193  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
194 
195  balanced_mat = a;
196  Complex *p_balanced_mat = balanced_mat.fortran_vec ();
197  balanced_mat2 = b;
198  Complex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
199 
200  char job = balance_job[0];
201 
202  F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
203  n, F77_DBLE_CMPLX_ARG (p_balanced_mat), n, F77_DBLE_CMPLX_ARG (p_balanced_mat2),
204  n, ilo, ihi, plscale, prscale, pwork, info
205  F77_CHAR_ARG_LEN (1)));
206 
207  balancing_mat = Matrix (n, n, 0.0);
208  balancing_mat2 = Matrix (n, n, 0.0);
209  for (octave_idx_type i = 0; i < n; i++)
210  {
211  octave_quit ();
212  balancing_mat.elem (i ,i) = 1.0;
213  balancing_mat2.elem (i ,i) = 1.0;
214  }
215 
216  double *p_balancing_mat = balancing_mat.fortran_vec ();
217  double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
218 
219  // first left
220  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
221  F77_CONST_CHAR_ARG2 ("L", 1),
222  n, ilo, ihi, plscale, prscale,
223  n, p_balancing_mat, n, info
224  F77_CHAR_ARG_LEN (1)
225  F77_CHAR_ARG_LEN (1)));
226 
227  // then right
228  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
229  F77_CONST_CHAR_ARG2 ("R", 1),
230  n, ilo, ihi, plscale, prscale,
231  n, p_balancing_mat2, n, info
232  F77_CHAR_ARG_LEN (1)
233  F77_CHAR_ARG_LEN (1)));
234 
235  return info;
236  }
237 
238  template <>
239  octave_idx_type
240  gepbalance<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
241  const FloatComplexMatrix& b,
242  const std::string& balance_job)
243  {
244  octave_idx_type n = a.cols ();
245 
246  if (a.rows () != n)
247  {
248  (*current_liboctave_error_handler)
249  ("FloatComplexGEPBALANCE requires square matrix");
250  return -1;
251  }
252 
253  if (a.dims () != b.dims ())
254  octave::err_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols());
255 
256  octave_idx_type info;
257  octave_idx_type ilo;
258  octave_idx_type ihi;
259 
260  OCTAVE_LOCAL_BUFFER (float, plscale, n);
261  OCTAVE_LOCAL_BUFFER (float, prscale, n);
262  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
263 
264  balanced_mat = a;
265  FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
266  balanced_mat2 = b;
267  FloatComplex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
268 
269  char job = balance_job[0];
270 
271  F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
272  n, F77_CMPLX_ARG (p_balanced_mat), n, F77_CMPLX_ARG (p_balanced_mat2),
273  n, ilo, ihi, plscale, prscale, pwork, info
274  F77_CHAR_ARG_LEN (1)));
275 
276  balancing_mat = FloatMatrix (n, n, 0.0);
277  balancing_mat2 = FloatMatrix (n, n, 0.0);
278  for (octave_idx_type i = 0; i < n; i++)
279  {
280  octave_quit ();
281  balancing_mat.elem (i ,i) = 1.0;
282  balancing_mat2.elem (i ,i) = 1.0;
283  }
284 
285  float *p_balancing_mat = balancing_mat.fortran_vec ();
286  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
287 
288  // first left
289  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
290  F77_CONST_CHAR_ARG2 ("L", 1),
291  n, ilo, ihi, plscale, prscale,
292  n, p_balancing_mat, n, info
293  F77_CHAR_ARG_LEN (1)
294  F77_CHAR_ARG_LEN (1)));
295 
296  // then right
297  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
298  F77_CONST_CHAR_ARG2 ("R", 1),
299  n, ilo, ihi, plscale, prscale,
300  n, p_balancing_mat2, n, info
301  F77_CHAR_ARG_LEN (1)
302  F77_CHAR_ARG_LEN (1)));
303 
304  return info;
305  }
306 
307  // Instantiations we need.
308 
309  template class gepbalance<Matrix>;
310 
311  template class gepbalance<FloatMatrix>;
312 
313  template class gepbalance<ComplexMatrix>;
314 
315  template class gepbalance<FloatComplexMatrix>;
316  }
317 }
octave
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
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:345
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:401
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:398
Array::dims
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:439
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:37
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
i
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave::math::gepbalance::init
octave_idx_type init(const T &a, const T &b, const std::string &job)
b
b
Definition: cellfun.cc:398
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
octave_idx_type
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:31
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:409
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:854
Generated on Wed May 10 2017 15:42:58 for GNU Octave by   doxygen 1.8.8