GNU Octave  4.4.1
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-2018 John W. Eaton
4 Copyright (C) 2008 Jaroslav Hajek
5 
6 This file is part of Octave.
7 
8 Octave is free software: you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "CMatrix.h"
29 #include "aepbalance.h"
30 #include "dColVector.h"
31 #include "dMatrix.h"
32 #include "fCMatrix.h"
33 #include "fColVector.h"
34 #include "fMatrix.h"
35 #include "lo-error.h"
36 #include "lo-lapack-proto.h"
37 
38 static inline char
39 get_job (bool noperm, bool noscal)
40 {
41  return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
42 }
43 
44 namespace octave
45 {
46  namespace math
47  {
48  template <>
49  aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal)
50  : balanced_mat (a), scale (), ilo (), ihi (),
51  job (get_job (noperm, noscal))
52  {
53  F77_INT n = to_f77_int (a.cols ());
54 
55  if (a.rows () != n)
57  ("aepbalance: requires square matrix");
58 
59  scale = ColumnVector (n);
60 
61  F77_INT info, t_ilo, t_ihi;
62 
63  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
64  balanced_mat.fortran_vec (), n,
65  t_ilo, t_ihi, scale.fortran_vec (), info
66  F77_CHAR_ARG_LEN (1)));
67 
68  ilo = t_ilo;
69  ihi = t_ihi;
70  }
71 
72  template <>
73  Matrix
75  {
76  F77_INT n = to_f77_int (balanced_mat.rows ());
77 
78  Matrix balancing_mat (n, n, 0.0);
79  for (F77_INT i = 0; i < n; i++)
80  balancing_mat.elem (i ,i) = 1.0;
81 
82  F77_INT info;
83  F77_INT t_ilo = to_f77_int (ilo);
84  F77_INT t_ihi = to_f77_int (ihi);
85 
86  char side = 'R';
87 
88  F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
89  F77_CONST_CHAR_ARG2 (&side, 1),
90  n, t_ilo, t_ihi, scale.data (), n,
91  balancing_mat.fortran_vec (), n, info
92  F77_CHAR_ARG_LEN (1)
93  F77_CHAR_ARG_LEN (1)));
94 
95  return balancing_mat;
96  }
97 
98  template <>
100  bool noscal)
101  : balanced_mat (a), scale (), ilo (), ihi (),
102  job (get_job (noperm, noscal))
103  {
104  F77_INT n = to_f77_int (a.cols ());
105 
106  if (a.rows () != n)
108  ("aepbalance: requires square matrix");
109 
110  scale = FloatColumnVector (n);
111 
112  F77_INT info, t_ilo, t_ihi;
113 
114  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
115  balanced_mat.fortran_vec (), n, t_ilo,
116  t_ihi, scale.fortran_vec (), info
117  F77_CHAR_ARG_LEN (1)));
118 
119  ilo = t_ilo;
120  ihi = t_ihi;
121  }
122 
123  template <>
126  {
127  F77_INT n = to_f77_int (balanced_mat.rows ());
128 
129  FloatMatrix balancing_mat (n, n, 0.0);
130  for (F77_INT i = 0; i < n; i++)
131  balancing_mat.elem (i ,i) = 1.0;
132 
133  F77_INT info;
134  F77_INT t_ilo = to_f77_int (ilo);
135  F77_INT t_ihi = to_f77_int (ihi);
136 
137  char side = 'R';
138 
139  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
140  F77_CONST_CHAR_ARG2 (&side, 1),
141  n, t_ilo, t_ihi, scale.data (), n,
142  balancing_mat.fortran_vec (), n, info
143  F77_CHAR_ARG_LEN (1)
144  F77_CHAR_ARG_LEN (1)));
145 
146  return balancing_mat;
147  }
148 
149  template <>
151  bool noscal)
152  : balanced_mat (a), scale (), ilo (), ihi (),
153  job (get_job (noperm, noscal))
154  {
155  F77_INT n = to_f77_int (a.cols ());
156 
157  if (a.rows () != n)
159  ("aepbalance: requires square matrix");
160 
161  scale = ColumnVector (n);
162 
163  F77_INT info, t_ilo, t_ihi;
164 
165  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
166  F77_DBLE_CMPLX_ARG (balanced_mat.fortran_vec ()),
167  n, t_ilo, t_ihi, scale.fortran_vec (), info
168  F77_CHAR_ARG_LEN (1)));
169 
170  ilo = t_ilo;
171  ihi = t_ihi;
172  }
173 
174  template <>
177  {
178  F77_INT n = to_f77_int (balanced_mat.rows ());
179 
180  ComplexMatrix balancing_mat (n, n, 0.0);
181  for (F77_INT i = 0; i < n; i++)
182  balancing_mat.elem (i, i) = 1.0;
183 
184  F77_INT info;
185  F77_INT t_ilo = to_f77_int (ilo);
186  F77_INT t_ihi = to_f77_int (ihi);
187 
188  char side = 'R';
189 
190  F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
191  F77_CONST_CHAR_ARG2 (&side, 1),
192  n, t_ilo, t_ihi, scale.data (), n,
193  F77_DBLE_CMPLX_ARG (balancing_mat.fortran_vec ()),
194  n, info
195  F77_CHAR_ARG_LEN (1)
196  F77_CHAR_ARG_LEN (1)));
197 
198  return balancing_mat;
199  }
200 
201  template <>
203  bool noperm, bool noscal)
204  : balanced_mat (a), scale (), ilo (), ihi (),
205  job (get_job (noperm, noscal))
206  {
207  F77_INT n = to_f77_int (a.cols ());
208 
209  if (a.rows () != n)
211  ("aepbalance: requires square matrix");
212 
213  scale = FloatColumnVector (n);
214 
215  F77_INT info, t_ilo, t_ihi;
216 
217  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
218  F77_CMPLX_ARG (balanced_mat.fortran_vec ()),
219  n, t_ilo, t_ihi, scale.fortran_vec (), info
220  F77_CHAR_ARG_LEN (1)));
221 
222  ilo = t_ilo;
223  ihi = t_ihi;
224  }
225 
226  template <>
229  {
230  F77_INT n = to_f77_int (balanced_mat.rows ());
231 
232  FloatComplexMatrix balancing_mat (n, n, 0.0);
233  for (F77_INT i = 0; i < n; i++)
234  balancing_mat.elem (i, i) = 1.0;
235 
236  F77_INT info;
237  F77_INT t_ilo = to_f77_int (ilo);
238  F77_INT t_ihi = to_f77_int (ihi);
239 
240  char side = 'R';
241 
242  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
243  F77_CONST_CHAR_ARG2 (&side, 1),
244  n, t_ilo, t_ihi, scale.data (), n,
245  F77_CMPLX_ARG (balancing_mat.fortran_vec ()),
246  n, info
247  F77_CHAR_ARG_LEN (1)
248  F77_CHAR_ARG_LEN (1)));
249 
250  return balancing_mat;
251  }
252 
253  // Instantiations we need.
254 
255  template class aepbalance<Matrix>;
256 
257  template class aepbalance<FloatMatrix>;
258 
259  template class aepbalance<ComplexMatrix>;
260 
261  template class aepbalance<FloatComplexMatrix>;
262  }
263 }
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
static char get_job(bool noperm, bool noscal)
Definition: aepbalance.cc:39
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:41
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:400
aepbalance(const Matrix &a, bool noperm, bool noscal)
Definition: aepbalance.cc:49
Definition: dMatrix.h:36
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
Matrix balancing_matrix(void) const
Definition: aepbalance.cc:74
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5442
for i
Definition: data.cc:5264