# HG changeset patch
# User jwe
# Date 906694642 0
# Node ID 5edc539c2f80bd8a4516ef469344e29169d1a2f0
# Parent  936b9ae8f7d245387e3161d3c00987f995d3f0bf
[project @ 1998-09-25 03:37:22 by jwe]

diff --git a/src/DLD-FUNCTIONS/qz.cc b/src/DLD-FUNCTIONS/qz.cc
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: ***
+*/