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