Mercurial > hg > octave-max

changeset 3185:9580887dd160

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1998-09-26 02:45:55 by jwe]
author jwe
date Sat, 26 Sep 1998 02:45:59 +0000
parents 3988763ec9d3
children edaa9a2d3d9c
files ChangeLog ROADMAP kpathsea/ChangeLog kpathsea/elt-dirs.c kpathsea/pathsearch.h liboctave/ChangeLog liboctave/Makefile.in liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/mx-ext.h liboctave/oct-env.cc liboctave/pathsearch.cc liboctave/pathsearch.h src/ChangeLog src/DLD-FUNCTIONS/balance.cc src/DLD-FUNCTIONS/qz.cc src/Makefile.in src/defaults.cc src/help.cc src/load-save.cc
diffstat 20 files changed, 1065 insertions(+), 840 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,7 @@
 Thu Sep 24 13:51:03 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* configure.in (AC_OUTPUT): Add libcruft/ordered-qz to the list.
+	Delete libcruft/balgen and libcruft/eispack from the list.
 
 Mon Aug 31 12:07:02 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
--- a/ROADMAP
+++ b/ROADMAP
@@ -14,16 +14,15 @@
   kpathsea      -- Karl Berry's path searching library
 
   libcruft      -- various numerical libraries (mostly Fortran)
-    balgen         * balancing subroutines for eigenvalue computations
     blas           * basic linear algebra subroutines
     dassl          * differential-algebraic system solver
-    eispack        * subroutines for eigenvalue problems
     fftpack        * subroutines for fast fourier transforms
     lapack         * linear algebra package
     linpack        * linear algebra package
     minpack        * nonlinear equation solver
     misc           * miscellaneous utilities
     odepack        * odinary differential equation solver
+    ordered-qz     * code for ordering eigenvalues for QZ factorization
     quadpack       * subroutines for numerical integration
     ranlib         * random number generators
     slatec-fn      * various special function subroutines
--- a/kpathsea/ChangeLog
+++ b/kpathsea/ChangeLog
@@ -1,3 +1,8 @@
+Fri Sep 25 13:28:54 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* elt-dirs.c (kpse_clear_dir_cache): New function.
+	* pathsearch.h (kpse_clear_dir_cache): Provide declaration.
+
 Mon Aug 31 12:07:02 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* config.sub: Accept armv4 everywhere arm is allowed.
--- a/kpathsea/elt-dirs.c
+++ b/kpathsea/elt-dirs.c
@@ -68,6 +68,33 @@
 static cache_entry *the_cache = NULL;
 static unsigned cache_length = 0;
 
+void
+kpse_clear_dir_cache P1H(void)
+{
+  while (cache_length > 0)
+    {
+      str_llist_type elt = *the_cache[--cache_length].value;
+
+      while (elt)
+	{
+	  str_llist_type next = STR_LLIST_NEXT (*elt);
+
+	  string s = STR_LLIST (*elt);
+
+	  if (s)
+	    free (s);
+
+	  free (elt);
+
+	  elt = next;
+	}
+    }
+
+  if (the_cache)
+    free (the_cache);
+
+  the_cache = NULL;
+}
 
 /* Associate KEY with VALUE.  We implement the cache as a simple linear
    list, since it's unlikely to ever be more than a dozen or so elements
--- a/kpathsea/pathsearch.h
+++ b/kpathsea/pathsearch.h
@@ -72,5 +72,8 @@
    all the filenames (or NULL if none), instead of taking the first.  */
 extern string *kpse_all_path_search P2H(const_string path, const_string name);
 
+/* Clear the directory cache.  */
+extern void kpse_clear_dir_cache P1H(void);
+
 #endif /* not KPATHSEA_PATHSEARCH_H */
 
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,5 +1,26 @@
+Fri Sep 25 14:26:44 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* oct-env.cc (octave_env::do_get_home_directory): 
+ 	If HOME can't be found, set it to "/".
+	(octave_env::do_get_user_name)
+	If user name can't be found, set it to "unknown".
+	(octave_env::do_get_host_name)
+	If host name can't be found, set it to "unknown".
+
+	* pathsearch.h (dir_path::rehash): New function.
+	* pathsearch.cc (dir_path::init): Clear kpathsea's internal
+	diretcory cache before doing initialization.
+
 Thu Sep 24 13:23:25 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* dMatrix.cc (Qzval): Delete.
