GNU Octave  8.1.0
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-2023 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)
OCTAVE_END_NAMESPACE
OCTAVE_END_NAMESPACE(octave)
nn
static F77_INT nn
Definition: DASPK.cc:66
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
OCTINTERP_API 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:1054
error
void error(const char *fmt,...)
Definition: error.cc:979
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:211