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