+	(qzhes, qzit, qzval): Delete F77_FCN declarations.
+	* dMatrix.h (Qzval): Delete declaration.
+
+	* dbleGEPBAL.h, dbleGEPBAL.cc: Delete.
+	* Makefile.in (MATRIX_INC, MATRIX_SRC): Delete them from the lists.
+	* mx-ext.h: Don't include dbleGEPBAL.
+
 	* lo-ieee.cc (octave_ieee_init): For now, use X_CAST instead of
 	static_cast.
 
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -31,7 +31,7 @@
 	CmplxCHOL.h CmplxDET.h CmplxHESS.h CmplxLU.h CmplxQR.h \
 	CmplxQRP.h CmplxSCHUR.h CmplxSVD.h EIG.h boolMatrix.h chMatrix.h \
 	dColVector.h dDiagMatrix.h dMatrix.h dRowVector.h dbleAEPBAL.h \
-	dbleCHOL.h dbleDET.h dbleGEPBAL.h dbleHESS.h dbleLU.h dbleQR.h \
+	dbleCHOL.h dbleDET.h dbleHESS.h dbleLU.h dbleQR.h \
 	dbleQRP.h dbleSCHUR.h dbleSVD.h
 
 MX_OP_INC := mx-cdm-cm.h mx-cdm-cs.h mx-cdm-dm.h mx-cdm-m.h \
@@ -65,9 +65,8 @@
 	CmplxHESS.cc CmplxLU.cc CmplxQR.cc CmplxQRP.cc CmplxSCHUR.cc \
 	CmplxSVD.cc EIG.cc MArray-misc.cc boolMatrix.cc chMatrix.cc \
 	dColVector.cc dDiagMatrix.cc dMatrix.cc dRowVector.cc \
-	dbleAEPBAL.cc dbleCHOL.cc dbleDET.cc dbleGEPBAL.cc \
-	dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc dbleSCHUR.cc \
-	dbleSVD.cc
+	dbleAEPBAL.cc dbleCHOL.cc dbleDET.cc dbleHESS.cc dbleLU.cc \
+	dbleQR.cc dbleQRP.cc dbleSCHUR.cc dbleSVD.cc
 
 MX_OP_SRC := mx-cdm-cm.cc mx-cdm-cs.cc mx-cdm-dm.cc mx-cdm-m.cc \
 	mx-cdm-s.cc mx-cm-cdm.cc mx-cm-dm.cc mx-cm-m.cc mx-cm-s.cc \
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -99,17 +99,6 @@
   int F77_FCN (xdlange, XDLANGE) (const char*, const int&,
 				  const int&, const double*,
 				  const int&, double*, 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.
@@ -2935,94 +2924,6 @@
   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 = alfi.fortran_vec ();
-
-	      ColumnVector beta (a_nr);
-	      double *pbeta = beta.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(i) != 0)
-		      cnt++;
-
-		  ComplexColumnVector cx (cnt);
-
-		  cnt = 0;
-
-		  for (int i = 0; i < a_nr; i++)
-		    {
-		      if (beta(i) != 0)
-			{
-			  // Finite generalized eigenvalue.
-
-			  cx(cnt++) = Complex (alfr(i), alfi(i)) / beta(i);
-			}
-		    }
-
-		  retval = cx;
-		}
-	      else
-		(*current_liboctave_error_handler)
-		  ("qzval: trouble in qzit, info = %d", info);
-	    }
-	}
-      else
-	gripe_nonconformant ("qzval", a_nr, a_nc, b_nr, b_nc);
-    }
-  else
-    (*current_liboctave_error_handler) ("qzval: square matrices required");
-
-  return retval;
-}
-
 // matrix by matrix -> matrix operations
 
 Matrix
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -234,8 +234,6 @@
 
 extern Matrix Sylvester (const Matrix&, const Matrix&, const Matrix&);
 
