GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
aepbalance.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2022 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 "aepbalance.h"
32#include "dColVector.h"
33#include "dMatrix.h"
34#include "fCMatrix.h"
35#include "fColVector.h"
36#include "fMatrix.h"
37#include "lo-error.h"
38#include "lo-lapack-proto.h"
39
40namespace octave
41{
42 static inline char
43 get_job (bool noperm, bool noscal)
44 {
45 return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
46 }
47
48 namespace math
49 {
50 template <>
51 OCTAVE_API
52 aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal)
53 : m_balanced_mat (a), m_scale (), m_ilo (), m_ihi (),
54 m_job (get_job (noperm, noscal))
55 {
56 F77_INT n = to_f77_int (a.cols ());
57
58 if (a.rows () != n)
59 (*current_liboctave_error_handler)
60 ("aepbalance: requires square matrix");
61
62 m_scale = ColumnVector (n);
63
64 F77_INT info, t_ilo, t_ihi;
65
66 F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&m_job, 1), n,
67 m_balanced_mat.fortran_vec (), n,
68 t_ilo, t_ihi, m_scale.fortran_vec (), info
69 F77_CHAR_ARG_LEN (1)));
70
71 m_ilo = t_ilo;
72 m_ihi = t_ihi;
73 }
74
75 template <>
76 OCTAVE_API Matrix
77 aepbalance<Matrix>::balancing_matrix (void) const
78 {
79 F77_INT n = to_f77_int (m_balanced_mat.rows ());
80
81 Matrix balancing_mat (n, n, 0.0);
82 for (F77_INT i = 0; i < n; i++)
83 balancing_mat.elem (i, i) = 1.0;
84
85 F77_INT info;
86 F77_INT t_ilo = to_f77_int (m_ilo);
87 F77_INT t_ihi = to_f77_int (m_ihi);
88
89 char side = 'R';
90
91 F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&m_job, 1),
92 F77_CONST_CHAR_ARG2 (&side, 1),
93 n, t_ilo, t_ihi, m_scale.data (), n,
94 balancing_mat.fortran_vec (), n, info
95 F77_CHAR_ARG_LEN (1)
96 F77_CHAR_ARG_LEN (1)));
97
98 return balancing_mat;
99 }
100
101 template <>
102 OCTAVE_API
103 aepbalance<FloatMatrix>::aepbalance (const FloatMatrix& a, bool noperm,
104 bool noscal)
105 : m_balanced_mat (a), m_scale (), m_ilo (), m_ihi (),
106 m_job (get_job (noperm, noscal))
107 {
108 F77_INT n = to_f77_int (a.cols ());
109
110 if (a.rows () != n)
111 (*current_liboctave_error_handler)
112 ("aepbalance: requires square matrix");
113
114 m_scale = FloatColumnVector (n);
115
116 F77_INT info, t_ilo, t_ihi;
117
118 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&m_job, 1), n,
119 m_balanced_mat.fortran_vec (), n, t_ilo,
120 t_ihi, m_scale.fortran_vec (), info
121 F77_CHAR_ARG_LEN (1)));
122
123 m_ilo = t_ilo;
124 m_ihi = t_ihi;
125 }
126
127 template <>
128 OCTAVE_API FloatMatrix
129 aepbalance<FloatMatrix>::balancing_matrix (void) const
130 {
131 F77_INT n = to_f77_int (m_balanced_mat.rows ());
132
133 FloatMatrix balancing_mat (n, n, 0.0);
134 for (F77_INT i = 0; i < n; i++)
135 balancing_mat.elem (i, i) = 1.0;
136
137 F77_INT info;
138 F77_INT t_ilo = to_f77_int (m_ilo);
139 F77_INT t_ihi = to_f77_int (m_ihi);
140
141 char side = 'R';
142
143 F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&m_job, 1),
144 F77_CONST_CHAR_ARG2 (&side, 1),
145 n, t_ilo, t_ihi, m_scale.data (), n,
146 balancing_mat.fortran_vec (), n, info
147 F77_CHAR_ARG_LEN (1)
148 F77_CHAR_ARG_LEN (1)));
149
150 return balancing_mat;
151 }
152
153 template <>
154 OCTAVE_API
155 aepbalance<ComplexMatrix>::aepbalance (const ComplexMatrix& a, bool noperm,
156 bool noscal)
157 : m_balanced_mat (a), m_scale (), m_ilo (), m_ihi (),
158 m_job (get_job (noperm, noscal))
159 {
160 F77_INT n = to_f77_int (a.cols ());
161
162 if (a.rows () != n)
163 (*current_liboctave_error_handler)
164 ("aepbalance: requires square matrix");
165
166 m_scale = ColumnVector (n);
167
168 F77_INT info, t_ilo, t_ihi;
169
170 F77_XFCN (zgebal, ZGEBAL,
171 (F77_CONST_CHAR_ARG2 (&m_job, 1), n,
172 F77_DBLE_CMPLX_ARG (m_balanced_mat.fortran_vec ()),
173 n, t_ilo, t_ihi, m_scale.fortran_vec (), info
174 F77_CHAR_ARG_LEN (1)));
175
176 m_ilo = t_ilo;
177 m_ihi = t_ihi;
178 }
179
180 template <>
181 OCTAVE_API ComplexMatrix
182 aepbalance<ComplexMatrix>::balancing_matrix (void) const
183 {
184 F77_INT n = to_f77_int (m_balanced_mat.rows ());
185
186 ComplexMatrix balancing_mat (n, n, 0.0);
187 for (F77_INT i = 0; i < n; i++)
188 balancing_mat.elem (i, i) = 1.0;
189
190 F77_INT info;
191 F77_INT t_ilo = to_f77_int (m_ilo);
192 F77_INT t_ihi = to_f77_int (m_ihi);
193
194 char side = 'R';
195
196 F77_XFCN (zgebak, ZGEBAK,
197 (F77_CONST_CHAR_ARG2 (&m_job, 1),
198 F77_CONST_CHAR_ARG2 (&side, 1),
199 n, t_ilo, t_ihi, m_scale.data (), n,
200 F77_DBLE_CMPLX_ARG (balancing_mat.fortran_vec ()),
201 n, info
202 F77_CHAR_ARG_LEN (1)
203 F77_CHAR_ARG_LEN (1)));
204
205 return balancing_mat;
206 }
207
208 template <>
209 OCTAVE_API
210 aepbalance<FloatComplexMatrix>::aepbalance (const FloatComplexMatrix& a,
211 bool noperm, bool noscal)
212 : m_balanced_mat (a), m_scale (), m_ilo (), m_ihi (),
213 m_job (get_job (noperm, noscal))
214 {
215 F77_INT n = to_f77_int (a.cols ());
216
217 if (a.rows () != n)
218 (*current_liboctave_error_handler)
219 ("aepbalance: requires square matrix");
220
221 m_scale = FloatColumnVector (n);
222
223 F77_INT info, t_ilo, t_ihi;
224
225 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&m_job, 1), n,
226 F77_CMPLX_ARG (m_balanced_mat.fortran_vec ()),
227 n, t_ilo, t_ihi, m_scale.fortran_vec (), info
228 F77_CHAR_ARG_LEN (1)));
229
230 m_ilo = t_ilo;
231 m_ihi = t_ihi;
232 }
233
234 template <>
235 OCTAVE_API FloatComplexMatrix
236 aepbalance<FloatComplexMatrix>::balancing_matrix (void) const
237 {
238 F77_INT n = to_f77_int (m_balanced_mat.rows ());
239
240 FloatComplexMatrix balancing_mat (n, n, 0.0);
241 for (F77_INT i = 0; i < n; i++)
242 balancing_mat.elem (i, i) = 1.0;
243
244 F77_INT info;
245 F77_INT t_ilo = to_f77_int (m_ilo);
246 F77_INT t_ihi = to_f77_int (m_ihi);
247
248 char side = 'R';
249
250 F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&m_job, 1),
251 F77_CONST_CHAR_ARG2 (&side, 1),
252 n, t_ilo, t_ihi, m_scale.data (), n,
253 F77_CMPLX_ARG (balancing_mat.fortran_vec ()),
254 n, info
255 F77_CHAR_ARG_LEN (1)
256 F77_CHAR_ARG_LEN (1)));
257
258 return balancing_mat;
259 }
260
261 // Instantiations we need.
262
263 template class aepbalance<Matrix>;
264
265 template class aepbalance<FloatMatrix>;
266
267 template class aepbalance<ComplexMatrix>;
268
269 template class aepbalance<FloatComplexMatrix>;
270 }
271}
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:457
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:534
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:449
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
Matrix
Definition: dMatrix.h:42
aepbalance::aepbalance
OCTAVE_API aepbalance(const Matrix &a, bool noperm, bool noscal)
Definition: aepbalance.cc:52
aepbalance::balancing_matrix
OCTAVE_API Matrix balancing_matrix(void) const
Definition: aepbalance.cc:77
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
FloatColumnVector
class OCTAVE_API FloatColumnVector
Definition: mx-fwd.h:47
ColumnVector
class OCTAVE_API ColumnVector
Definition: mx-fwd.h:45
octave::get_job
static char get_job(bool noperm, bool noscal)
Definition: aepbalance.cc:43