home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / src / f-balance.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  7KB  |  290 lines
  1. // f-balance.cc                                           -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1993, 1994, 1995 John W. Eaton
  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 2, or (at your option) any
  11. 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, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. // Written by A. S. Hodel <scotte@eng.auburn.edu>
  25.  
  26. #ifdef HAVE_CONFIG_H
  27. #include "config.h"
  28. #endif
  29.  
  30. #include "dMatrix.h"
  31. #include "CMatrix.h"
  32. #include "dbleAEPBAL.h"
  33. #include "CmplxAEPBAL.h"
  34. #include "dbleAEPBAL.h"
  35. #include "CmplxAEPBAL.h"
  36. #include "dbleGEPBAL.h"
  37.  
  38. #include "tree-const.h"
  39. #include "user-prefs.h"
  40. #include "gripes.h"
  41. #include "error.h"
  42. #include "utils.h"
  43. #include "help.h"
  44. #include "defun-dld.h"
  45.  
  46. DEFUN_DLD_BUILTIN ("balance", Fbalance, Sbalance, 4, 4,
  47.   "AA = balance (A [, OPT]) or [[DD,] AA] =  balance (A [, OPT])\n\
  48. \n\
  49. generalized eigenvalue problem:\n\
  50. \n\
  51.   [cc, dd, aa, bb] = balance (a, b [, opt])\n\
  52. \n\
  53. where OPT is an optional single character argument as follows: \n\
  54. \n\
  55.   N: no balancing; arguments copied, transformation(s) set to identity\n\
  56.   P: permute argument(s) to isolate eigenvalues where possible\n\
  57.   S: scale to improve accuracy of computed eigenvalues\n\
  58.   B: (default) permute and scale, in that order.  Rows/columns\n\
  59.      of a (and b) that are isolated by permutation are not scaled\n\
  60. \n\
  61. [DD, AA] = balance (A, OPT) returns aa = dd\a*dd,\n\
  62. \n\
  63. [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)")
  64. {
  65.   Octave_object retval;
  66.  
  67.   int nargin = args.length ();
  68.  
  69.   if (nargin < 1 || nargin > 3 || nargout < 0 || nargout > 4)
  70.     {
  71.       print_usage ("balance");
  72.       return retval;
  73.     }
  74.  
  75.   char *bal_job;
  76.   int my_nargin;        // # args w/o optional string arg
  77.  
  78. // Determine if balancing option is listed.  Set my_nargin to the
  79. // number of matrix inputs.
  80.  
  81.   if (args(nargin-1).is_string ())
  82.     {
  83.       bal_job = args(nargin-1).string_value ();
  84.       my_nargin = nargin-1;
  85.     }
  86.   else
  87.     {
  88.       bal_job = "B";
  89.       my_nargin = nargin;
  90.     }
  91.  
  92.   tree_constant arg_a = args(0);
  93.  
  94.   int a_nr = arg_a.rows ();
  95.   int a_nc = arg_a.columns ();
  96.  
  97. // Check argument 1 dimensions.
  98.  
  99.   int arg_is_empty = empty_arg ("balance", a_nr, a_nc);
  100.  
  101.   if (arg_is_empty < 0)
  102.     return retval;
  103.   if (arg_is_empty > 0)
  104.     return Octave_object (2, Matrix ());
  105.  
  106.   if (a_nr != a_nc)
  107.     {
  108.       gripe_square_matrix_required ("balance");
  109.       return retval;
  110.     }
  111.  
  112. // Extract argument 1 parameter for both AEP and GEP.
  113.  
  114.   Matrix aa;
  115.   ComplexMatrix caa;
  116.   if (arg_a.is_complex_type ())
  117.     caa = arg_a.complex_matrix_value ();
  118.   else
  119.     aa = arg_a.matrix_value ();
  120.  
  121.   if (error_state)
  122.     return retval;
  123.  
  124. // Treat AEP/ GEP cases.
  125.  
  126.   switch (my_nargin)
  127.     {
  128.     case 1:
  129.       
  130. // Algebraic eigenvalue problem.
  131.  
  132.       if (arg_a.is_complex_type ())
  133.     {
  134.       ComplexAEPBALANCE result (caa, bal_job);
  135.  
  136.       if (nargout == 0 || nargout == 1)
  137.         retval(0) = result.balanced_matrix ();
  138.       else
  139.         {
  140.           retval(1) = result.balanced_matrix ();
  141.           retval(0) = result.balancing_matrix ();
  142.         }
  143.     }
  144.       else
  145.     {
  146.       AEPBALANCE result (aa, bal_job);
  147.  
  148.       if (nargout == 0 || nargout == 1)
  149.         retval(0) = result.balanced_matrix ();
  150.       else
  151.         {
  152.           retval(1) = result.balanced_matrix ();
  153.           retval(0) = result.balancing_matrix ();
  154.         }
  155.     }
  156.       break;
  157.  
  158.     case 2:
  159.       {
  160. // Generalized eigenvalue problem.
  161.  
  162. // 1st we have to check argument 2 dimensions and type...
  163.  
  164.     tree_constant arg_b = args(1);
  165.  
  166.     int b_nr = arg_b.rows ();
  167.     int b_nc = arg_b.columns ();
  168.       
  169. // Check argument 2 dimensions -- must match arg 1.
  170.  
  171.     if (b_nr != b_nc || b_nr != a_nr)
  172.       {
  173.         gripe_nonconformant ();
  174.         return retval;
  175.       }
  176.       
  177. // Now, extract the second matrix...
  178. // Extract argument 1 parameter for both AEP and GEP.
  179.  
  180.     Matrix bb;
  181.     ComplexMatrix cbb;
  182.     if (arg_b.is_complex_type ())
  183.       cbb = arg_b.complex_matrix_value ();
  184.     else
  185.       bb = arg_b.matrix_value ();
  186.  
  187.     if (error_state)
  188.       return retval;
  189.  
  190. // Both matrices loaded, now let's check what kind of arithmetic:
  191.  
  192.     if (arg_a.is_complex_type () || arg_b.is_complex_type ())
  193.       {
  194.         if (arg_a.is_real_type ())
  195.           caa = aa;
  196.  
  197.         if (arg_b.is_real_type ())
  198.           cbb = bb;
  199.  
  200. // Compute magnitudes of elements for balancing purposes.
  201. // Surely there's a function I can call someplace!
  202.  
  203.         for (int i = 0; i < a_nr; i++)
  204.           for (int j = 0; j < a_nc; j++)
  205.         {
  206.           aa.elem (i, j) = abs (caa.elem (i, j));
  207.           bb.elem (i, j) = abs (cbb.elem (i, j));
  208.         }
  209.       }
  210.  
  211.     GEPBALANCE result (aa, bb, bal_job);
  212.  
  213.     if (arg_a.is_complex_type () || arg_b.is_complex_type ())
  214.       {
  215.         caa = result.left_balancing_matrix () * caa
  216.           * result.right_balancing_matrix ();
  217.  
  218.         cbb = result.left_balancing_matrix () * cbb
  219.           * result.right_balancing_matrix ();
  220.  
  221.         switch (nargout)
  222.           {
  223.           case 0:
  224.           case 1:
  225.         warning ("balance: should use two output arguments");
  226.         retval(0) = caa;
  227.         break;
  228.  
  229.           case 2:
  230.         retval(1) = cbb;
  231.         retval(0) = caa;
  232.         break;
  233.  
  234.           case 4:
  235.         retval(3) = cbb;
  236.         retval(2) = caa;
  237.         retval(1) = result.right_balancing_matrix ();
  238.         retval(0) = result.left_balancing_matrix ();
  239.         break;
  240.  
  241.           default:
  242.         error ("balance: invalid number of output arguments");
  243.         break;
  244.           }
  245.       }
  246.     else
  247.       {
  248.         switch (nargout)
  249.           {
  250.           case 0:
  251.           case 1:
  252.         warning ("balance: should use two output arguments");
  253.         retval(0) = result.balanced_a_matrix ();
  254.         break;
  255.  
  256.           case 2:
  257.         retval(1) = result.balanced_b_matrix ();
  258.         retval(0) = result.balanced_a_matrix ();
  259.         break;
  260.  
  261.           case 4:
  262.         retval(3) = result.balanced_b_matrix ();
  263.         retval(2) = result.balanced_a_matrix ();
  264.         retval(1) = result.right_balancing_matrix ();
  265.         retval(0) = result.left_balancing_matrix ();
  266.         break;
  267.  
  268.           default:
  269.         error ("balance: invalid number of output arguments");
  270.         break;
  271.           }
  272.       }
  273.       }
  274.       break;
  275.  
  276.     default:
  277.       error ("balance requires one (AEP) or two (GEP) numeric arguments");
  278.       break;
  279.     }
  280.  
  281.   return retval;
  282. }
  283.  
  284. /*
  285. ;;; Local Variables: ***
  286. ;;; mode: C++ ***
  287. ;;; page-delimiter: "^/\\*" ***
  288. ;;; End: ***
  289. */
  290.