-extern ComplexColumnVector Qzval (const Matrix& a, const Matrix& b);
-
 extern Matrix operator * (const Matrix& a, const Matrix& b);
 
 MS_CMP_OP_DECLS (Matrix, double)
--- a/liboctave/mx-ext.h
+++ b/liboctave/mx-ext.h
@@ -33,12 +33,6 @@
 #include "dbleDET.h"
 #include "CmplxDET.h"
 
-// Result of a GEP Balance operation
-// Note: currenlty only do balancing on real data.  Complex balancing
-// done on magnitudes of complex data.
-
-#include "dbleGEPBAL.h"
-
 // Result of a Cholesky Factorization
 
 #include "dbleCHOL.h"
--- a/liboctave/oct-env.cc
+++ b/liboctave/oct-env.cc
@@ -364,7 +364,7 @@
     {
       octave_passwd pw = octave_passwd::getpwuid (octave_syscalls::getuid ());
 
-      hd = pw ? pw.dir () : string ("I have no home!");
+      hd = pw ? pw.dir () : string ("/");
     }
 
   return hd;
@@ -380,7 +380,7 @@
     {
       octave_passwd pw = octave_passwd::getpwuid (octave_syscalls::getuid ());
 
-      user_name = pw ? string ("I have no name!") : pw.name ();
+      user_name = pw ? string ("unknown") : pw.name ();
     }
 
   return user_name;
@@ -398,7 +398,7 @@
 
       int status = gethostname (hostname, 255);
 
-      host_name = (status < 0) ? "I have no host!" : hostname;
+      host_name = (status < 0) ? "unknown" : hostname;
     }
 
   return host_name;
--- a/liboctave/pathsearch.cc
+++ b/liboctave/pathsearch.cc
@@ -159,6 +159,8 @@
       octave_kpathsea_initialized = true;
     }
 
