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-qzval.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  4KB  |  189 lines
  1. // f-qzval.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 <float.h>
  31.  
  32. #include "dMatrix.h"
  33. #include "dColVector.h"
  34. #include "CColVector.h"
  35. #include "f77-uscore.h"
  36.  
  37. #include "tree-const.h"
  38. #include "user-prefs.h"
  39. #include "gripes.h"
  40. #include "error.h"
  41. #include "utils.h"
  42. #include "help.h"
  43. #include "defun-dld.h"
  44.  
  45. extern "C"
  46. {
  47.   int F77_FCN (qzhes) (const int*, const int*, double*, double*, const
  48.                long*, double*);
  49.  
  50.   int F77_FCN (qzit) (const int*, const int*, double*, double*, const
  51.               double*, const long*, double*, int*);
  52.  
  53.   int F77_FCN (qzval) (const int*, const int*, double*, double*,
  54.                double*, double*, double*, const long*, double*);
  55. }
  56.  
  57. DEFUN_DLD_BUILTIN ("qzval", Fqzval, Sqzval, 3, 1,
  58.   "X = qzval (A, B)\n\
  59. \n\
  60. compute generalized eigenvalues of the matrix pencil (A - lambda B).\n\
  61. A and B must be real matrices.")
  62. {
  63.   Octave_object retval;
  64.  
  65.   int nargin = args.length ();
  66.  
  67.   if (nargin != 2 || nargout > 1)
  68.     {
  69.       print_usage ("qzval");
  70.       return retval;
  71.     }
  72.  
  73.   tree_constant arg_a = args(0);
  74.   tree_constant arg_b = args(1);
  75.  
  76.   int a_nr = arg_a.rows();
  77.   int a_nc = arg_a.columns();
  78.  
  79.   int b_nr = arg_b.rows();
  80.   int b_nc = arg_b.columns();
  81.  
  82.   int arg_a_is_empty = empty_arg ("qzval", a_nr, a_nc);
  83.   int arg_b_is_empty = empty_arg ("qzval", b_nr, b_nc);
  84.  
  85.   if (arg_a_is_empty > 0 && arg_b_is_empty > 0)
  86.     return Matrix ();
  87.   else if (arg_a_is_empty || arg_b_is_empty)
  88.     return retval;
  89.  
  90. // Arguments are not empty, so check for correct dimensions.
  91.  
  92.   if (a_nr != a_nc || b_nr != b_nc)
  93.     {
  94.       gripe_square_matrix_required ("qzval: first two parameters:");
  95.       return retval;
  96.     }
  97.  
  98.   if (a_nr != b_nr)
  99.     {
  100.       gripe_nonconformant ();
  101.       return retval;
  102.     }
  103.   
  104. // Dimensions look o.k., let's solve the problem.
  105.  
  106.   if (arg_a.is_complex_type () || arg_b.is_complex_type ())
  107.     {
  108.       error ("qzval: cannot yet do complex matrix arguments\n");
  109.       return retval;
  110.     }
  111.  
  112. // Do everything in real arithmetic.
  113.  
  114.   Matrix jnk (a_nr, a_nr, 0.0);
  115.  
  116.   ColumnVector alfr (a_nr);
  117.   ColumnVector alfi (a_nr);
  118.   ColumnVector beta (a_nr);
  119.  
  120.   long matz = 0;
  121.   int info;
  122.  
  123. // XXX FIXME ??? XXX
  124.   double eps = DBL_EPSILON;
  125.  
  126.   Matrix ca = arg_a.matrix_value ();
  127.  
  128.   if (error_state)
  129.     return retval;
  130.  
  131.   Matrix cb = arg_b.matrix_value ();
  132.  
  133.   if (error_state)
  134.     return retval;
  135.  
  136. // Use EISPACK qz functions.
  137.  
  138.   F77_FCN (qzhes) (&a_nr, &a_nr, ca.fortran_vec (),
  139.            cb.fortran_vec (), &matz, jnk.fortran_vec ());
  140.  
  141.   F77_FCN (qzit) (&a_nr, &a_nr, ca.fortran_vec (),
  142.           cb.fortran_vec (), &eps, &matz,
  143.           jnk.fortran_vec (), &info);  
  144.  
  145.   if (info)
  146.     error ("qzval: trouble in qzit, info = %d", info);
  147.  
  148.   F77_FCN (qzval) (&a_nr, &a_nr, ca.fortran_vec (),
  149.            cb.fortran_vec (), alfr.fortran_vec (),
  150.            alfi.fortran_vec (), beta.fortran_vec (),
  151.            &matz, jnk.fortran_vec ());
  152.  
  153. // Count and extract finite generalized eigenvalues.
  154.  
  155.   int i;
  156.   int cnt = 0;
  157.  
  158.   Complex Im (0, 1);
  159.  
  160.   for (i = 0; i < a_nr; i++)
  161.     if (beta (i) != 0)
  162.       cnt++;
  163.  
  164.   ComplexColumnVector cx (cnt, 0.0);
  165.  
  166.   for (i = 0; i < a_nr; i++)
  167.     {
  168.       if (beta (i) != 0)
  169.     {
  170.  
  171. // Finite generalized eigenvalue.
  172.  
  173.       cnt--;
  174.       cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i);
  175.     }
  176.     }
  177.  
  178.   retval = cx;
  179.  
  180.   return retval;
  181. }
  182.  
  183. /*
  184. ;;; Local Variables: ***
  185. ;;; mode: C++ ***
  186. ;;; page-delimiter: "^/\\*" ***
  187. ;;; End: ***
  188. */
  189.