GNU Octave  4.2.1
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
aepbalance.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2017 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 the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 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 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <string>
29 
30 #include "CColVector.h"
31 #include "CMatrix.h"
32 #include "aepbalance.h"
33 #include "dColVector.h"
34 #include "dMatrix.h"
35 #include "fCColVector.h"
36 #include "fCMatrix.h"
37 #include "fColVector.h"
38 #include "fMatrix.h"
39 #include "lo-lapack-proto.h"
40 
41 static inline char
42 get_job (bool noperm, bool noscal)
43 {
44  return noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
45 }
46 
47 namespace octave
48 {
49  namespace math
50  {
51  template <>
52  aepbalance<Matrix>::aepbalance (const Matrix& a, bool noperm, bool noscal)
53  : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
54  {
55  octave_idx_type n = a.cols ();
56 
57  if (a.rows () != n)
58  (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
59 
60  scale = ColumnVector (n);
61 
62  octave_idx_type info;
63 
64  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
65  balanced_mat.fortran_vec (), n, ilo, ihi,
66  scale.fortran_vec (), info
67  F77_CHAR_ARG_LEN (1)));
68  }
69 
70  template <>
71  Matrix
72  aepbalance<Matrix>::balancing_matrix (void) const
73  {
74  octave_idx_type n = balanced_mat.rows ();
75  Matrix balancing_mat (n, n, 0.0);
76  for (octave_idx_type i = 0; i < n; i++)
77  balancing_mat.elem (i ,i) = 1.0;
78 
79  octave_idx_type info;
80 
81  char side = 'R';
82 
83  F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
84  F77_CONST_CHAR_ARG2 (&side, 1),
85  n, ilo, ihi, scale.data (), n,
86  balancing_mat.fortran_vec (), n, info
87  F77_CHAR_ARG_LEN (1)
88  F77_CHAR_ARG_LEN (1)));
89 
90  return balancing_mat;
91  }
92 
93  template <>
94  aepbalance<FloatMatrix>::aepbalance (const FloatMatrix& a, bool noperm,
95  bool noscal)
96  : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
97  {
98  octave_idx_type n = a.cols ();
99 
100  if (a.rows () != n)
101  (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
102 
103  scale = FloatColumnVector (n);
104 
105  octave_idx_type info;
106 
107  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
108  balanced_mat.fortran_vec (), n, ilo, ihi,
109  scale.fortran_vec (), info
110  F77_CHAR_ARG_LEN (1)));
111  }
112 
113  template <>
114  FloatMatrix
115  aepbalance<FloatMatrix>::balancing_matrix (void) const
116  {
117  octave_idx_type n = balanced_mat.rows ();
118  FloatMatrix balancing_mat (n, n, 0.0);
119  for (octave_idx_type i = 0; i < n; i++)
120  balancing_mat.elem (i ,i) = 1.0;
121 
122  octave_idx_type info;
123 
124  char side = 'R';
125 
126  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
127  F77_CONST_CHAR_ARG2 (&side, 1),
128  n, ilo, ihi, scale.data (), n,
129  balancing_mat.fortran_vec (), n, info
130  F77_CHAR_ARG_LEN (1)
131  F77_CHAR_ARG_LEN (1)));
132 
133  return balancing_mat;
134  }
135 
136  template <>
137  aepbalance<ComplexMatrix>::aepbalance (const ComplexMatrix& a, bool noperm,
138  bool noscal)
139  : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
140  {
141  octave_idx_type n = a.cols ();
142 
143  if (a.rows () != n)
144  (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
145 
146  scale = ColumnVector (n);
147 
148  octave_idx_type info;
149 
150  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
151  F77_DBLE_CMPLX_ARG (balanced_mat.fortran_vec ()), n, ilo, ihi,
152  scale.fortran_vec (), info
153  F77_CHAR_ARG_LEN (1)));
154  }
155 
156  template <>
157  ComplexMatrix
158  aepbalance<ComplexMatrix>::balancing_matrix (void) const
159  {
160  octave_idx_type n = balanced_mat.rows ();
161  ComplexMatrix balancing_mat (n, n, 0.0);
162  for (octave_idx_type i = 0; i < n; i++)
163  balancing_mat.elem (i, i) = 1.0;
164 
165  octave_idx_type info;
166 
167  char side = 'R';
168 
169  F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
170  F77_CONST_CHAR_ARG2 (&side, 1),
171  n, ilo, ihi, scale.data (), n,
172  F77_DBLE_CMPLX_ARG (balancing_mat.fortran_vec ()), n, info
173  F77_CHAR_ARG_LEN (1)
174  F77_CHAR_ARG_LEN (1)));
175 
176  return balancing_mat;
177  }
178 
179  template <>
180  aepbalance<FloatComplexMatrix>::aepbalance (const FloatComplexMatrix& a,
181  bool noperm, bool noscal)
182  : balanced_mat (a), scale (), ilo (), ihi (), job (get_job (noperm, noscal))
183  {
184  octave_idx_type n = a.cols ();
185 
186  if (a.rows () != n)
187  (*current_liboctave_error_handler) ("aepbalance: requires square matrix");
188 
189  scale = FloatColumnVector (n);
190 
191  octave_idx_type info;
192 
193  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), n,
194  F77_CMPLX_ARG (balanced_mat.fortran_vec ()), n, ilo, ihi,
195  scale.fortran_vec (), info
196  F77_CHAR_ARG_LEN (1)));
197  }
198 
199  template <>
200  FloatComplexMatrix
201  aepbalance<FloatComplexMatrix>::balancing_matrix (void) const
202  {
203  octave_idx_type n = balanced_mat.rows ();
204  FloatComplexMatrix balancing_mat (n, n, 0.0);
205  for (octave_idx_type i = 0; i < n; i++)
206  balancing_mat.elem (i, i) = 1.0;
207 
208  octave_idx_type info;
209 
210  char side = 'R';
211 
212  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
213  F77_CONST_CHAR_ARG2 (&side, 1),
214  n, ilo, ihi, scale.data (), n,
215  F77_CMPLX_ARG (balancing_mat.fortran_vec ()), n, info
216  F77_CHAR_ARG_LEN (1)
217  F77_CHAR_ARG_LEN (1)));
218 
219  return balancing_mat;
220  }
221 
222  // Instantiations we need.
223 
224  template class aepbalance<Matrix>;
225 
226  template class aepbalance<FloatMatrix>;
227 
228  template class aepbalance<ComplexMatrix>;
229 
230  template class aepbalance<FloatComplexMatrix>;
231  }
232 }
octave
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
get_job
static char get_job(bool noperm, bool noscal)
Definition: aepbalance.cc:42
Array::elem
T & elem(octave_idx_type n)
Definition: Array.h:482
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:401
a
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:398
aepbalance::aepbalance
aepbalance(const Matrix &a, bool noperm, bool noscal)
Definition: aepbalance.cc:52
Matrix
Definition: dMatrix.h:37
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
aepbalance::balancing_matrix
Matrix balancing_matrix(void) const
Definition: aepbalance.cc:72
i
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
scale
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
octave_idx_type
Array::fortran_vec
const T * fortran_vec(void) const
Definition: Array.h:584
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:409
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8