--- a/src/DLD-FUNCTIONS/qz.cc	Fri Sep 25 04:35:09 1998 +0000
+++ b/src/DLD-FUNCTIONS/qz.cc	Sat Sep 26 02:45:59 1998 +0000
@@ -47,61 +47,77 @@
 #include "oct-map.h"
 #include "ov.h"
 #include "pager.h"
-#if defined(DEBUG) || defined(DEBUG_SORT)
+#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);
+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 (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 (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 (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, long, 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 (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,
+				long, long, 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 );
+  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);
+  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 (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&,
+  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&);
 }
@@ -114,64 +130,79 @@
 //    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)
+static int
+fcrhp (const int& lsize, const double& alpha,
+       const double& beta, const double& s, const double&)
 {
-  if(lsize == 1)
+  if (lsize == 1)
     return (alpha*beta >= 0 ? 1 : -1);
-  else 
+  else
     return (s >= 0 ? 1 : -1);
 }
-static int fin(const int& lsize, const double& alpha, 
-	const double& beta, const double& s, const double& p)
+
+static int
+fin (const int& lsize, const double& alpha,
+     const double& beta, const double&, const double& p)
 {
   int retval;
-  if(lsize == 1)
-    retval = (fabs(alpha) < fabs(beta) ? 1 : -1);
-  else 
-    retval = (fabs(p) < 1 ? 1 : -1);
 
-  #ifdef DEBUG
+  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
+#endif
+
   return retval;
 }
-static int folhp(const int& lsize, const double& alpha, 
-	const double& beta, const double& s, const double& p)
+
+static int
+folhp (const int& lsize, const double& alpha,
+       const double& beta, const double& s, const double&)
 {
-  if(lsize == 1)
+  if (lsize == 1)
     return (alpha*beta < 0 ? 1 : -1);
-  else 
+  else
     return (s < 0 ? 1 : -1);
 }
-static int fout(const int& lsize, const double& alpha, 
-	const double& beta, const double& s, const double& p)
+
+static int
+fout (const int& lsize, const double& alpha,
+      const double& beta, const double&, const double& p)
 {
-  if(lsize == 1)
-    return (fabs(alpha) >= fabs(beta) ? 1 : -1);
-  else 
-    return (fabs(p) >= 1 ? 1 : -1);
+  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\
+"Usage:\n\
+
+  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\
+\n\
 Generalized eigenvalue problem    A x = s B x \n\
-
+\n\
 Form [1]: Computes the generalized eigenvalues lambda of (A - sB).\n\
+\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\
+          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\
+\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\
+          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\
+     opt: for ordering eigenvalues of the GEP pencil.  The leading  block\n\
+          of the revised pencil contains all eigenvalues that satisfy:\n\
+\n\
           \"N\" = unordered (default) \n\
           \"S\" = small: leading block has all |lambda| <=1 \n\
           \"B\" = big: leading block has all |lambda >= 1 \n\
@@ -183,586 +214,705 @@
 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")
+See also: balance, dare, eig, schur")
 {
   octave_value_list retval;
   int nargin = args.length ();
 
-  #ifdef DEBUG
+#ifdef DEBUG
   cout << "qz: nargin = " << nargin << ", nargout = " << nargout << endl;
-  #endif
+#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");
-  }
+  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");
+      return retval;
+    }
 
-  #ifdef DEBUG
+#ifdef DEBUG
   cout << "qz: determine ordering option" << endl;
-  #endif
+#endif
 
   // Determine ordering option
   string ord_job;
   static double safmin;
-  if(nargin == 2)
+
+  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");
+  else if (!args(2).is_string ())
+    {
+      error ("qz: argument 3 must be a string");
+      return retval;
+    }
+  else
+    {
+      ord_job = args(2).string_value ();
 
-    // 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
+      if (ord_job[0] != 'N'
+	  && ord_job[0] != 'S'
+	  && ord_job[0] != 'B'
+	  && ord_job[0] != '+'
+	  && ord_job[0] != '-')
+	{
+	  error ("qz: illegal order option");
+	  return retval;
+	}
+
+      // overflow constant required by dlag2
+      F77_FCN (xdlamch, XDLAMCH) ("S", safmin, 1L);
 
-    // 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
+#ifdef DEBUG_EIG
+      cout << "qz: initial value of safmin=" << setiosflags (ios::scientific)
+	   << safmin << endl;
+#endif
 
-      F77_XFCN ( xdlamch, XDLAMCH, ("E", safmin, 1L));
+      // 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
 
-      #ifdef DEBUG_EIG
-      cout << "qz: safmin set to " << setiosflags(ios::scientific) 
-	<< safmin << endl;
-      #endif
+	  F77_FCN (xdlamch, XDLAMCH) ("E", safmin, 1L);
+
+#ifdef DEBUG_EIG
+	  cout << "qz: safmin set to " << setiosflags (ios::scientific)
+	       << safmin << endl;
+#endif
+	}
     }