+  kpse_clear_dir_cache ();
+
   char *tmp = kpse_path_expand (p_orig.c_str ());
   if (tmp)
     {
--- a/liboctave/pathsearch.h
+++ b/liboctave/pathsearch.h
@@ -71,6 +71,12 @@
 
   static void set_program_name (const string&);
 
+  void rehash (void)
+    {
+      initialized = false;
+      init ();
+    }
+
 private:
 
   // The colon separated list that we were given.
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,28 @@
+Fri Sep 25 11:50:44 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* load-save.cc (write_header): Rename from write_binary_header.
+	Also write header for Octave ASCII files.
+
+	* load-save.cc (Fsave): Implement -append option.
+
+	* defaults.cc (Frehash): New function.
+
+Fri Sep 25 11:50:44 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* help.cc (help_from_info): Improve error message in case that
+	info doesn't work.
+
 Thu Sep 24 10:48:12 1998  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* Makefile.in (DLD_XSRC): Replace qzval.cc with qz.cc
+
+	* DLD-FUNCTIONS/balance.cc: Update from A. S. Hodel
+	<scotte@eng.auburn.edu>.
+
+	* DLD-FUNCTIONS/qz.cc: New file.
+
+	* DLD-FUNCTIONS/qzval.cc: Delete.
+
 	* parse.y (plot_command1): Don't allow it to be empty.
 	(plot_command): Handle simple `PLOT' and `PLOT ranges' as special
 	cases here.
--- a/src/DLD-FUNCTIONS/balance.cc
+++ b/src/DLD-FUNCTIONS/balance.cc
@@ -42,24 +42,27 @@
 
 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( zggbal, ZGGBAL) ( const char* JOB,  const int& N,
-        Complex* A, const int& LDA, Complex* B, const int& LDB,
-        int& ILO, int & IHI, double* LSCALE,
-        double* RSCALE, double* WORK, int& INFO, long );
+  int F77_FCN (zggbal, ZGGBAL) (const char* JOB, const int& N,
+				Complex* A, const int& LDA, Complex* B,
+				const int& LDB, int& ILO, int& IHI,
+				double* LSCALE, double* RSCALE,
+				double* WORK, int& INFO, long);
 }
 
 DEFUN_DLD (balance, args, nargout,
-  "AA = balance (A [, OPT]) or [[DD,] AA] =  balance (A [, OPT])\n\
+  "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\
 \n\
 generalized eigenvalue problem:\n\
 \n\
@@ -77,7 +80,6 @@
 \n\
 [CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)")
 {
-
   octave_value_list retval;
 
   int nargin = args.length ();
@@ -89,7 +91,7 @@
     }
 
   // determine if it's AEP or GEP
-  int AEPcase = (nargin == 1 ? 1 : args(1).is_string() );
+  int AEPcase = nargin == 1 ? 1 : args(1).is_string ();
   string bal_job;
 
   // problem dimension
@@ -99,181 +101,224 @@
 
   if (arg_is_empty < 0)
     return retval;
+
   if (arg_is_empty > 0)
     return octave_value_list (2, Matrix ());
 
   if (nn != args(0).columns())
-  {
-    gripe_square_matrix_required ("balance");
-    return retval;
-  }
+    {
+      gripe_square_matrix_required ("balance");
+      return retval;
+    }
 
   // Extract argument 1 parameter for both AEP and GEP.
   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;
+
+  if (args(0).is_complex_type ())
+    caa = args(0).complex_matrix_value ();
+  else
+    aa = args(0).matrix_value ();
+
+  if (error_state)
+    return retval;
 
   // Treat AEP/GEP cases.
-  if(AEPcase)
-  {   
-    // Algebraic eigenvalue problem.
-    if(nargin == 1)
-      bal_job = "B";
-    else if(args(1).is_string())
-      bal_job = args(1).string_value();
-    // the next line should never execute, but better safe than sorry.
-    else error("balance: AEP argument 2 must be a string");
+  if (AEPcase)
+    {  
+      // Algebraic eigenvalue problem.
 
-    // balance the AEP
-    if (args(0).is_complex_type ())
-    {
-      ComplexAEPBALANCE result (caa, bal_job);
-
-      if (nargout == 0 || nargout == 1)
-        retval(0) = result.balanced_matrix ();
+      if (nargin == 1)
+	bal_job = "B";
+      else if (args(1).is_string ())
+	bal_job = args(1).string_value ();
       else
-      {
-        retval(1) = result.balanced_matrix ();
-        retval(0) = result.balancing_matrix ();
-      }
-    }
-    else
-    {
-      AEPBALANCE result (aa, bal_job);
+	{
+	  error ("balance: AEP argument 2 must be a string");
+	  return retval;
+	}
+
+      // balance the AEP
+      if (args(0).is_complex_type ())
+	{
+	  ComplexAEPBALANCE result (caa, bal_job);
 
-      if (nargout == 0 || nargout == 1)
-        retval(0) = result.balanced_matrix ();
+	  if (nargout == 0 || nargout == 1)
+	    retval(0) = result.balanced_matrix ();
+	  else
+	    {
+	      retval(1) = result.balanced_matrix ();
+	      retval(0) = result.balancing_matrix ();
+	    }
+	}
       else
-      {
-        retval(1) = result.balanced_matrix ();
-        retval(0) = result.balancing_matrix ();
-       }
+	{
+	  AEPBALANCE result (aa, bal_job);
+
+	  if (nargout == 0 || nargout == 1)
+	    retval(0) = result.balanced_matrix ();
+	  else
+	    {
+	      retval(1) = result.balanced_matrix ();
+	      retval(0) = result.balancing_matrix ();
+	    }
+	}
     }
-  }
-  //
-  // end of AEP case, now do GEP case
   else
-  {
-    // Generalized eigenvalue problem.
-    if(nargin == 2)
-      bal_job = "B";
-    else if(args(2).is_string())
-      bal_job = args(2).string_value();
-    else error("balance: GEP argument 3 must be a string");
-
-    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:
-    // first, declare variables used in both the real and complex case
-    int ilo, ihi, info;
-    RowVector lscale(nn), rscale(nn), work(6*nn);
-    char job = bal_job[0];
-    static int complex_case = (args(0).is_complex_type() 
-                       || args(1).is_complex_type());
-
-    // now balance
-    if (complex_case)
-    {
-      if (args(0).is_real_type ()) caa = aa;
-      if (args(1).is_real_type ()) cbb = bb;
-  
-      F77_XFCN( zggbal, ZGGBAL, ( &job, nn, caa.fortran_vec(), nn,
-          cbb.fortran_vec(), nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), work.fortran_vec(), info, 1L));
-    }
-    else          // real matrices case
     {
-      F77_XFCN( dggbal, DGGBAL, (&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 balance GEP");
-    }
-      
-    // Since we just want the balancing matrices, we can use dggbal
-    // for both the real and complex cases;
-    Matrix Pl(nn,nn), Pr(nn,nn);
-    for(int ii=0; ii < nn ; ii++)
-      for( int jj=0; jj < nn ; jj++)
-        Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0);
+      // Generalized eigenvalue problem.
+      if (nargin == 2)
+	bal_job = "B";
+      else if (args(2).is_string ())
+	bal_job = args(2).string_value ();
+      else
+	{
+	  error ("balance: GEP argument 3 must be a string");
+	  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 ())
+	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:
+      // first, declare variables used in both the real and complex case
+
+      int ilo, ihi, info;
+      RowVector lscale(nn), rscale(nn), work(6*nn);
+      char job = bal_job[0];
+
+      static int complex_case
+	= (args(0).is_complex_type () || args(1).is_complex_type ());
+
+      // now balance
+      if (complex_case)
+	{
+	  if (args(0).is_real_type ())
+	    caa = aa;
+
+	  if (args(1).is_real_type ())
+	    cbb = bb;
   
-    // left first
-    F77_XFCN( dggbak, DGGBAK, (&job, "L",
-          nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), nn, Pl.fortran_vec(),
-          nn, info, 1L, 1L));
-      
-    if(f77_exception_encountered)
-      (*current_liboctave_error_handler) 
-        ("unrecoverable error in balance GEP(L)");
-      
-    // then right
-    F77_XFCN(dggbak, DGGBAK, (&job, "R",
-          nn, ilo, ihi, lscale.fortran_vec(),
-          rscale.fortran_vec(), nn, Pr.fortran_vec(),
-          nn, info, 1L, 1L));
-    if(f77_exception_encountered)
-      (*current_liboctave_error_handler) 
-        ("unrecoverable error in balance GEP(R)");
+	  F77_XFCN (zggbal, ZGGBAL,
+		    (&job, nn, caa.fortran_vec(), nn,
+		     cbb.fortran_vec(), nn, ilo, ihi,
+		     lscale.fortran_vec(), rscale.fortran_vec(),
+		     work.fortran_vec(), info, 1L));
+
+	  if (f77_exception_encountered)
+	    {
+	      error ("unrecoverable error in balance GEP");
+	      return retval;
+	    }
+	}
+      else
+	{
+	  // real matrices case
 
-    switch (nargout)
-    {
-    case 0:
-    case 1:
-      warning ("balance: used GEP, should have two output arguments");
-      if(complex_case)
-        retval(0) = caa;
-      else
-        retval(0) = aa;
-      break;
+	  F77_XFCN (dggbal, DGGBAL,
+		    (&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 balance GEP");
+	      return retval;
+	    }
+	}
+      
+      // Since we just want the balancing matrices, we can use dggbal
+      // for both the real and complex cases.
+
+      Matrix Pl(nn,nn), Pr(nn,nn);
+
+      for (int ii = 0; ii < nn; ii++)
+	for (int jj = 0; jj < nn; jj++)
+	  Pl(ii,jj) = Pr(ii,jj) = (ii == jj ? 1.0 : 0.0);
+  
+      // left first
+      F77_XFCN (dggbak, DGGBAK,
+		(&job, "L", nn, ilo, ihi, lscale.fortran_vec(),
+		 rscale.fortran_vec(), nn, Pl.fortran_vec(),
+		 nn, info, 1L, 1L));
+      
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in balance GEP(L)");
+	  return retval;
+	}
+      
+      // then right
+      F77_XFCN (dggbak, DGGBAK,
+		(&job, "R", nn, ilo, ihi, lscale.fortran_vec(),
+		 rscale.fortran_vec(), nn, Pr.fortran_vec(),
+		 nn, info, 1L, 1L));
+
+      if (f77_exception_encountered)
+	{
+	  error ("unrecoverable error in balance GEP(R)");
+	  return retval;
+	}
 
-    case 2:
-      if(complex_case)
-      {
-        retval(1) = cbb;
-        retval(0) = caa;
-      }
-      else
-      {
-        retval(1) = bb;
-        retval(0) = aa;
-      }
-      break;
+      switch (nargout)
+	{
+	case 0:
+	case 1:
+	  warning ("balance: used GEP, should have two output arguments");
+	  if (complex_case)
+	    retval(0) = caa;
+	  else
+	    retval(0) = aa;
+	  break;
 
-    case 4:
-      if(complex_case)
-      {
-        retval(3) = cbb;
-        retval(2) = caa;
-      }
-      else
-      {
-        retval(3) = bb;
-        retval(2) = aa;
-      }
-      retval(1) = Pr;
-      retval(0) = Pl.transpose();  // so that aa_bal = cc*aa*dd, etc.
-      break;
+	case 2:
+	  if (complex_case)
+	    {
+	      retval(1) = cbb;
+	      retval(0) = caa;
+	    }
+	  else
+	    {
+	      retval(1) = bb;
+	      retval(0) = aa;
+	    }
+	  break;
 
-    default:
-      error ("balance: invalid number of output arguments");
-      break;
+	case 4:
+	  if (complex_case)
+	    {
+	      retval(3) = cbb;
+	      retval(2) = caa;
+	    }
+	  else
+	    {
+	      retval(3) = bb;
+	      retval(2) = aa;
+	    }
+	  retval(1) = Pr;
+	  retval(0) = Pl.transpose ();  // so that aa_bal = cc*aa*dd, etc.
+	  break;
+
+	default:
+	  error ("balance: invalid number of output arguments");
+	  break;
+	}
     }
-  }
+
   return retval;
 }
 
