GNU Octave  8.1.0
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-2023 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include "CMatrix.h"
31 #include "dMatrix.h"
32 #include "fCMatrix.h"
33 #include "fMatrix.h"
34 #include "gepbalance.h"
35 #include "lo-array-errwarn.h"
36 #include "lo-error.h"
37 #include "lo-lapack-proto.h"
38 #include "oct-locbuf.h"
39 #include "quit.h"
40 
42 
44 
45 template <>
48  const std::string& balance_job)
49 {
50  F77_INT n = to_f77_int (a.cols ());
51 
52  if (a.rows () != n)
53  (*current_liboctave_error_handler)
54  ("GEPBALANCE requires square matrix");
55 
56  if (a.dims () != b.dims ())
57  err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
58 
59  F77_INT info;
60  F77_INT ilo;
61  F77_INT ihi;
62 
63  OCTAVE_LOCAL_BUFFER (double, plscale, n);
64  OCTAVE_LOCAL_BUFFER (double, prscale, n);
65  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
66 
67  m_balanced_mat = a;
68  double *p_balanced_mat = m_balanced_mat.fortran_vec ();
69  m_balanced_mat2 = b;
70  double *p_balanced_mat2 = m_balanced_mat2.fortran_vec ();
71 
72  char job = balance_job[0];
73 
74  F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
75  n, p_balanced_mat, n, p_balanced_mat2,
76  n, ilo, ihi, plscale, prscale, pwork, info
77  F77_CHAR_ARG_LEN (1)));
78 
79  m_balancing_mat = Matrix (n, n, 0.0);
80  m_balancing_mat2 = Matrix (n, n, 0.0);
81  for (F77_INT i = 0; i < n; i++)
82  {
83  octave_quit ();
84  m_balancing_mat.elem (i, i) = 1.0;
85  m_balancing_mat2.elem (i, i) = 1.0;
86  }
87 
88  double *p_balancing_mat = m_balancing_mat.fortran_vec ();
89  double *p_balancing_mat2 = m_balancing_mat2.fortran_vec ();
90 
91  // first left
92  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
93  F77_CONST_CHAR_ARG2 ("L", 1),
94  n, ilo, ihi, plscale, prscale,
95  n, p_balancing_mat, n, info
96  F77_CHAR_ARG_LEN (1)
97  F77_CHAR_ARG_LEN (1)));
98 
99  // then right
100  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
101  F77_CONST_CHAR_ARG2 ("R", 1),
102  n, ilo, ihi, plscale, prscale,
103  n, p_balancing_mat2, n, info
104  F77_CHAR_ARG_LEN (1)
105  F77_CHAR_ARG_LEN (1)));
106 
107  return info;
108 }
109 
110 template <>
113  const std::string& balance_job)
114 {
115  F77_INT n = to_f77_int (a.cols ());
116 
117  if (a.rows () != n)
118  (*current_liboctave_error_handler)
119  ("FloatGEPBALANCE requires square matrix");
120 
121  if (a.dims () != b.dims ())
122  err_nonconformant ("FloatGEPBALANCE",
123  n, n, b.rows(), b.cols());
124 
125  F77_INT info;
126  F77_INT ilo;
127  F77_INT ihi;
128 
129  OCTAVE_LOCAL_BUFFER (float, plscale, n);
130  OCTAVE_LOCAL_BUFFER (float, prscale, n);
131  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
132 
133  m_balanced_mat = a;
134  float *p_balanced_mat = m_balanced_mat.fortran_vec ();
135  m_balanced_mat2 = b;
136  float *p_balanced_mat2 = m_balanced_mat2.fortran_vec ();
137 
138  char job = balance_job[0];
139 
140  F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
141  n, p_balanced_mat, n, p_balanced_mat2,
142  n, ilo, ihi, plscale, prscale, pwork, info
143  F77_CHAR_ARG_LEN (1)));
144 
145  m_balancing_mat = FloatMatrix (n, n, 0.0);
146  m_balancing_mat2 = FloatMatrix (n, n, 0.0);
147  for (F77_INT i = 0; i < n; i++)
148  {
149  octave_quit ();
150  m_balancing_mat.elem (i, i) = 1.0;
151  m_balancing_mat2.elem (i, i) = 1.0;
152  }
153 
154  float *p_balancing_mat = m_balancing_mat.fortran_vec ();
155  float *p_balancing_mat2 = m_balancing_mat2.fortran_vec ();
156 
157  // first left
158  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
159  F77_CONST_CHAR_ARG2 ("L", 1),
160  n, ilo, ihi, plscale, prscale,
161  n, p_balancing_mat, n, info
162  F77_CHAR_ARG_LEN (1)
163  F77_CHAR_ARG_LEN (1)));
164 
165  // then right
166  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
167  F77_CONST_CHAR_ARG2 ("R", 1),
168  n, ilo, ihi, plscale, prscale,
169  n, p_balancing_mat2, n, info
170  F77_CHAR_ARG_LEN (1)
171  F77_CHAR_ARG_LEN (1)));
172 
173  return info;
174 }
175 
176 template <>
179  const ComplexMatrix& b,
180  const std::string& balance_job)
181 {
182  F77_INT n = to_f77_int (a.cols ());
183 
184  if (a.rows () != n)
185  (*current_liboctave_error_handler)
186  ("ComplexGEPBALANCE requires square matrix");
187 
188  if (a.dims () != b.dims ())
189  err_nonconformant ("ComplexGEPBALANCE",
190  n, n, b.rows(), b.cols());
191 
192  F77_INT info;
193  F77_INT ilo;
194  F77_INT ihi;
195 
196  OCTAVE_LOCAL_BUFFER (double, plscale, n);
197  OCTAVE_LOCAL_BUFFER (double, prscale, n);
198  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
199 
200  m_balanced_mat = a;
201  Complex *p_balanced_mat = m_balanced_mat.fortran_vec ();
202  m_balanced_mat2 = b;
203  Complex *p_balanced_mat2 = m_balanced_mat2.fortran_vec ();
204 
205  char job = balance_job[0];
206 
207  F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
208  n, F77_DBLE_CMPLX_ARG (p_balanced_mat),
209  n, F77_DBLE_CMPLX_ARG (p_balanced_mat2),
210  n, ilo, ihi, plscale, prscale, pwork, info
211  F77_CHAR_ARG_LEN (1)));
212 
213  m_balancing_mat = Matrix (n, n, 0.0);
214  m_balancing_mat2 = Matrix (n, n, 0.0);
215  for (F77_INT i = 0; i < n; i++)
216  {
217  octave_quit ();
218  m_balancing_mat.elem (i, i) = 1.0;
219  m_balancing_mat2.elem (i, i) = 1.0;
220  }
221 
222  double *p_balancing_mat = m_balancing_mat.fortran_vec ();
223  double *p_balancing_mat2 = m_balancing_mat2.fortran_vec ();
224 
225  // first left
226  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
227  F77_CONST_CHAR_ARG2 ("L", 1),
228  n, ilo, ihi, plscale, prscale,
229  n, p_balancing_mat, n, info
230  F77_CHAR_ARG_LEN (1)
231  F77_CHAR_ARG_LEN (1)));
232 
233  // then right
234  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
235  F77_CONST_CHAR_ARG2 ("R", 1),
236  n, ilo, ihi, plscale, prscale,
237  n, p_balancing_mat2, n, info
238  F77_CHAR_ARG_LEN (1)
239  F77_CHAR_ARG_LEN (1)));
240 
241  return info;
242 }
243 
244 template <>
247  const FloatComplexMatrix& b,
248  const std::string& balance_job)
249 {
250  F77_INT n = to_f77_int (a.cols ());
251 
252  if (a.rows () != n)
253  {
254  (*current_liboctave_error_handler)
255  ("FloatComplexGEPBALANCE requires square matrix");
256  return -1;
257  }
258 
259  if (a.dims () != b.dims ())
260  err_nonconformant ("FloatComplexGEPBALANCE",
261  n, n, b.rows(), b.cols());
262 
263  F77_INT info;
264  F77_INT ilo;
265  F77_INT ihi;
266 
267  OCTAVE_LOCAL_BUFFER (float, plscale, n);
268  OCTAVE_LOCAL_BUFFER (float, prscale, n);
269  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
270 
271  m_balanced_mat = a;
272  FloatComplex *p_balanced_mat = m_balanced_mat.fortran_vec ();
273  m_balanced_mat2 = b;
274  FloatComplex *p_balanced_mat2 = m_balanced_mat2.fortran_vec ();
275 
276  char job = balance_job[0];
277 
278  F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
279  n, F77_CMPLX_ARG (p_balanced_mat),
280  n, F77_CMPLX_ARG (p_balanced_mat2),
281  n, ilo, ihi, plscale, prscale, pwork, info
282  F77_CHAR_ARG_LEN (1)));
283 
284  m_balancing_mat = FloatMatrix (n, n, 0.0);
285  m_balancing_mat2 = FloatMatrix (n, n, 0.0);
286  for (F77_INT i = 0; i < n; i++)
287  {
288  octave_quit ();
289  m_balancing_mat.elem (i, i) = 1.0;
290  m_balancing_mat2.elem (i, i) = 1.0;
291  }
292 
293  float *p_balancing_mat = m_balancing_mat.fortran_vec ();
294  float *p_balancing_mat2 = m_balancing_mat2.fortran_vec ();
295 
296  // first left
297  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
298  F77_CONST_CHAR_ARG2 ("L", 1),
299  n, ilo, ihi, plscale, prscale,
300  n, p_balancing_mat, n, info
301  F77_CHAR_ARG_LEN (1)
302  F77_CHAR_ARG_LEN (1)));
303 
304  // then right
305  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
306  F77_CONST_CHAR_ARG2 ("R", 1),
307  n, ilo, ihi, plscale, prscale,
308  n, p_balancing_mat2, n, info
309  F77_CHAR_ARG_LEN (1)
310  F77_CHAR_ARG_LEN (1)));
311 
312  return info;
313 }
314 
315 // Instantiations we need.
316 
317 template class gepbalance<Matrix>;
318 
319 template class gepbalance<FloatMatrix>;
320 
321 template class gepbalance<ComplexMatrix>;
322 
323 template class gepbalance<FloatComplexMatrix>;
324 
OCTAVE_END_NAMESPACE(octave)
OCTARRAY_OVERRIDABLE_FUNC_API const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:503
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
Definition: Array.h:459
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type cols(void) const
Definition: Array.h:469
Definition: dMatrix.h:42
OCTAVE_API octave_idx_type init(const T &a, const T &b, const std::string &job)
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
#define OCTAVE_API
Definition: main.in.cc:55
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
class OCTAVE_API FloatMatrix
Definition: mx-fwd.h:33
octave_idx_type n
Definition: mx-inlines.cc:753
std::complex< double > Complex
Definition: oct-cmplx.h:33
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44