-  }
 
-  #ifdef DEBUG
+#ifdef DEBUG
   cout << "qz: check argument 1" << endl;
-  #endif
+#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());
+  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;
-  }
+    {
+      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;
-  }
+    {
+      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 ()) 
+
+  if (args(0).is_complex_type ())
     caa = args(0).complex_matrix_value ();
-  else                            
+  else
     aa = args(0).matrix_value ();
-  if (error_state) 
+
+  if (error_state)
     return retval;
 
-  #ifdef DEBUG
+#ifdef DEBUG
   cout << "qz: check argument 2" << endl;
-  #endif
+#endif
 
   // Extract argument 2 (bb, or cbb if complex)
-  if( (nn != args(1).columns()) || (nn != args(1).rows() ))
-  {
-    gripe_nonconformant ();
-    return retval;
-  }
+  if ((nn != args(1).columns ()) || (nn != args(1).rows ()))
+    {
+      gripe_nonconformant ();
+      return retval;
+    }
+
   Matrix bb;
   ComplexMatrix cbb;
-  if (args(1).is_complex_type ()) 
+
+  if (args(1).is_complex_type ())
     cbb = args(1).complex_matrix_value ();
   else
     bb = args(1).matrix_value ();
-  if (error_state) 
+
+  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());
+
+  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.");
+  if (nargin == 3 && complex_case)
+    {
+      error ("qz: cannot re-order complex qz decomposition.");
+      return retval;
+    }
 
   // 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);
+  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');
+  char compq = (nargout >= 3 ? 'V' : 'N');
+  char 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++)
+  // 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 
+  // 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");
+  if (complex_case)
+    {
+      error ("Complex case not implemented yet");
+      return retval;
+    }
   else
-  {
-    #ifdef DEBUG
-    if(compq == 'V')
-      cout << "qz: performing balancing; QQ=" << endl << QQ << endl;
-    #endif
+    {
+#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)");
-  }
+      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)
+	{
+	  error ("unrecoverable error in qz (bal)");
+	  return retval;
+	}
+    }
 
   // 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));
+
+  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
+#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)");
+    if (f77_exception_encountered)
+      {
+	error ("unrecoverable error in qz (bal-L)");
+	return retval;
+      }
   }
 
   // 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));
+  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
+#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)");
-  }
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in qz (bal-R)");
+	  return retval;
+	}
+    }
 
   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");
-  }
+    {
+      // 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");
+      return retval;
+    }
   else  	// real matrices case
-  {
-    #ifdef DEBUG
-    cout << "qz: peforming qr decomposition of bb" << endl;
-    #endif
+    {
+#ifdef DEBUG
+      cout << "qz: peforming qr decomposition of bb" << endl;
+#endif
 
-    // compute the QR factorization of bb
-    QR bqr(bb);
+      // compute the QR factorization of bb
+      QR bqr (bb);
+
+#ifdef DEBUG
+      cout << "qz: qr (bb) done; now peforming qz decomposition" << endl;
+#endif
 
-    #ifdef DEBUG
-    cout << "qz: qr(bb) done; now peforming qz decomposition" << endl;
-    #endif
+      bb = bqr.R ();
+
+#ifdef DEBUG
+      cout << "qz: extracted bb" << 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;
 
-    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')
+	cout << "QQ =" << QQ << endl;
+#endif
 
-    if(compq == 'V') QQ = QQ*bqr.Q();
+      if (compq == 'V')
+	QQ = QQ*bqr.Q ();
 
-    #ifdef DEBUG
-    cout << "qz: precursors done..." << endl;
-    #endif
+#ifdef DEBUG
+      cout << "qz: precursors done..." << endl;
+#endif
 
+#ifdef DEBUG
+      cout << "qz: compq = " << compq << ", compz = " << compz << 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));
 
