GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
gepbalance.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2025 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
41OCTAVE_BEGIN_NAMESPACE(octave)
42
43OCTAVE_BEGIN_NAMESPACE(math)
44
45template <>
46OCTAVE_API octave_idx_type
47gepbalance<Matrix>::init (const Matrix& a, const Matrix& b,
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.rwdata ();
69 m_balanced_mat2 = b;
70 double *p_balanced_mat2 = m_balanced_mat2.rwdata ();
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.rwdata ();
89 double *p_balancing_mat2 = m_balancing_mat2.rwdata ();
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
110template <>
111OCTAVE_API octave_idx_type
112gepbalance<FloatMatrix>::init (const FloatMatrix& a, const FloatMatrix& b,
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.rwdata ();
135 m_balanced_mat2 = b;
136 float *p_balanced_mat2 = m_balanced_mat2.rwdata ();
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.rwdata ();
155 float *p_balancing_mat2 = m_balancing_mat2.rwdata ();
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
176template <>
177OCTAVE_API octave_idx_type
178gepbalance<ComplexMatrix>::init (const ComplexMatrix& a,
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.rwdata ();
202 m_balanced_mat2 = b;
203 Complex *p_balanced_mat2 = m_balanced_mat2.rwdata ();
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.rwdata ();
223 double *p_balancing_mat2 = m_balancing_mat2.rwdata ();
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
244template <>
245OCTAVE_API octave_idx_type
246gepbalance<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
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.rwdata ();
273 m_balanced_mat2 = b;
274 FloatComplex *p_balanced_mat2 = m_balanced_mat2.rwdata ();
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.rwdata ();
294 float *p_balancing_mat2 = m_balancing_mat2.rwdata ();
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
317template class gepbalance<Matrix>;
318
319template class gepbalance<FloatMatrix>;
320
321template class gepbalance<ComplexMatrix>;
322
323template class gepbalance<FloatComplexMatrix>;
324
325OCTAVE_END_NAMESPACE(math)
326OCTAVE_END_NAMESPACE(octave)
Array::dims
const dim_vector & dims() const
Return a const-reference so that dims ()(i) works efficiently.
Definition Array.h:507
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
octave_idx_type
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
err_nonconformant
void err_nonconformant()
Definition errwarn.cc:95
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
OCTAVE_API
#define OCTAVE_API
Definition main.in.cc:55
Complex
std::complex< double > Complex
Definition oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44