Mercurial > hg > octave-max
changeset 3183:5edc539c2f80
[project @ 1998-09-25 03:37:22 by jwe]
author | jwe |
---|---|
date | Fri, 25 Sep 1998 03:37:22 +0000 |
parents | 936b9ae8f7d2 |
children | 3988763ec9d3 |
files | src/DLD-FUNCTIONS/qz.cc |
diffstat | 1 files changed, 773 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/qz.cc @@ -0,0 +1,773 @@ +/* + +Copyright (C) 1998 A. S. Hodel + +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, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Generalized eigenvalue balancing via LAPACK +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#undef DEBUG +#undef DEBUG_SORT +#undef DEBUG_EIG + +#include "config.h" + +#include <cfloat> +#include <iostream.h> +#include <iomanip.h> +#include <math.h> +#include <string.h> +#include <strstream.h> + +#include "CmplxQRP.h" +#include "dbleQR.h" +#include "defun-dld.h" +#include "error.h" +#include "f77-fcn.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "oct-map.h" +#include "ov.h" +#include "pager.h" +#if defined(DEBUG) || defined(DEBUG_SORT) +#include "pr-output.h" +#endif +#include "symtab.h" +#include "utils.h" +#include "variables.h" + +typedef int (*sort_function) (const int& LSIZE, const double& ALPHA, + const double& BETA, const double& S, const double& P); + +extern "C" +{ + int F77_FCN( dggbal, DGGBAL) (const char* JOB, const int& N, + double* A, const int& LDA, double* B, const int& LDB, + int& ILO, int & IHI, double* LSCALE, + double* RSCALE, double* WORK, int& INFO, long ); + + int F77_FCN( dggbak, DGGBAK) (const char* JOB, const char* SIDE, + const int& N, const int& ILO, const int& IHI, + double* LSCALE, double* RSCALE, int& M, + double* V, const int& LDV, int& INFO, long, long); + + int F77_FCN( dgghrd, DGGHRD) ( const char* COMPQ, const char* COMPZ, + const int& N, const int& ILO, const int& IHI, double* A, const int& LDA, + double* B, const int& LDB, double* Q, const int& LDQ, double* Z, + const int& LDZ, int& INFO, const long, const long); + + int F77_FCN( dhgeqz, DHGEQZ) ( const char* JOB, const char* COMPQ, + const char* COMPZ, const int& N, const int& ILO, const int& IHI, + double* A, const int& LDA, double* B, const int& LDB, + double* ALPHAR, double* ALPHAI, double* BETA, double* Q, + const int& LDQ, double* Z, const int& LDZ, double* WORK, + const int& LWORK, int& INFO, const long, const long, const long ); + + int F77_FCN( dlag2, DLAG2) ( double* A, const int& LDA, double* B, + const int& LDB, const double& SAFMIN, double& SCALE1, double& SCALE2, + double& WR1, double& WR2, double& WI ); + + // Van Dooren's code (netlib.org: toms/590) for reordering + // GEP. Only processes Z, not Q. + int F77_FCN( dsubsp, DSUBSP) ( const int& NMAX, const int& N, double* A, + double* B, double* Z, sort_function, const double& EPS, + int& NDIM, int& FAIL, int* IND); + + // documentation for DTGEVC incorrectly states that VR, VL are + // complex*16; they are declared in DTGEVC as double precision + // (probably a cut and paste problem fro ZTGEVC) + int F77_FCN( dtgevc, DTGEVC) ( const char* SIDE, const char* HOWMNY, + int* SELECT, const int& N, double* A, const int& LDA, double* B, + const int& LDB, double* VL, const int& LDVL, double* VR, + const int& LDVR, const int& MM, int& M, double* WORK, int& INFO, + long, long ); + + int F77_FCN ( xdlamch, XDLAMCH) (const char* cmach, double& retval, long); + int F77_FCN ( xdlange, XDLANGE) (const char*, const int&, + const int&, const double*, + const int&, double*, double&); +} + +// fcrhp, fin, fout, folhp: +// routines for ordering of generalized eigenvalues +// return 1 if test is passed, 0 otherwise +// fin: |lambda| < 1 +// fout: |lambda| >= 1 +// fcrhp: real(lambda) >= 0 +// folhp: real(lambda) < 0 + +static int fcrhp(const int& lsize, const double& alpha, + const double& beta, const double& s, const double& p) +{ + if(lsize == 1) + return (alpha*beta >= 0 ? 1 : -1); + else + return (s >= 0 ? 1 : -1); +} +static int fin(const int& lsize, const double& alpha, + const double& beta, const double& s, const double& p) +{ + int retval; + if(lsize == 1) + retval = (fabs(alpha) < fabs(beta) ? 1 : -1); + else + retval = (fabs(p) < 1 ? 1 : -1); + + #ifdef DEBUG + cout << "qz: fin: retval=" << retval << endl; + #endif + return retval; +} +static int folhp(const int& lsize, const double& alpha, + const double& beta, const double& s, const double& p) +{ + if(lsize == 1) + return (alpha*beta < 0 ? 1 : -1); + else + return (s < 0 ? 1 : -1); +} +static int fout(const int& lsize, const double& alpha, + const double& beta, const double& s, const double& p) +{ + if(lsize == 1) + return (fabs(alpha) >= fabs(beta) ? 1 : -1); + else + return (fabs(p) >= 1 ? 1 : -1); +} + +DEFUN_DLD (qz, args, nargout, +"Usage: lambda = qz(A,B) form [1]\n\ + [AA,BB,Q,Z{,V,W,lambda}] = qz(A,B) form [2]\n\ + [AA,BB,Z{,lambda}] = qz(A,B,opt) form [3]\n\ +Generalized eigenvalue problem A x = s B x \n\ + +Form [1]: Computes the generalized eigenvalues lambda of (A - sB).\n\ +Form [2]: Computes qz decomposition, generalized eigenvectors, and \n\ + generalized eigenvalues of (A - sB)\n\ + A V = B V diag(lambda)\n\ + W' A = diag(lambda) W' B\n\ + AA = Q'*A*Z, BB = Q'*B*Z with Q, Z orthogonal (unitary)= I\n\ +Form [3]: As in form [2], but allows ordering of generalized eigenpairs\n\ + for (e.g.) solution of discrete time algebraic Riccati equations.\n\ + Form 3 is not available for complex matrices and does not compute\n\ + the generalized eigenvectors V, W, nor the orthogonal matrix Q.\n\ +\n\ + opt: for ordering eigenvalues of the GEP pencil. The leading block\n\ + of the revised pencil contains all eigenvalues that satisfy:\n\ + \"N\" = unordered (default) \n\ + \"S\" = small: leading block has all |lambda| <=1 \n\ + \"B\" = big: leading block has all |lambda >= 1 \n\ + \"-\" = negative real part: leading block has all eigenvalues\n\ + in the open left half-plant\n\ + \"+\" = nonnegative real part: leading block has all eigenvalues\n\ + in the closed right half-plane\n\ +\n\ +Note: Permutation balancing is performed, but not scaling (see balance)\n\ + Order of output arguments was selected for compatibility with MATLAB\n\ +\n\ +See also: balance, dare, eig, schur\n") +{ + octave_value_list retval; + int nargin = args.length (); + + #ifdef DEBUG + cout << "qz: nargin = " << nargin << ", nargout = " << nargout << endl; + #endif + + if (nargin < 2 || nargin > 3 || nargout > 7 ) + { + print_usage ("qz"); + return retval; + } + else if(nargin == 3 && (nargout < 3 || nargout > 4)) + { + error("qz: Illegal number of output arguments for form [3] call"); + } + + #ifdef DEBUG + cout << "qz: determine ordering option" << endl; + #endif + + // Determine ordering option + string ord_job; + static double safmin; + if(nargin == 2) + ord_job = "N"; + else if( !args(2).is_string() ) + error("qz: argument 3 must be a string"); + else + { + ord_job = args(2).string_value(); + if(ord_job[0] != 'N' && ord_job[0] != 'S' && ord_job[0] != 'B' + && ord_job[0] != '+' && ord_job[0] != '-') + error("qz: illegal order option"); + + // overflow constant required by dlag2 + F77_XFCN ( xdlamch, XDLAMCH, ("S", safmin, 1L)); + + #ifdef DEBUG_EIG + cout << "qz: initial value of safmin=" << setiosflags(ios::scientific) + << safmin << endl; + #endif + + // some machines (e.g., DEC alpha) get safmin = 0; + // for these, use eps instead to avoid problems in dlag2 + if(safmin == 0) + { + #ifdef DEBUG_EIG + cout << "qz: DANGER WILL ROBINSON: safmin is 0!" << endl; + #endif + + F77_XFCN ( xdlamch, XDLAMCH, ("E", safmin, 1L)); + + #ifdef DEBUG_EIG + cout << "qz: safmin set to " << setiosflags(ios::scientific) + << safmin << endl; + #endif + } + } + + #ifdef DEBUG + cout << "qz: check argument 1" << endl; + #endif + + // Argument 1: check if it's o.k. dimensioned + int nn = args(0).rows(); + + #ifdef DEBUG + cout << "argument 1 dimensions: (" << nn << "," << args(0).columns() << ")" + << endl; + #endif + int arg_is_empty = empty_arg ("qz", nn, args(0).columns()); + if (arg_is_empty < 0) + { + gripe_empty_arg("qz: parameter 1",0); + return retval; + } + else if (arg_is_empty > 0) + { + gripe_empty_arg("qz: parameter 1; continuing",0); + return octave_value_list (2, Matrix ()); + } + else if (args(0).columns() != nn) + { + gripe_square_matrix_required ("qz"); + return retval; + } + + // Argument 1: dimensions look good; get the value + Matrix aa; + ComplexMatrix caa; + if (args(0).is_complex_type ()) + caa = args(0).complex_matrix_value (); + else + aa = args(0).matrix_value (); + if (error_state) + return retval; + + #ifdef DEBUG + cout << "qz: check argument 2" << endl; + #endif + + // Extract argument 2 (bb, or cbb if complex) + if( (nn != args(1).columns()) || (nn != args(1).rows() )) + { + gripe_nonconformant (); + return retval; + } + Matrix bb; + ComplexMatrix cbb; + if (args(1).is_complex_type ()) + cbb = args(1).complex_matrix_value (); + else + bb = args(1).matrix_value (); + if (error_state) + return retval; + + // Both matrices loaded, now let's check what kind of arithmetic: + //declared static to avoid compiler warnings about long jumps, vforks. + static int complex_case + = (args(0).is_complex_type() || args(1).is_complex_type()); + + if(nargin == 3 && complex_case) + error("qz: cannot re-order complex qz decomposition."); + + // first, declare variables used in both the real and complex case + Matrix QQ(nn,nn), ZZ(nn,nn), VR(nn,nn), VL(nn,nn); + RowVector alphar(nn), alphai(nn), betar(nn); + + ComplexMatrix CQ(nn,nn), CZ(nn,nn),CVR(nn,nn),CVL(nn,nn); + int ilo, ihi, info; + char compq = (nargout >= 3 ? 'V' : 'N'), + compz = (nargout >= 4 ? 'V' : 'N'); + + // initialize Q,Z to identity if we need either of them + if(compq == 'V' || compz == 'V') + for(int ii=0; ii < nn ; ii++) + for( int jj=0; jj < nn ; jj++) + QQ(ii,jj) = ZZ(ii,jj) = (ii == jj ? 1.0 : 0.0); + + // always perform permutation balancing + char bal_job = 'P'; + RowVector lscale(nn), rscale(nn), work(6*nn); + + if(complex_case) + error("Complex case not implemented yet"); + else + { + #ifdef DEBUG + if(compq == 'V') + cout << "qz: performing balancing; QQ=" << endl << QQ << endl; + #endif + + F77_XFCN( dggbal, DGGBAL, (&bal_job, nn, aa.fortran_vec(), + nn, bb.fortran_vec() , nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), work.fortran_vec(), info , 1L)); + if(f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in qz(bal)"); + } + + // Since we just want the balancing matrices, we can use dggbal + // for both the real and complex cases; + // left first + if(compq == 'V') + { + F77_XFCN( dggbak, DGGBAK, (&bal_job, "L", + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, QQ.fortran_vec(), + nn, info, 1L, 1L)); + + #ifdef DEBUG + if(compq == 'V') cout << "qz: balancing done; QQ=" << endl << QQ << endl; + #endif + + if(f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-L)"); + } + + // then right + if(compz == 'V') + { + F77_XFCN(dggbak, DGGBAK, (&bal_job, "R", + nn, ilo, ihi, lscale.fortran_vec(), + rscale.fortran_vec(), nn, ZZ.fortran_vec(), + nn, info, 1L, 1L)); + + #ifdef DEBUG + if(compz == 'V') cout << "qz: balancing done; ZZ=" << endl << ZZ << endl; + #endif + + if(f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in qz(bal-R)"); + } + + static char qz_job; + qz_job = (nargout < 2 ? 'E' : 'S'); + if (complex_case) + { + // complex case + if (args(0).is_real_type ()) caa = aa; + if (args(1).is_real_type ()) cbb = bb; + if(compq == 'V') CQ = QQ; + if(compz == 'V') CZ = ZZ; + error("complex case not done yet"); + } + else // real matrices case + { + #ifdef DEBUG + cout << "qz: peforming qr decomposition of bb" << endl; + #endif + + // compute the QR factorization of bb + QR bqr(bb); + + #ifdef DEBUG + cout << "qz: qr(bb) done; now peforming qz decomposition" << endl; + #endif + + bb = bqr.R(); + #ifdef DEBUG + cout << "qz: extracted bb" << endl; + #endif + + aa = (bqr.Q()).transpose()*aa; + #ifdef DEBUG + cout << "qz: updated aa " << endl; + cout << "bqr.Q () = " << endl << bqr.Q () << endl; + if(compq == 'V') cout << "QQ =" << QQ << endl; + #endif + + if(compq == 'V') QQ = QQ*bqr.Q(); + + #ifdef DEBUG + cout << "qz: precursors done..." << endl; + #endif + + + #ifdef DEBUG + cout << "qz: compq = " << compq << ", compz = " << compz << endl; + #endif + + // reduce to generalized hessenberg form + F77_XFCN( dgghrd, DGGHRD, (&compq, &compz, nn, ilo, ihi, aa.fortran_vec(), + nn, bb.fortran_vec(), nn, QQ.fortran_vec(), nn, ZZ.fortran_vec(), + nn, info,1L,1L)); + if(f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in qz(dgghrd)"); + + // check if just computing generalized eigenvalues or if we're + // actually computing the decomposition + + // reduce to generalized Schur form + F77_XFCN( dhgeqz, DHGEQZ, ( &qz_job, &compq, &compz, nn, ilo, ihi, + aa.fortran_vec(), nn, bb.fortran_vec(), nn, alphar.fortran_vec(), + alphai.fortran_vec(), betar.fortran_vec(), QQ.fortran_vec(), + nn, ZZ.fortran_vec(), nn, work.fortran_vec(), nn, info, 1L, 1L, 1L)); + if(f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in qz(dhgeqz)"); + + } + + // order the QZ decomposition? + if(ord_job[0] != 'N') + { + if(complex_case) // probably not needed, but better be safe + error("qz: cannot re-order complex qz decomposition."); + + else + { + #ifdef DEBUG_SORT + cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl; + #endif + + // declared static to avoid vfork/long jump compiler complaints + static sort_function sort_test; + sort_test = NULL; + + switch(ord_job[0]) + { + case 'S': + sort_test = &fin; + break; + case 'B': + sort_test = &fout; + break; + case '+': + sort_test = &fcrhp; + break; + case '-': + sort_test = &folhp; + break; + default: // this should never happen + error("qz: illegal order option"); + } + + int ndim, fail, ind[nn]; + double inf_norm; + F77_XFCN (xdlange, XDLANGE, ("I", nn, nn, aa.fortran_vec (), nn, + work.fortran_vec (), inf_norm)); + + double eps = DBL_EPSILON*inf_norm*nn; + + #ifdef DEBUG_SORT + cout << "qz: calling dsubsp: aa=" << endl; + octave_print_internal(cout,aa,0); + cout << endl << "bb=" << endl; + octave_print_internal(cout,bb,0); + if(compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal(cout,ZZ,0); + } + cout << endl; + cout << "alphar = " << endl; + octave_print_internal(cout,(Matrix) alphar,0); + cout << endl << "alphai = " << endl; + octave_print_internal(cout,(Matrix) alphai,0); + cout << endl << "beta = " << endl; + octave_print_internal(cout,(Matrix) betar,0); + cout << endl; + #endif + + F77_XFCN( dsubsp, DSUBSP, (nn,nn,aa.fortran_vec(), bb.fortran_vec(), + ZZ.fortran_vec(), sort_test, eps, ndim, fail, ind)); + + #ifdef DEBUG + cout << "qz: back from dsubsp: aa=" << endl; + octave_print_internal(cout,aa,0); + cout << endl << "bb=" << endl; + octave_print_internal(cout,bb,0); + if(compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal(cout,ZZ,0); + } + cout << endl; + #endif + + // manually update alphar, alphai, betar + static int jj; + jj=0; + while(jj < nn) + { + #ifdef DEBUG_EIG + cout << "computing gen eig #" << jj << endl; + #endif + + static int zcnt; // number of zeros in this block + if(jj == (nn-1)) + zcnt = 1; + else if(aa(jj+1,jj) == 0) + zcnt = 1; + else zcnt = 2; + + if(zcnt == 1) // real zero + { + #ifdef DEBUG_EIG + cout << " single gen eig:" << endl; + cout << " alphar(" << jj << ") = " << aa(jj,jj) << endl; + cout << " betar( " << jj << ") = " << bb(jj,jj) << endl; + cout << " alphai(" << jj << ") = 0" << endl; + #endif + + alphar(jj) = aa(jj,jj); + alphai(jj) = 0; + betar(jj) = bb(jj,jj); + } + else // complex conjugate pair + { + #ifdef DEBUG_EIG + cout << "qz: calling dlag2:" << endl; + cout << "safmin=" << setiosflags(ios::scientific) << safmin << endl; + for(int idr = jj ; idr <= jj+1 ; idr++) + { + for(int idc = jj ; idc <= jj+1 ; idc++) + { + cout << "aa(" << idr << "," << idc << ")=" + << aa(idr,idc) << endl; + cout << "bb(" << idr << "," << idc << ")=" + << bb(idr,idc) << endl; + } + } + #endif + double scale1, scale2, wr1, wr2, wi; + F77_XFCN( dlag2, DLAG2, ( &aa(jj,jj), nn, &bb(jj,jj), nn, safmin, + scale1, scale2, wr1, wr2, wi)); + + #ifdef DEBUG_EIG + cout << "dlag2 returns: scale1=" << scale1 + << "\tscale2=" << scale2 << endl + << "\twr1=" << wr1 << "\twr2=" << wr2 + << "\twi=" << wi << endl; + #endif + // just to be safe, check if it's a real pair + if(wi == 0) + { + alphar(jj) = wr1; + alphai(jj) = 0; + betar(jj) = scale1; + alphar(jj+1) = wr2; + alphai(jj+1) = 0; + betar(jj+1) = scale2; + } + else + { + alphar(jj) = alphar(jj+1)=wr1; + alphai(jj) = -(alphai(jj+1) = wi); + betar(jj) = betar(jj+1) = scale1; + } + } + + jj += zcnt; // advance past this block + } + + #ifdef DEBUG_SORT + cout << "qz: back from dsubsp: aa=" << endl; + octave_print_internal(cout,aa,0); + cout << endl << "bb=" << endl; + octave_print_internal(cout,bb,0); + if(compz == 'V') + { + cout << endl << "ZZ=" << endl; + octave_print_internal(cout,ZZ,0); + } + cout << endl << "qz: ndim=" << ndim << endl << "fail=" << fail << endl; + cout << "alphar = " << endl; + octave_print_internal(cout,(Matrix) alphar,0); + cout << endl << "alphai = " << endl; + octave_print_internal(cout,(Matrix) alphai,0); + cout << endl << "beta = " << endl; + octave_print_internal(cout,(Matrix) betar,0); + cout << endl; + #endif + } + } + + // compute generalized eigenvalues? + ComplexColumnVector gev; + if(nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4)) + { + if(complex_case) + error("complex case not yet implemented"); + else + { + #ifdef DEBUG + cout << "qz: computing generalized eigenvalues" << endl; + #endif + + // return finite generalized eigenvalues + int ii, cnt = 0; + for( ii=0 ; ii < nn ; ii++) + if(betar(ii) != 0) + cnt++; + ComplexColumnVector tmp(cnt); + for( ii=0 ; ii < nn ; ii++) + if(betar(ii) != 0) + tmp(ii) = Complex(alphar(ii), alphai(ii))/betar(ii); + gev = tmp; + } + } + + // right, left eigenvector matrices + if(nargout >= 5) + { + char side = (nargout == 5 ? 'R' : 'B'), // which side to compute? + howmny = 'B'; // compute all of them and backtransform + int *select = NULL; // dummy pointer; select is not used. + int m; + + if(complex_case) + error("complex type not yet implemented"); + else + { + #ifdef DEBUG + cout << "qz: computing generalized eigenvectors" << endl; + #endif + + VL = QQ; + VR = ZZ; + + F77_XFCN( dtgevc, DTGEVC, ( &side, &howmny, select, nn, aa.fortran_vec(), + nn, bb.fortran_vec(), nn, VL.fortran_vec(), nn, VR.fortran_vec(), + nn, nn, m, work.fortran_vec(), info, 1L, 1L )); + if(f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in qz(dtgevc)"); + + // now construct the complex form of VV, WW + int jj = 0; + while(jj < nn) + { + // see if real or complex eigenvalue + int cinc = 2; // column increment; assume complex eigenvalue + if(jj == (nn-1)) + cinc = 1; // single column + else if(aa(jj+1,jj) == 0) + cinc = 1; + + // now copy the eigenvector(s) to CVR, CVL + if(cinc == 1) + { + int ii; + for(ii = 0; ii < nn ; ii++) + CVR(ii,jj) = VR(ii,jj); + if(side == 'B') + for(ii = 0; ii < nn ; ii++) + CVL(ii,jj) = VL(ii,jj); + } + else // double column; complex vector + { + int ii; + for(ii = 0; ii < nn ; ii++) + { + CVR(ii,jj) = Complex(VR(ii,jj),VR(ii,jj+1)); + CVR(ii,jj+1) = Complex(VR(ii,jj),-VR(ii,jj+1)); + } + if(side == 'B') + for(ii = 0; ii < nn ; ii++) + { + CVL(ii,jj) = Complex(VL(ii,jj),VL(ii,jj+1)); + CVL(ii,jj+1) = Complex(VL(ii,jj),-VL(ii,jj+1)); + } + } + jj += cinc; // advance to next eigenvectors (if any) + } + } + } + + switch(nargout) + { + case 7: + retval(6) = gev; + case 6: // return eigenvectors + retval(5) = CVL; + case 5: // return eigenvectors + retval(4) = CVR; + case 4: + if(nargin == 3) + { + #ifdef DEBUG + cout << "qz: sort: retval(3) = gev = " << endl; + octave_print_internal(cout,gev); + cout << endl; + #endif + retval(3) = gev; + } + else retval(3) = ZZ; + case 3: + if(nargin == 3) + retval(2) = ZZ; + else + retval(2) = QQ; + case 2: + #ifdef DEBUG + cout << "qz: retval(1) = bb = " << endl; + octave_print_internal(cout,bb,0); + cout << endl << "qz: retval(0) = aa = " <<endl; + octave_print_internal(cout,aa,0); + cout << endl; + #endif + retval(1) = bb; + retval(0) = aa; + break; + case 1: + case 0: + #ifdef DEBUG + cout << "qz: retval(0) = gev = " << gev << endl; + #endif + retval(0) = gev; + break; + default: + error("qz: too many return arguments. Sorry. "); + } + + #ifdef DEBUG + cout << "qz: exiting (at long last)" << endl; + #endif + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/