GNU Octave 7.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-2022 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
47OCTAVE_NAMESPACE_BEGIN
48
49DEFUN (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
57Balance the matrix @var{A} to reduce numerical errors in future
58calculations.
59
60Compute @code{@var{AA} = @var{DD} \ @var{A} * @var{DD}} in which @var{AA}
61is 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
63matrix and @var{D} is a diagonal matrix of powers of two. This allows the
64equilibration to be computed without round-off. Results of eigenvalue
65calculation are typically improved by balancing first.
66
67If two output values are requested, @code{balance} returns
68the diagonal @var{D} and the permutation @var{P} separately as vectors.
69In this case, @code{@var{DD} = eye(n)(:,@var{P}) * diag (@var{D})}, where
70@math{n} is the matrix size.
71
72If 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}},
74in which @var{AA} and @var{BB} have nonzero elements of approximately the
75same magnitude and @var{CC} and @var{DD} are permuted diagonal matrices as
76in @var{DD} for the algebraic eigenvalue problem.
77
78The eigenvalue balancing option @var{opt} may be one of:
79
80@table @asis
81@item @qcode{"noperm"}, @qcode{"S"}
82Scale only; do not permute.
83
84@item @qcode{"noscal"}, @qcode{"P"}
85Permute only; do not scale.
86@end table
87
88Algebraic eigenvalue balancing uses standard @sc{lapack} routines.
89
90Generalized 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
377OCTAVE_NAMESPACE_END
nn
static F77_INT nn
Definition: DASPK.cc:66
Matrix
Definition: dMatrix.h:42
octave_idx_type
octave_value_list
Definition: ovl.h:44
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:1055
error
void error(const char *fmt,...)
Definition: error.cc:980
err_square_matrix_required
void err_square_matrix_required(const char *fcn, const char *name)
Definition: errwarn.cc:122
octave::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