Mercurial > hg > octave-lojdl
changeset 1959:9fae6fc592f2
[project @ 1996-02-15 04:54:32 by jwe]
author | jwe |
---|---|
date | Thu, 15 Feb 1996 04:54:32 +0000 |
parents | 9ca852da0017 |
children | 285a7f683a4c |
files | liboctave/dMatrix.cc liboctave/dMatrix.h |
diffstat | 2 files changed, 105 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -96,6 +96,17 @@ double F77_FCN (dlange, DLANGE) (const char*, const int&, const int&, const double*, const int&, double*); + + int F77_FCN (qzhes, QZHES) (const int&, const int&, double*, + double*, const long&, double*); + + int F77_FCN (qzit, QZIT) (const int&, const int&, double*, double*, + const double&, const long&, double*, + int&); + + int F77_FCN (qzval, QZVAL) (const int&, const int&, double*, + double*, double*, double*, double*, + const long&, double*); } // Matrix class. @@ -2684,6 +2695,96 @@ return retval; } +ComplexColumnVector +Qzval (const Matrix& a, const Matrix& b) +{ + ComplexColumnVector retval; + + int a_nr = a.rows(); + int a_nc = a.cols(); + + int b_nr = b.rows(); + int b_nc = b.cols(); + + if (a_nr == a_nc) + { + if (a_nr == b_nr && a_nc == b_nc) + { + if (a_nr != 0) + { + Matrix jnk (a_nr, a_nr, 0.0); + double *pjnk = jnk.fortran_vec (); + + ColumnVector alfr (a_nr); + double *palfr = alfr.fortran_vec (); + + ColumnVector alfi (a_nr); + double *palfi = alfr.fortran_vec (); + + ColumnVector beta (a_nr); + double *pbeta = alfr.fortran_vec (); + + Matrix atmp = a; + double *pa = atmp.fortran_vec (); + + Matrix btmp = b; + double *pb = btmp.fortran_vec (); + + long matz = 0; + int info; + + // XXX FIXME ??? XXX + double eps = DBL_EPSILON; + + F77_FCN (qzhes, QZHES) (a_nr, a_nr, pa, pb, matz, pjnk); + + F77_FCN (qzit, QZIT) (a_nr, a_nr, pa, pb, eps, matz, pjnk, info); + + if (! info) + { + F77_FCN (qzval, QZVAL) (a_nr, a_nr, pa, pb, palfr, + palfi, pbeta, matz, pjnk); + + // Count and extract finite generalized eigenvalues. + + int cnt = 0; + + for (int i = 0; i < a_nr; i++) + if (beta.elem (i) != 0) + cnt++; + + ComplexColumnVector cx (cnt, 0.0); + + Complex Im (0, 1); + + for (int i = 0; i < a_nr; i++) + { + if (beta.elem (i) != 0) + { + // Finite generalized eigenvalue. + + cnt--; + cx.elem (cnt) = (alfr.elem (i) + Im * alfi.elem (i)) + / beta.elem (i); + } + } + + retval = cx; + } + else + (*current_liboctave_error_handler) + ("qzval: trouble in qzit, info = %d", info); + } + } + else + (*current_liboctave_error_handler) ("nonconformant matrices"); + } + else + (*current_liboctave_error_handler) ("qzval: square matrices required"); + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -242,9 +242,11 @@ Matrix (double *d, int r, int c) : MArray2<double> (d, r, c) { } }; -Matrix Givens (double, double); +extern Matrix Givens (double, double); -Matrix Sylvester (const Matrix&, const Matrix&, const Matrix&); +extern Matrix Sylvester (const Matrix&, const Matrix&, const Matrix&); + +extern ComplexColumnVector Qzval (const Matrix& a, const Matrix& b); #endif