--- a/src/DLD-FUNCTIONS/qz.cc
+++ b/src/DLD-FUNCTIONS/qz.cc
@@ -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;
 }
 
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -43,7 +43,7 @@
 	det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \
 	fsolve.cc gammainc.cc getgrent.cc getpwent.cc getrusage.cc \
 	givens.cc hess.cc ifft.cc ifft2.cc inv.cc log.cc lpsolve.cc \
-	lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qzval.cc rand.cc \
+	lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qz.cc rand.cc \
 	schur.cc sort.cc svd.cc syl.cc time.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- a/src/defaults.cc
+++ b/src/defaults.cc
@@ -46,6 +46,7 @@
 #include "file-ops.h"
 #include "gripes.h"
 #include "help.h"
+#include "oct-obj.h"
 #include "ov.h"
 #include "toplev.h"
 #include "variables.h"
@@ -427,6 +428,16 @@
     "Octave version");
 }
 
+DEFUN (rehash, , ,
+  "rehash (): reinitialize LOADPATH directory cache")
+{
+  octave_value_list retval;
+
+  Vload_path_dir_path.rehash ();
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/src/help.cc
+++ b/src/help.cc
@@ -640,7 +640,8 @@
 		}
 	      else
 		{
-		  error ("help: unable to find info!");
+		  error ("help: unable to find info");
+		  error ("help: you need info 2.18 or later (texinfo 3.12)");
 		  break;
 		}
 	    }
