GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
balance.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2018 John W. Eaton
4 Copyright (C) 2008-2009 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 // Author: A. S. Hodel <scotte@eng.auburn.edu>
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <string>
31 
32 #include "CMatrix.h"
33 #include "aepbalance.h"
34 #include "dMatrix.h"
35 #include "fCMatrix.h"
36 #include "fMatrix.h"
37 #include "gepbalance.h"
38 #include "quit.h"
39 
40 #include "defun.h"
41 #include "error.h"
42 #include "f77-fcn.h"
43 #include "errwarn.h"
44 #include "ovl.h"
45 #include "utils.h"
46 
47 DEFUN (balance, args, nargout,
48  doc: /* -*- texinfo -*-
49 @deftypefn {} {@var{AA} =} balance (@var{A})
50 @deftypefnx {} {@var{AA} =} balance (@var{A}, @var{opt})
51 @deftypefnx {} {[@var{DD}, @var{AA}] =} balance (@var{A}, @var{opt})
52 @deftypefnx {} {[@var{D}, @var{P}, @var{AA}] =} balance (@var{A}, @var{opt})
53 @deftypefnx {} {[@var{CC}, @var{DD}, @var{AA}, @var{BB}] =} balance (@var{A}, @var{B}, @var{opt})
54 
55 Balance the matrix @var{A} to reduce numerical errors in future
56 calculations.
57 
58 Compute @code{@var{AA} = @var{DD} \ @var{A} * @var{DD}} in which @var{AA}
59 is a matrix whose row and column norms are roughly equal in magnitude, and
60 @code{@var{DD} = @var{P} * @var{D}}, in which @var{P} is a permutation
61 matrix and @var{D} is a diagonal matrix of powers of two. This allows the
62 equilibration to be computed without round-off. Results of eigenvalue
63 calculation are typically improved by balancing first.
64 
65 If two output values are requested, @code{balance} returns
66 the diagonal @var{D} and the permutation @var{P} separately as vectors.
67 In this case, @code{@var{DD} = eye(n)(:,@var{P}) * diag (@var{D})}, where
68 @math{n} is the matrix size.
69 
70 If four output values are requested, compute @code{@var{AA} =
71 @var{CC}*@var{A}*@var{DD}} and @code{@var{BB} = @var{CC}*@var{B}*@var{DD}},
72 in which @var{AA} and @var{BB} have nonzero elements of approximately the
73 same magnitude and @var{CC} and @var{DD} are permuted diagonal matrices as
74 in @var{DD} for the algebraic eigenvalue problem.
75 
76 The eigenvalue balancing option @var{opt} may be one of:
77 
78 @table @asis
79 @item @qcode{"noperm"}, @qcode{"S"}
80 Scale only; do not permute.
81 
82 @item @qcode{"noscal"}, @qcode{"P"}
83 Permute only; do not scale.
84 @end table
85 
86 Algebraic eigenvalue balancing uses standard @sc{lapack} routines.
87 
88 Generalized eigenvalue problem balancing uses Ward's algorithm
89 (SIAM Journal on Scientific and Statistical Computing, 1981).
90 @end deftypefn */)
91 {
92  int nargin = args.length ();
93 
94  if (nargin < 1 || nargin > 3 || nargout < 0)
95  print_usage ();
96 
97  octave_value_list retval;
98 
99  // determine if it's AEP or GEP
100  bool AEPcase = nargin == 1 || args(1).is_string ();
101 
102  // problem dimension
103  octave_idx_type nn = args(0).rows ();
104 
105  if (nn != args(0).columns ())
106  err_square_matrix_required ("balance", "A");
107 
108  bool isfloat = args(0).is_single_type ()
109  || (! AEPcase && args(1).is_single_type ());
110 
111  bool complex_case = args(0).iscomplex ()
112  || (! AEPcase && args(1).iscomplex ());
113 
114  // Extract argument 1 parameter for both AEP and GEP.
115  Matrix aa;
116  ComplexMatrix caa;
117  FloatMatrix faa;
118  FloatComplexMatrix fcaa;
119 
120  if (isfloat)
121  {
122  if (complex_case)
123  fcaa = args(0).float_complex_matrix_value ();
124  else
125  faa = args(0).float_matrix_value ();
126  }
127  else
128  {
129  if (complex_case)
130  caa = args(0).complex_matrix_value ();
131  else
132  aa = args(0).matrix_value ();
133  }
134 
135  // Treat AEP/GEP cases.
136  if (AEPcase)
137  {
138  // Algebraic eigenvalue problem.
139  bool noperm = false;
140  bool noscal = false;
141  if (nargin > 1)
142  {
143  std::string a1s = args(1).string_value ();
144  noperm = a1s == "noperm" || a1s == "S";
145  noscal = a1s == "noscal" || a1s == "P";
146  }
147 
148  // balance the AEP
149  if (isfloat)
150  {
151  if (complex_case)
152  {
153  octave::math::aepbalance<FloatComplexMatrix> result (fcaa, noperm, noscal);
154 
155  if (nargout == 0 || nargout == 1)
156  retval = ovl (result.balanced_matrix ());
157  else if (nargout == 2)
158  retval = ovl (result.balancing_matrix (),
159  result.balanced_matrix ());
160  else
161  retval = ovl (result.scaling_vector (),
162  result.permuting_vector (),
163  result.balanced_matrix ());
164  }
165  else
166  {
167  octave::math::aepbalance<FloatMatrix> result (faa, noperm, noscal);
168 
169  if (nargout == 0 || nargout == 1)
170  retval = ovl (result.balanced_matrix ());
171  else if (nargout == 2)
172  retval = ovl (result.balancing_matrix (),
173  result.balanced_matrix ());
174  else
175  retval = ovl (result.scaling_vector (),
176  result.permuting_vector (),
177  result.balanced_matrix ());
178  }
179  }
180  else
181  {
182  if (complex_case)
183  {
184  octave::math::aepbalance<ComplexMatrix> result (caa, noperm, noscal);
185 
186  if (nargout == 0 || nargout == 1)
187  retval = ovl (result.balanced_matrix ());
188  else if (nargout == 2)
189  retval = ovl (result.balancing_matrix (),
190  result.balanced_matrix ());
191  else
192  retval = ovl (result.scaling_vector (),
193  result.permuting_vector (),
194  result.balanced_matrix ());
195  }
196  else
197  {
198  octave::math::aepbalance<Matrix> result (aa, noperm, noscal);
199 
200  if (nargout == 0 || nargout == 1)
201  retval = ovl (result.balanced_matrix ());
202  else if (nargout == 2)
203  retval = ovl (result.balancing_matrix (),
204  result.balanced_matrix ());
205  else
206  retval = ovl (result.scaling_vector (),
207  result.permuting_vector (),
208  result.balanced_matrix ());
209  }
210  }
211  }
212  else
213  {
214  std::string bal_job;
215  if (nargout == 1)
216  warning ("balance: used GEP, should have two output arguments");
217 
218  // Generalized eigenvalue problem.
219  if (nargin == 2)
220  bal_job = 'B';
221  else
222  bal_job = args(2).xstring_value ("balance: OPT argument must be a string");
223 
224  if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
225  err_nonconformant ();
226 
227  Matrix bb;
228  ComplexMatrix cbb;
229  FloatMatrix fbb;
230  FloatComplexMatrix fcbb;
231 
232  if (isfloat)
233  {
234  if (complex_case)
235  fcbb = args(1).float_complex_matrix_value ();
236  else
237  fbb = args(1).float_matrix_value ();
238  }
239  else
240  {
241  if (complex_case)
242  cbb = args(1).complex_matrix_value ();
243  else
244  bb = args(1).matrix_value ();
245  }
246 
247  // balance the GEP
248  if (isfloat)
249  {
250  if (complex_case)
251  {
252  octave::math::gepbalance<FloatComplexMatrix> result (fcaa, fcbb, bal_job);
253 
254  switch (nargout)
255  {
256  case 4:
257  retval(3) = result.balanced_matrix2 ();
258  OCTAVE_FALLTHROUGH;
259 
260  case 3:
261  retval(2) = result.balanced_matrix ();
262  retval(1) = result.balancing_matrix2 ();
263  retval(0) = result.balancing_matrix ();
264  break;
265 
266  case 2:
267  retval(1) = result.balancing_matrix2 ();
268  OCTAVE_FALLTHROUGH;
269 
270  case 1:
271  retval(0) = result.balancing_matrix ();
272  break;
273 
274  default:
275  error ("balance: invalid number of output arguments");
276  break;
277  }
278  }
279  else
280  {
281  octave::math::gepbalance<FloatMatrix> result (faa, fbb, bal_job);
282 
283  switch (nargout)
284  {
285  case 4:
286  retval(3) = result.balanced_matrix2 ();
287  OCTAVE_FALLTHROUGH;
288 
289  case 3:
290  retval(2) = result.balanced_matrix ();
291  retval(1) = result.balancing_matrix2 ();
292  retval(0) = result.balancing_matrix ();
293  break;
294 
295  case 2:
296  retval(1) = result.balancing_matrix2 ();
297  OCTAVE_FALLTHROUGH;
298 
299  case 1:
300  retval(0) = result.balancing_matrix ();
301  break;
302 
303  default:
304  error ("balance: invalid number of output arguments");
305  break;
306  }
307  }
308  }
309  else
310  {
311  if (complex_case)
312  {
313  octave::math::gepbalance<ComplexMatrix> result (caa, cbb, bal_job);
314 
315  switch (nargout)
316  {
317  case 4:
318  retval(3) = result.balanced_matrix2 ();
319  OCTAVE_FALLTHROUGH;
320 
321  case 3:
322  retval(2) = result.balanced_matrix ();
323  retval(1) = result.balancing_matrix2 ();
324  retval(0) = result.balancing_matrix ();
325  break;
326 
327  case 2:
328  retval(1) = result.balancing_matrix2 ();
329  OCTAVE_FALLTHROUGH;
330 
331  case 1:
332  retval(0) = result.balancing_matrix ();
333  break;
334 
335  default:
336  error ("balance: invalid number of output arguments");
337  break;
338  }
339  }
340  else
341  {
342  octave::math::gepbalance<Matrix> result (aa, bb, bal_job);
343 
344  switch (nargout)
345  {
346  case 4:
347  retval(3) = result.balanced_matrix2 ();
348  OCTAVE_FALLTHROUGH;
349 
350  case 3:
351  retval(2) = result.balanced_matrix ();
352  retval(1) = result.balancing_matrix2 ();
353  retval(0) = result.balancing_matrix ();
354  break;
355 
356  case 2:
357  retval(1) = result.balancing_matrix2 ();
358  OCTAVE_FALLTHROUGH;
359 
360  case 1:
361  retval(0) = result.balancing_matrix ();
362  break;
363 
364  default:
365  error ("balance: invalid number of output arguments");
366  break;
367  }
368  }
369  }
370  }
371 
372  return retval;
373 }
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:53
error
void error(const char *fmt,...)
Definition: error.cc:578
nn
static F77_INT nn
Definition: DASPK.cc:62
err_square_matrix_required
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:118
octave_value_list
Definition: ovl.h:39
nargout
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:997
isfloat
OCTAVE_EXPORT octave_value_list isfloat
Definition: data.cc:3212
octave::err_nonconformant
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Definition: lo-array-errwarn.cc:61
retval
octave_value retval
Definition: data.cc:6246
Matrix
Definition: dMatrix.h:36
result
With real return the complex result
Definition: data.cc:3260
warning
void warning(const char *fmt,...)
Definition: error.cc:801
ovl
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).isinteger())
octave_idx_type
nargin
args.length() nargin
Definition: file-io.cc:589
octave_value::is_string
bool is_string(void) const
Definition: ov.h:577
string
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:888