-    // 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)");
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in qz (dgghrd)");
+	  return retval;
+	}
 
-    // check if just computing generalized eigenvalues or if we're
-    // actually computing the decomposition
+      // 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)");
+      // 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)
+	{
+	  error ("unrecoverable error in qz (dhgeqz)");
+	  return retval;
+	}
+    }
 
   // 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
+  if (ord_job[0] != 'N')
     {
-      #ifdef DEBUG_SORT
-      cout << "qz: ordering eigenvalues: ord_job = " << ord_job[0] << endl;
-      #endif
+      if (complex_case)
+	{
+	  // probably not needed, but better be safe
+	  error ("qz: cannot re-order complex qz decomposition.");
+	  return retval;
+	}
+      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));
+	  // declared static to avoid vfork/long jump compiler complaints
+	  static sort_function sort_test;
+	  sort_test = NULL;
 
-      #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
+	  switch (ord_job[0])
+	    {
+	    case 'S':
+	      sort_test = &fin;
+	      break;
 
-      // 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
+	    case 'B':
+	      sort_test = &fout;
+	      break;
 
-          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));
+	    case '+':
+	      sort_test = &fcrhp;
+	      break;
 
-          #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
+	    case '-':
+	      sort_test = &folhp;
+	      break;
+
+	    default:
+	      // illegal order option (should never happen, since we
+	      // checked the options at the top).
+	      panic_impossible ();
+	      break;
       }
 
-      #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
+	  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;
+		    }
+		}
+
+	      // advance past this block
+	      jj += zcnt;
+	    }
+
+#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
+
+  if (nargout < 2 || nargout == 7 || (nargin == 3 && nargout == 4))
     {
-      #ifdef DEBUG
-      cout << "qz: computing generalized eigenvalues" << endl;
-      #endif
+      if (complex_case)
+	{
+	  error ("complex case not yet implemented");
+	  return retval;
+	}
+      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;
+	  // return finite generalized eigenvalues
+	  int cnt = 0;
+
+	  for (int ii = 0; ii < nn; ii++)
+	    if (betar(ii) != 0)
+	      cnt++;
+
+	  ComplexColumnVector tmp(cnt);
+
+	  for (int 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
+  if (nargout >= 5)
     {
-      #ifdef DEBUG
-      cout << "qz: computing  generalized eigenvectors" << endl;
-      #endif
+      char side = (nargout == 5 ? 'R' : 'B');	// which side to compute?
+      char 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");
+	  return retval;
+	}
+      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)
+	    {
+	      error ("unrecoverable error in qz (dtgevc)");
+	      return retval;
+	    }
 
-      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)
+		{
+		  for (int ii = 0; ii < nn; ii++)
+		    CVR(ii,jj) = VR(ii,jj);
+
+		  if (side == 'B')
+		    for (int ii = 0; ii < nn; ii++)
+		      CVL(ii,jj) = VL(ii,jj);
+		}
+	      else
+		{
+		  // double column; complex vector
+
+		  for (int 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));
+		    }
 
-      // 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)
-      }
+		  if (side == 'B')
+		    for (int 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));
+		      }
+		}
+
+	      // advance to next eigenvectors (if any)
+	      jj += cinc;
+	    }
+	}
     }
+
+  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.");
+      break;
   }
 
-  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
 
-  #ifdef DEBUG
-  cout << "qz: exiting (at long last)" << endl;
-  #endif
-  
   return retval;
 }