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