GNU Octave 11.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-2026 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 "lapack-proto.h"
36#include "lo-array-errwarn.h"
37#include "oct-error.h"
38#include "oct-locbuf.h"
39#include "quit.h"
40
41OCTAVE_BEGIN_NAMESPACE(octave)
42OCTAVE_BEGIN_NAMESPACE(math)
43
44template <>
45OCTAVE_API octave_idx_type
46gepbalance<Matrix>::init (const Matrix& a, const Matrix& b,
47 const std::string& balance_job)
48{
49 F77_INT n = to_f77_int (a.cols ());
50
51 if (a.rows () != n)
52 (*current_liboctave_error_handler)
53 ("GEPBALANCE requires square matrix");
54
55 if (a.dims () != b.dims ())
56 err_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
57
58 F77_INT info;
59 F77_INT ilo;
60 F77_INT ihi;
61
62 OCTAVE_LOCAL_BUFFER (double, plscale, n);
63 OCTAVE_LOCAL_BUFFER (double, prscale, n);
64 OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
65
66 m_balanced_mat = a;
67 double *p_balanced_mat = m_balanced_mat.rwdata ();
68 m_balanced_mat2 = b;
69 double *p_balanced_mat2 = m_balanced_mat2.rwdata ();
70
71 char job = balance_job[0];
72
73 F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
74 n, p_balanced_mat, n, p_balanced_mat2,
75 n, ilo, ihi, plscale, prscale, pwork, info
76 F77_CHAR_ARG_LEN (1)));
77
78 m_balancing_mat = Matrix (n, n, 0.0);
79 m_balancing_mat2 = Matrix (n, n, 0.0);
80 for (F77_INT i = 0; i < n; i++)
81 {
82 octave_quit ();
83 m_balancing_mat.elem (i, i) = 1.0;
84 m_balancing_mat2.elem (i, i) = 1.0;
85 }
86
87 double *p_balancing_mat = m_balancing_mat.rwdata ();
88 double *p_balancing_mat2 = m_balancing_mat2.rwdata ();
89
90 // first left
91 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
92 F77_CONST_CHAR_ARG2 ("L", 1),
93 n, ilo, ihi, plscale, prscale,
94 n, p_balancing_mat, n, info
95 F77_CHAR_ARG_LEN (1)
96 F77_CHAR_ARG_LEN (1)));
97
98 // then right
99 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
100 F77_CONST_CHAR_ARG2 ("R", 1),
101 n, ilo, ihi, plscale, prscale,
102 n, p_balancing_mat2, n, info
103 F77_CHAR_ARG_LEN (1)
104 F77_CHAR_ARG_LEN (1)));
105
106 return info;
107}
108
109template <>
110OCTAVE_API octave_idx_type
111gepbalance<FloatMatrix>::init (const FloatMatrix& a, const FloatMatrix& b,
112 const std::string& balance_job)
113{
114 F77_INT n = to_f77_int (a.cols ());
115
116 if (a.rows () != n)
117 (*current_liboctave_error_handler)
118 ("FloatGEPBALANCE requires square matrix");
119
120 if (a.dims () != b.dims ())
121 err_nonconformant ("FloatGEPBALANCE",
122 n, n, b.rows(), b.cols());
123
124 F77_INT info;
125 F77_INT ilo;
126 F77_INT ihi;
127
128 OCTAVE_LOCAL_BUFFER (float, plscale, n);
129 OCTAVE_LOCAL_BUFFER (float, prscale, n);
130 OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
131
132 m_balanced_mat = a;
133 float *p_balanced_mat = m_balanced_mat.rwdata ();
134 m_balanced_mat2 = b;
135 float *p_balanced_mat2 = m_balanced_mat2.rwdata ();
136
137 char job = balance_job[0];
138
139 F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
140 n, p_balanced_mat, n, p_balanced_mat2,
141 n, ilo, ihi, plscale, prscale, pwork, info
142 F77_CHAR_ARG_LEN (1)));
143
144 m_balancing_mat = FloatMatrix (n, n, 0.0);
145 m_balancing_mat2 = FloatMatrix (n, n, 0.0);
146 for (F77_INT i = 0; i < n; i++)
147 {
148 octave_quit ();
149 m_balancing_mat.elem (i, i) = 1.0;
150 m_balancing_mat2.elem (i, i) = 1.0;
151 }
152
153 float *p_balancing_mat = m_balancing_mat.rwdata ();
154 float *p_balancing_mat2 = m_balancing_mat2.rwdata ();
155
156 // first left
157 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
158 F77_CONST_CHAR_ARG2 ("L", 1),
159 n, ilo, ihi, plscale, prscale,
160 n, p_balancing_mat, n, info
161 F77_CHAR_ARG_LEN (1)
162 F77_CHAR_ARG_LEN (1)));
163
164 // then right
165 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
166 F77_CONST_CHAR_ARG2 ("R", 1),
167 n, ilo, ihi, plscale, prscale,
168 n, p_balancing_mat2, n, info
169 F77_CHAR_ARG_LEN (1)
170 F77_CHAR_ARG_LEN (1)));
171
172 return info;
173}
174
175template <>
176OCTAVE_API octave_idx_type
177gepbalance<ComplexMatrix>::init (const ComplexMatrix& a,
178 const ComplexMatrix& b,
179 const std::string& balance_job)
180{
181 F77_INT n = to_f77_int (a.cols ());
182
183 if (a.rows () != n)
184 (*current_liboctave_error_handler)
185 ("ComplexGEPBALANCE requires square matrix");
186
187 if (a.dims () != b.dims ())
188 err_nonconformant ("ComplexGEPBALANCE",
189 n, n, b.rows(), b.cols());
190
191 F77_INT info;
192 F77_INT ilo;
193 F77_INT ihi;
194
195 OCTAVE_LOCAL_BUFFER (double, plscale, n);
196 OCTAVE_LOCAL_BUFFER (double, prscale, n);
197 OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
198
199 m_balanced_mat = a;
200 Complex *p_balanced_mat = m_balanced_mat.rwdata ();
201 m_balanced_mat2 = b;
202 Complex *p_balanced_mat2 = m_balanced_mat2.rwdata ();
203
204 char job = balance_job[0];
205
206 F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
207 n, F77_DBLE_CMPLX_ARG (p_balanced_mat),
208 n, F77_DBLE_CMPLX_ARG (p_balanced_mat2),
209 n, ilo, ihi, plscale, prscale, pwork, info
210 F77_CHAR_ARG_LEN (1)));
211
212 m_balancing_mat = Matrix (n, n, 0.0);
213 m_balancing_mat2 = Matrix (n, n, 0.0);
214 for (F77_INT i = 0; i < n; i++)
215 {
216 octave_quit ();
217 m_balancing_mat.elem (i, i) = 1.0;
218 m_balancing_mat2.elem (i, i) = 1.0;
219 }
220
221 double *p_balancing_mat = m_balancing_mat.rwdata ();
222 double *p_balancing_mat2 = m_balancing_mat2.rwdata ();
223
224 // first left
225 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
226 F77_CONST_CHAR_ARG2 ("L", 1),
227 n, ilo, ihi, plscale, prscale,
228 n, p_balancing_mat, n, info
229 F77_CHAR_ARG_LEN (1)
230 F77_CHAR_ARG_LEN (1)));
231
232 // then right
233 F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
234 F77_CONST_CHAR_ARG2 ("R", 1),
235 n, ilo, ihi, plscale, prscale,
236 n, p_balancing_mat2, n, info
237 F77_CHAR_ARG_LEN (1)
238 F77_CHAR_ARG_LEN (1)));
239
240 return info;
241}
242
243template <>
244OCTAVE_API octave_idx_type
245gepbalance<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
246 const FloatComplexMatrix& b,
247 const std::string& balance_job)
248{
249 F77_INT n = to_f77_int (a.cols ());
250
251 if (a.rows () != n)
252 {
253 (*current_liboctave_error_handler)
254 ("FloatComplexGEPBALANCE requires square matrix");
255 return -1;
256 }
257
258 if (a.dims () != b.dims ())
259 err_nonconformant ("FloatComplexGEPBALANCE",
260 n, n, b.rows(), b.cols());
261
262 F77_INT info;
263 F77_INT ilo;
264 F77_INT ihi;
265
266 OCTAVE_LOCAL_BUFFER (float, plscale, n);
267 OCTAVE_LOCAL_BUFFER (float, prscale, n);
268 OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
269
270 m_balanced_mat = a;
271 FloatComplex *p_balanced_mat = m_balanced_mat.rwdata ();
272 m_balanced_mat2 = b;
273 FloatComplex *p_balanced_mat2 = m_balanced_mat2.rwdata ();
274
275 char job = balance_job[0];
276
277 F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
278 n, F77_CMPLX_ARG (p_balanced_mat),
279 n, F77_CMPLX_ARG (p_balanced_mat2),
280 n, ilo, ihi, plscale, prscale, pwork, info
281 F77_CHAR_ARG_LEN (1)));
282
283 m_balancing_mat = FloatMatrix (n, n, 0.0);
284 m_balancing_mat2 = FloatMatrix (n, n, 0.0);
285 for (F77_INT i = 0; i < n; i++)
286 {
287 octave_quit ();
288 m_balancing_mat.elem (i, i) = 1.0;
289 m_balancing_mat2.elem (i, i) = 1.0;
290 }
291
292 float *p_balancing_mat = m_balancing_mat.rwdata ();
293 float *p_balancing_mat2 = m_balancing_mat2.rwdata ();
294
295 // first left
296 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
297 F77_CONST_CHAR_ARG2 ("L", 1),
298 n, ilo, ihi, plscale, prscale,
299 n, p_balancing_mat, n, info
300 F77_CHAR_ARG_LEN (1)
301 F77_CHAR_ARG_LEN (1)));
302
303 // then right
304 F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
305 F77_CONST_CHAR_ARG2 ("R", 1),
306 n, ilo, ihi, plscale, prscale,
307 n, p_balancing_mat2, n, info
308 F77_CHAR_ARG_LEN (1)
309 F77_CHAR_ARG_LEN (1)));
310
311 return info;
312}
313
314// Instantiations we need.
315
316template class gepbalance<Matrix>;
317
318template class gepbalance<FloatMatrix>;
319
320template class gepbalance<ComplexMatrix>;
321
322template class gepbalance<FloatComplexMatrix>;
323
324OCTAVE_END_NAMESPACE(math)
325OCTAVE_END_NAMESPACE(octave)
Array::dims
const dim_vector & dims() const
Return a const-reference so that dims ()(i) works efficiently.
Definition Array-base.h:529
Array::rows
octave_idx_type rows() const
Definition Array-base.h:485
Array::cols
octave_idx_type cols() const
Definition Array-base.h:495
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1790
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