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