Mercurial > hg > octave-terminal
changeset 44:f3215b07c171
[project @ 1993-08-10 22:33:20 by jwe]
Initial revision
author | jwe |
---|---|
date | Tue, 10 Aug 1993 22:33:20 +0000 |
parents | e399beacf758 |
children | 53764ee8f358 |
files | src/qzval.cc |
diffstat | 1 files changed, 165 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/src/qzval.cc @@ -0,0 +1,165 @@ +// tc-qzval.cc -*- C++ -*- +/* + +Copyright (C) 1993 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef __GNUG__ +#pragma implementation +#endif + +#include "Matrix.h" + +#include "tree-const.h" +#include "user-prefs.h" +#include "gripes.h" +#include "error.h" + +int F77_FCN (qzhes) (const int*, const int*, double*, double*, const + long*, double*); + +int F77_FCN (qzit) (const int*, const int*, double*, double*, const + double*, const long*, double*, int*); + +int F77_FCN (qzval) (const int*, const int*, double*, double*, + double*, double*, double*, const long*, double*); + +// XXX FIXME XXX +extern int empty_arg (tree_constant&); +extern tree_constant* empty_tree (int , char*); +extern ComplexMatrix ComplexMatrixLoad (tree_constant&); + +#ifdef WITH_DLD +tree_constant * +builtin_qzvalue_2 (tree_constant *args, int nargin, int nargout) +{ + return qzvalue (args, nargin, nargout); +} +#endif + +tree_constant * +qzvalue (tree_constant *args, int nargin, int nargout) +{ + tree_constant *retval = NULL_TREE_CONST; + + tree_constant arga = args[1].make_numeric (); + tree_constant argb = args[2].make_numeric(); + + if (empty_arg (arga) || empty_arg (argb)) + retval = empty_tree (nargout, "qzvalue"); + else + { + +// Arguments are not empty, so check for correct dimensions. + + int a_rows = arga.rows(); + int a_cols = arga.columns(); + int b_rows = argb.rows(); + int b_cols = argb.columns(); + + if ((a_rows != a_cols) || (b_rows != b_cols)) + { + gripe_square_matrix_required ("qzvalue: first two parameters:"); + return retval; + } + else if (a_rows != b_rows) + { + gripe_nonconformant (); + return retval; + } + +// Dimensions look o.k., let's solve the problem. + + retval = new tree_constant[nargout+1]; + + if (arga.is_complex_type () || argb.is_complex_type ()) + error ("qzvalue: cannot yet do complex matrix arguments\n"); + else + { + +// Do everything in real arithmetic. + + Matrix jnk (a_rows, a_rows, 0.0); + + ColumnVector alfr (a_rows); + ColumnVector alfi (a_rows); + ColumnVector beta (a_rows); + + long matz = 0; + int info; + +// XXX FIXME ??? XXX + double eps = DBL_EPSILON; + + Matrix ca = arga.to_matrix (); + Matrix cb = argb.to_matrix (); + +// Use EISPACK qz functions. + + F77_FCN (qzhes) (&a_rows, &a_rows, ca.fortran_vec (), + cb.fortran_vec (), &matz, jnk.fortran_vec ()); + + F77_FCN (qzit) (&a_rows, &a_rows, ca.fortran_vec (), + cb.fortran_vec (), &eps, &matz, + jnk.fortran_vec (), &info); + + if (info) + error ("qzvalue: trouble in qzit, info = %d", info); + + F77_FCN (qzval) (&a_rows, &a_rows, ca.fortran_vec (), + cb.fortran_vec (), alfr.fortran_vec (), + alfi.fortran_vec (), beta.fortran_vec (), + &matz, jnk.fortran_vec ()); + +// Count and extract finite generalized eigenvalues. + + int i, cnt; + Complex Im (0, 1); + for (i = 0, cnt = 0; i < a_rows; i++) + if (beta (i) != 0) + cnt++; + + ComplexColumnVector cx (cnt, 0.0); + + for (i = 0; i < a_rows; i++) + { + if (beta (i) != 0) + { + +// Finite generalized eigenvalue. + + cnt--; + cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i); + } + } + retval[0] = tree_constant (cx); + } + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/