--- a/src/load-save.cc
+++ b/src/load-save.cc
@@ -27,6 +27,7 @@
 #include <cfloat>
 #include <cstring>
 #include <cctype>
+#include <ctime>
 
 #include <string>
 
@@ -41,6 +42,7 @@
 #include "glob-match.h"
 #include "lo-mappers.h"
 #include "mach-info.h"
+#include "oct-env.h"
 #include "str-vec.h"
 
 #include "defun.h"
@@ -55,6 +57,7 @@
 #include "unwind-prot.h"
 #include "utils.h"
 #include "variables.h"
+#include "version.h"
 
 // The default output format.  May be one of "binary", "text", or
 // "mat-binary".
@@ -2366,28 +2369,53 @@
 }
 
 static void
-write_binary_header (ostream& os, load_save_format format)
+write_header (ostream& os, load_save_format format)
 {
-  if (format == LS_BINARY)
+  switch (format)
     {
-      os << (oct_mach_info::words_big_endian ()
-	     ? "Octave-1-B" : "Octave-1-L");
-
-      oct_mach_info::float_format flt_fmt =
-	oct_mach_info::native_float_format ();
-
-      char tmp = (char) float_format_to_mopt_digit (flt_fmt);
-
-      os.write (&tmp, 1);
+    case LS_BINARY:
+      {
+	os << (oct_mach_info::words_big_endian ()
+	       ? "Octave-1-B" : "Octave-1-L");
+
+	oct_mach_info::float_format flt_fmt =
+	  oct_mach_info::native_float_format ();
+
+	char tmp = (char) float_format_to_mopt_digit (flt_fmt);
+
+	os.write (&tmp, 1);
+      }
+    break;
+
+    case LS_ASCII:
+      {
+	time_t now = time (0);
+
+	string time_string = asctime (gmtime (&now));
+	time_string = time_string.substr (0, time_string.length () - 1);
+
+	os << "# Created by Octave " OCTAVE_VERSION ", "
+	   << time_string
+	   << " <"
+	   << octave_env::get_user_name ()
+	   << "@"
+	   << octave_env::get_host_name ()
+	   << ">" << "\n";
+      }
+    break;
+
+    default:
+      break;
     }
 }
 
 static void
 save_vars (const string_vector& argv, int argv_idx, int argc,
 	   ostream& os, bool save_builtins, load_save_format fmt,
-	   bool save_as_floats) 
+	   bool save_as_floats, bool write_header_info)
 {
-  write_binary_header (os, fmt);
+  if (write_header_info)
+    write_header (os, fmt);
 
   if (argv_idx == argc)
     {
@@ -2424,7 +2452,7 @@
 
   if (file)
     {
-      save_vars (string_vector (), 0, 0, file, false, format, false);
+      save_vars (string_vector (), 0, 0, file, false, format, false, true);
       message (0, "save to `%s' complete", fname);
     }
   else
@@ -2432,7 +2460,7 @@
 }
 
 DEFUN_TEXT (save, args, ,
-  "save [-ascii] [-binary] [-float-binary] [-mat-binary] \n\
+  "save [-append] [-ascii] [-binary] [-float-binary] [-mat-binary] \n\
      [-save-builtins] file [pattern ...]\n\
 \n\
 save variables in a file")
@@ -2455,9 +2483,15 @@
 
   load_save_format format = get_default_save_format ();
 
+  bool append = false;
+
   int i;
   for (i = 1; i < argc; i++)
     {
+      if (argv[i] == "-append")
+	{
+	  append = true;
+	}
       if (argv[i] == "-ascii" || argv[i] == "-a")
 	{
 	  format = LS_ASCII;
@@ -2503,7 +2537,7 @@
       // in a octave_value (string)?
 
       save_vars (argv, i, argc, octave_stdout, save_builtins, format,
-		 save_as_floats);
+		 save_as_floats, true);
     }
 
   // Guard against things like `save a*', which are probably mistakes...
@@ -2519,16 +2553,21 @@
 
       i++;
 
-      unsigned mode = ios::out|ios::trunc;
+      unsigned mode = ios::out;
       if (format == LS_BINARY || format == LS_MAT_BINARY)
 	mode |= ios::bin;
 
+      mode |= append ? ios::ate : ios::trunc;
+
       ofstream file (fname.c_str (), mode);
 
       if (file)
 	{
+	  bool write_header_info
+	    = ((file.rdbuf ())->seekoff (0, ios::cur) == 0);
+
 	  save_vars (argv, i, argc, file, save_builtins, format,
-		     save_as_floats);
+		     save_as_floats, write_header_info);
 	}
       else
 	{