# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1459288028 -7200
#      Tue Mar 29 23:47:08 2016 +0200
# Node ID 40f344f2ad2e168812b76d1d49ea5671245e1ae0
# Parent  a7669bade52fbbe92969c3d465eeaf5210df1543
add preliminary balancing to eig

* EIG.cc: changing dgeev to dgeevx and zgeev to zgeevx
* EIG.h: new parameter: balance
* fEIG.cc: changing sgeev to sgeevx and cgeev to cgeevx
* fEIG.h: new parameter: balance
* eig.cc: add option to choose balance

diff -r a7669bade52f -r 40f344f2ad2e libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Sat Mar 12 19:41:03 2016 +0100
+++ b/libinterp/corefcn/eig.cc	Tue Mar 29 23:47:08 2016 +0200
@@ -75,26 +75,40 @@ The eigenvalues returned by @code{eig} a
   if (! arg_a.is_double_type () && ! arg_a.is_single_type ())
     err_wrong_type_arg ("eig", arg_a);
 
+  bool balance = true;
   if (nargin == 2)
     {
-      arg_b = args(1);
-      nr_b = arg_b.rows ();
-      nc_b = arg_b.columns ();
+      if(args(1).is_string ())
+        {
+          std::string a1s = args(1).string_value ();
+          if(a1s != "balance" && a1s != "nobalance")
+            error ("eig: unexpected second input: %s", a1s.c_str ());
+          balance = a1s == "balance";
+        }
+      else
+        {
+          arg_b = args(1);
+          nr_b = arg_b.rows ();
+          nc_b = arg_b.columns ();
 
-      arg_is_empty = empty_arg ("eig", nr_b, nc_b);
-      if (arg_is_empty < 0)
-        return retval;
-      else if (arg_is_empty > 0)
-        return ovl (2, Matrix ());
+          arg_is_empty = empty_arg ("eig", nr_b, nc_b);
+          if (arg_is_empty < 0)
+            return retval;
+          else if (arg_is_empty > 0)
+            return ovl (2, Matrix ());
 
-      if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
-        err_wrong_type_arg ("eig", arg_b);
+          if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
+            err_wrong_type_arg ("eig", arg_b);
+        }
     }
 
   if (nr_a != nc_a)
     err_square_matrix_required ("eig", "A");
-
-  if (nargin == 2 && nr_b != nc_b)
+  
+  // determine if it's AEP or GEP
+  bool AEPcase = nargin == 1 || args(1).is_string ();
+  
+  if (!AEPcase && nr_b != nc_b)
     err_square_matrix_required ("eig", "B");
 
   Matrix tmp_a, tmp_b;
@@ -106,22 +120,22 @@ The eigenvalues returned by @code{eig} a
     {
       FloatEIG result;
 
-      if (nargin == 1)
+      if (AEPcase)
         {
           if (arg_a.is_real_type ())
             {
               ftmp_a = arg_a.float_matrix_value ();
 
-              result = FloatEIG (ftmp_a, nargout > 1);
+              result = FloatEIG (ftmp_a, nargout > 1, balance);
             }
           else
             {
               fctmp_a = arg_a.float_complex_matrix_value ();
 
-              result = FloatEIG (fctmp_a, nargout > 1);
+              result = FloatEIG (fctmp_a, nargout > 1, balance);
             }
         }
-      else if (nargin == 2)
+      else
         {
           if (arg_a.is_real_type () && arg_b.is_real_type ())
             {
@@ -155,22 +169,22 @@ The eigenvalues returned by @code{eig} a
     {
       EIG result;
 
-      if (nargin == 1)
+      if (AEPcase)
         {
           if (arg_a.is_real_type ())
             {
               tmp_a = arg_a.matrix_value ();
 
-              result = EIG (tmp_a, nargout > 1);
+              result = EIG (tmp_a, nargout > 1, balance);
             }
           else
             {
               ctmp_a = arg_a.complex_matrix_value ();
 
-              result = EIG (ctmp_a, nargout > 1);
+              result = EIG (ctmp_a, nargout > 1, balance);
             }
         }
-      else if (nargin == 2)
+      else
         {
           if (arg_a.is_real_type () && arg_b.is_real_type ())
             {
diff -r a7669bade52f -r 40f344f2ad2e liboctave/numeric/EIG.cc
--- a/liboctave/numeric/EIG.cc	Sat Mar 12 19:41:03 2016 +0100
+++ b/liboctave/numeric/EIG.cc	Tue Mar 29 23:47:08 2016 +0200
@@ -32,26 +32,40 @@ along with Octave; see the file COPYING.
 extern "C"
 {
   F77_RET_T
-  F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL,
-                           F77_CONST_CHAR_ARG_DECL,
-                           const octave_idx_type&, double*,
-                           const octave_idx_type&, double*, double*,
-                           double*, const octave_idx_type&, double*,
-                           const octave_idx_type&, double*,
-                           const octave_idx_type&, octave_idx_type&
-                           F77_CHAR_ARG_LEN_DECL
-                           F77_CHAR_ARG_LEN_DECL);
+  F77_FUNC (dgeevx, DGEEVX) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*,
+                             double*, const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, double*, double&,
+                             double*, double*, double*,
+                             const octave_idx_type&, octave_idx_type*,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
-  F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL,
-                           F77_CONST_CHAR_ARG_DECL,
-                           const octave_idx_type&, Complex*,
-                           const octave_idx_type&, Complex*,
-                           Complex*, const octave_idx_type&, Complex*,
-                           const octave_idx_type&, Complex*,
-                           const octave_idx_type&, double*, octave_idx_type&
-                           F77_CHAR_ARG_LEN_DECL
-                           F77_CHAR_ARG_LEN_DECL);
+  F77_FUNC (zgeevx, ZGEEVX) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*,
+                             Complex*, const octave_idx_type&, Complex*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, double*, double&,
+                             double*, double*, Complex*,
+                             const octave_idx_type&, double*,
+                             octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
   F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL,
@@ -137,7 +151,7 @@ extern "C"
 }
 
 octave_idx_type
-EIG::init (const Matrix& a, bool calc_ev)
+EIG::init (const Matrix& a, bool calc_ev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -171,33 +185,61 @@ EIG::init (const Matrix& a, bool calc_ev
 
   double *dummy = 0;
   octave_idx_type idummy = 1;
-
-  F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pwr, pwi, dummy,
-                           idummy, pvr, n, &dummy_work, lwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+  
+  Array<double> scale (dim_vector (n, 1));
+  double *pscale = scale.fortran_vec ();
+  
+  double abnrm;
+  
+  Array<double> rconde (dim_vector (n, 1));
+  double *prconde = rconde.fortran_vec ();
+  
+  Array<double> rcondv (dim_vector (n, 1));
+  double *prcondv = rcondv.fortran_vec ();
+  
+  octave_idx_type dummy_iwork;
+        
+  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pwr, pwi, dummy,
+                             idummy, pvr, n, ilo, ihi, pscale,
+                             abnrm, prconde, prcondv, &dummy_work,
+                             lwork, &dummy_iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info != 0)
-    (*current_liboctave_error_handler) ("dgeev workspace query failed");
+    (*current_liboctave_error_handler) ("dgeevx workspace query failed");
 
   lwork = static_cast<octave_idx_type> (dummy_work);
   Array<double> work (dim_vector (lwork, 1));
   double *pwork = work.fortran_vec ();
-
-  F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pwr, pwi, dummy,
-                           idummy, pvr, n, pwork, lwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  
+  F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pwr, pwi, dummy,
+                             idummy, pvr, n, ilo, ihi, pscale,
+                             abnrm, prconde, prcondv, pwork,
+                             lwork, &dummy_iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
+    (*current_liboctave_error_handler) ("unrecoverable error in dgeevx");
 
   if (info > 0)
-    (*current_liboctave_error_handler) ("dgeev failed to converge");
+    (*current_liboctave_error_handler) ("dgeevx failed to converge");
 
   lambda.resize (n);
   octave_idx_type nvr = calc_ev ? n : 0;
@@ -284,7 +326,7 @@ EIG::symmetric_init (const Matrix& a, bo
 }
 
 octave_idx_type
-EIG::init (const ComplexMatrix& a, bool calc_ev)
+EIG::init (const ComplexMatrix& a, bool calc_ev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -320,32 +362,58 @@ EIG::init (const ComplexMatrix& a, bool 
   Complex *dummy = 0;
   octave_idx_type idummy = 1;
 
-  F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pw, dummy, idummy,
-                           pv, n, &dummy_work, lwork, prwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+  
+  Array<double> scale (dim_vector (n, 1));
+  double *pscale = scale.fortran_vec ();
+  
+  double abnrm;
+  
+  Array<double> rconde (dim_vector (n, 1));
+  double *prconde = rconde.fortran_vec ();
+  
+  Array<double> rcondv (dim_vector (n, 1));
+  double *prcondv = rcondv.fortran_vec ();
+   
+  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pw, dummy, idummy,
+                             pv, n, ilo, ihi, pscale, abnrm,
+                             prconde, prcondv,
+                             &dummy_work, lwork, prwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info != 0)
-    (*current_liboctave_error_handler) ("zgeev workspace query failed");
+    (*current_liboctave_error_handler) ("zgeevx workspace query failed");
 
   lwork = static_cast<octave_idx_type> (dummy_work.real ());
   Array<Complex> work (dim_vector (lwork, 1));
   Complex *pwork = work.fortran_vec ();
 
-  F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pw, dummy, idummy,
-                           pv, n, pwork, lwork, prwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pw, dummy, idummy,
+                             pv, n, ilo, ihi, pscale, abnrm,
+                             prconde, prcondv,
+                             pwork, lwork, prwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
+    (*current_liboctave_error_handler) ("unrecoverable error in zgeevx");
 
   if (info > 0)
-    (*current_liboctave_error_handler) ("zgeev failed to converge");
+    (*current_liboctave_error_handler) ("zgeevx failed to converge");
 
   lambda = w;
   v = vtmp;
diff -r a7669bade52f -r 40f344f2ad2e liboctave/numeric/EIG.h
--- a/liboctave/numeric/EIG.h	Sat Mar 12 19:41:03 2016 +0100
+++ b/liboctave/numeric/EIG.h	Tue Mar 29 23:47:08 2016 +0200
@@ -42,16 +42,18 @@ public:
 
   EIG (void) : lambda (), v () { }
 
-  EIG (const Matrix& a, bool calc_eigenvectors = true)
+  EIG (const Matrix& a, bool calc_eigenvectors = true,
+       bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_eigenvectors, balance);
   }
 
-  EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
+  EIG (const Matrix& a, octave_idx_type& info,
+       bool calc_eigenvectors = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_eigenvectors, balance);
   }
 
   EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true)
@@ -67,17 +69,18 @@ public:
     info = init (a, b, calc_eigenvectors);
   }
 
-  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
+  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true,
+       bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_eigenvectors, balance);
   }
 
   EIG (const ComplexMatrix& a, octave_idx_type& info,
-       bool calc_eigenvectors = true)
+       bool calc_eigenvectors = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_eigenvectors, balance);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
@@ -120,12 +123,14 @@ private:
   ComplexColumnVector lambda;
   ComplexMatrix v;
 
-  octave_idx_type init (const Matrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const Matrix& a, bool calc_eigenvectors,
+                        bool balance);
 
   octave_idx_type init (const Matrix& a, const Matrix& b,
                         bool calc_eigenvectors);
 
-  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors,
+                        bool balance);
 
   octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b,
                         bool calc_eigenvectors);
diff -r a7669bade52f -r 40f344f2ad2e liboctave/numeric/fEIG.cc
--- a/liboctave/numeric/fEIG.cc	Sat Mar 12 19:41:03 2016 +0100
+++ b/liboctave/numeric/fEIG.cc	Tue Mar 29 23:47:08 2016 +0200
@@ -32,26 +32,38 @@ along with Octave; see the file COPYING.
 extern "C"
 {
   F77_RET_T
-  F77_FUNC (sgeev, SGEEV) (F77_CONST_CHAR_ARG_DECL,
-                           F77_CONST_CHAR_ARG_DECL,
-                           const octave_idx_type&, float*,
-                           const octave_idx_type&, float*, float*, float*,
-                           const octave_idx_type&, float*,
-                           const octave_idx_type&, float*,
-                           const octave_idx_type&, octave_idx_type&
-                           F77_CHAR_ARG_LEN_DECL
-                           F77_CHAR_ARG_LEN_DECL);
-
+  F77_FUNC (sgeevx, SGEEVX) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, float*, float&, float*,
+                             float*, float*, const octave_idx_type&,
+                             octave_idx_type*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+                      
   F77_RET_T
-  F77_FUNC (cgeev, CGEEV) (F77_CONST_CHAR_ARG_DECL,
-                           F77_CONST_CHAR_ARG_DECL,
-                           const octave_idx_type&, FloatComplex*,
-                           const octave_idx_type&, FloatComplex*, FloatComplex*,
-                           const octave_idx_type&, FloatComplex*,
-                           const octave_idx_type&, FloatComplex*,
-                           const octave_idx_type&, float*, octave_idx_type&
-                           F77_CHAR_ARG_LEN_DECL
-                           F77_CHAR_ARG_LEN_DECL);
+  F77_FUNC (cgeevx, CGEEVX) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, float*, float&, float*,
+                             float*, FloatComplex*, const octave_idx_type&,
+                             float*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
   F77_FUNC (ssyev, SSYEV) (F77_CONST_CHAR_ARG_DECL,
@@ -134,7 +146,7 @@ extern "C"
 }
 
 octave_idx_type
-FloatEIG::init (const FloatMatrix& a, bool calc_ev)
+FloatEIG::init (const FloatMatrix& a, bool calc_ev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -169,32 +181,60 @@ FloatEIG::init (const FloatMatrix& a, bo
   float *dummy = 0;
   octave_idx_type idummy = 1;
 
-  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pwr, pwi, dummy,
-                           idummy, pvr, n, &dummy_work, lwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+  
+  Array<float> scale (dim_vector (n, 1));
+  float *pscale = scale.fortran_vec ();
+  
+  float abnrm;
+  
+  Array<float> rconde (dim_vector (n, 1));
+  float *prconde = rconde.fortran_vec ();
+  
+  Array<float> rcondv (dim_vector (n, 1));
+  float *prcondv = rcondv.fortran_vec ();
+  
+  octave_idx_type dummy_iwork;
+  
+  F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pwr, pwi, dummy,
+                             idummy, pvr, n,
+                             ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             &dummy_work, lwork, &dummy_iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));                        
 
   if (info != 0)
-    (*current_liboctave_error_handler) ("sgeev workspace query failed");
+    (*current_liboctave_error_handler) ("sgeevx workspace query failed");
 
   lwork = static_cast<octave_idx_type> (dummy_work);
   Array<float> work (dim_vector (lwork, 1));
   float *pwork = work.fortran_vec ();
-
-  F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pwr, pwi, dummy,
-                           idummy, pvr, n, pwork, lwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  
+  F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pwr, pwi, dummy,
+                             idummy, pvr, n,
+                             ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             pwork, lwork, &dummy_iwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in sgeev");
+    (*current_liboctave_error_handler) ("unrecoverable error in sgeevx");
 
   if (info > 0)
-    (*current_liboctave_error_handler) ("sgeev failed to converge");
+    (*current_liboctave_error_handler) ("sgeevx failed to converge");
 
   lambda.resize (n);
   v.resize (nvr, nvr);
@@ -280,7 +320,7 @@ FloatEIG::symmetric_init (const FloatMat
 }
 
 octave_idx_type
-FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev)
+FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -316,32 +356,56 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplex *dummy = 0;
   octave_idx_type idummy = 1;
 
-  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pw, dummy, idummy,
-                           pv, n, &dummy_work, lwork, prwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+  
+  Array<float> scale (dim_vector (n, 1));
+  float *pscale = scale.fortran_vec ();
+  
+  float abnrm;
+  
+  Array<float> rconde (dim_vector (n, 1));
+  float *prconde = rconde.fortran_vec ();
+  
+  Array<float> rcondv (dim_vector (n, 1));
+  float *prcondv = rcondv.fortran_vec ();
+  
+  F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pw, dummy, idummy,
+                             pv, n, ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             &dummy_work, lwork, prwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
 
   if (info != 0)
-    (*current_liboctave_error_handler) ("cgeev workspace query failed");
+    (*current_liboctave_error_handler) ("cgeevx workspace query failed");
 
   lwork = static_cast<octave_idx_type> (dummy_work.real ());
   Array<FloatComplex> work (dim_vector (lwork, 1));
   FloatComplex *pwork = work.fortran_vec ();
 
-  F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, tmp_data, n, pw, dummy, idummy,
-                           pv, n, pwork, lwork, prwork, info
-                           F77_CHAR_ARG_LEN (1)
-                           F77_CHAR_ARG_LEN (1)));
+  F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, tmp_data, n, pw, dummy, idummy,
+                             pv, n, ilo, ihi, pscale,abnrm, prconde, prcondv,
+                             pwork, lwork, prwork, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));  
 
   if (info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in cgeev");
+    (*current_liboctave_error_handler) ("unrecoverable error in cgeevx");
 
   if (info > 0)
-    (*current_liboctave_error_handler) ("cgeev failed to converge");
+    (*current_liboctave_error_handler) ("cgeevx failed to converge");
 
   lambda = w;
   v = vtmp;
diff -r a7669bade52f -r 40f344f2ad2e liboctave/numeric/fEIG.h
--- a/liboctave/numeric/fEIG.h	Sat Mar 12 19:41:03 2016 +0100
+++ b/liboctave/numeric/fEIG.h	Tue Mar 29 23:47:08 2016 +0200
@@ -43,17 +43,18 @@ public:
   FloatEIG (void)
     : lambda (), v () { }
 
-  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true)
+  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true,
+            bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_eigenvectors, balance);
   }
 
   FloatEIG (const FloatMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true)
+            bool calc_eigenvectors = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_eigenvectors, balance);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
@@ -63,24 +64,25 @@ public:
     init (a, b, calc_eigenvectors);
   }
 
-  FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info,
-            bool calc_eigenvectors = true)
+  FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
+            octave_idx_type& info, bool calc_eigenvectors = true)
     : lambda (), v ()
   {
     info = init (a, b, calc_eigenvectors);
   }
 
-  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true)
+  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true,
+            bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_eigenvectors, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true)
+            bool calc_eigenvectors = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_eigenvectors, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
@@ -122,10 +124,12 @@ private:
   FloatComplexColumnVector lambda;
   FloatComplexMatrix v;
 
-  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors,
+                        bool balance);
   octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b,
                         bool calc_eigenvectors);
-  octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors,
+                        bool balance);
   octave_idx_type init (const FloatComplexMatrix& a,
                         const FloatComplexMatrix& b, bool calc_eigenvectors);
 
# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1466106663 -7200
#      Thu Jun 16 21:51:03 2016 +0200
# Node ID 559b56d48ceb057179dac882ed54b1ac5bf55e23
# Parent  40f344f2ad2e168812b76d1d49ea5671245e1ae0
eig.cc: Add tests

diff -r 40f344f2ad2e -r 559b56d48ceb libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Tue Mar 29 23:47:08 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Thu Jun 16 21:51:03 2016 +0200
@@ -227,8 +227,16 @@ The eigenvalues returned by @code{eig} a
 %! assert (d, [-1, 0; 0, 3], sqrt (eps));
 %! assert (v, [-x, x; x, x], sqrt (eps));
 
+%!test
+%! [v, d] = eig ([1, 2; 2, 1], 'balance');
+%! x = 1 / sqrt (2);
+%! assert (d, [-1, 0; 0, 3], sqrt (eps));
+%! assert (v, [-x, x; x, x], sqrt (eps));
+
 %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single")))
 
+%!assert (eig (single ([1, 2; 2, 1]), 'balance'), single ([-1; 3]), sqrt (eps ("single")))
+
 %!test
 %! [v, d] = eig (single ([1, 2; 2, 1]));
 %! x = single (1 / sqrt (2));
@@ -236,6 +244,12 @@ The eigenvalues returned by @code{eig} a
 %! assert (v, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
+%! [v, d] = eig (single ([1, 2; 2, 1]), 'balance');
+%! x = single (1 / sqrt (2));
+%! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")));
+%! assert (v, [-x, x; x, x], sqrt (eps ("single")));
+
+%!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
 %! [v, d] = eig (A, B);
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
@@ -306,6 +320,14 @@ The eigenvalues returned by @code{eig} a
 %!error <EIG requires same size matrices> eig ([1, 2; 3, 4], 2)
 %!error <must be a square matrix> eig ([1, 2; 3, 4; 5, 6])
 %!error <wrong type argument> eig ("abcd")
-%!error <wrong type argument> eig ([1 2 ; 2 3], "abcd")
+%!error <eig: unexpected second input: abcd> eig ([1 2 ; 2 3], "abcd")
 %!error <wrong type argument> eig (false, [1 2 ; 2 3])
-*/
+%!test
+%! A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1];
+%! [V1, D1] = eig (A);
+%! [V2, D2] = eig (A, 'balance');
+%! [V3, D3] = eig (A, 'nobalance');
+%! assert (V1, V2);
+%! assert (D1, D2);
+%! assert (isequal (V2, V3), false);
+*/
\ No newline at end of file
 
 
# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1466461551 -7200
#      Tue Jun 21 00:25:51 2016 +0200
# Node ID 05420af0a41bae4f311bc1db9024aa645c73e2c9
# Parent  78ba6d182780423b994b94dc584d2f198a9b18c1
add left eigenvector computation

* EIG.cc: add calc_lev
* fEIG.cc: add calc_lev
* EIG.h: new matrix for left eigenvectors: w; rename eigenvectors to right_eigenvectors
* fEIG.h: new matrix for left eigenvectors: w; rename eigenvectors to right_eigenvectors
* eig.cc: return left eigenvectors
* xpow.cc: rename eigenvectors to right_eigenvectors
* __qp__.cc: rename eigenvectors to right_eigenvectors

diff -r 78ba6d182780 -r 05420af0a41b libinterp/corefcn/__qp__.cc
--- a/libinterp/corefcn/__qp__.cc	Thu Jun 16 22:06:03 2016 +0200
+++ b/libinterp/corefcn/__qp__.cc	Tue Jun 21 00:25:51 2016 +0200
@@ -148,7 +148,7 @@ qp (const Matrix& H, const ColumnVector&
     }
 
   ColumnVector eigenvalH = real (eigH.eigenvalues ());
-  Matrix eigenvecH = real (eigH.eigenvectors ());
+  Matrix eigenvecH = real (eigH.right_eigenvectors ());
   double minReal = eigenvalH.min ();
   octave_idx_type indminR = 0;
   for (octave_idx_type i = 0; i < n; i++)
@@ -296,7 +296,7 @@ qp (const Matrix& H, const ColumnVector&
                 }
 
               ColumnVector eigenvalrH = real (eigrH.eigenvalues ());
-              Matrix eigenvecrH = real (eigrH.eigenvectors ());
+              Matrix eigenvecrH = real (eigrH.right_eigenvectors ());
               double mRrH = eigenvalrH.min ();
               indminR = 0;
               for (octave_idx_type i = 0; i < n; i++)
diff -r 78ba6d182780 -r 05420af0a41b libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Thu Jun 16 22:06:03 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Tue Jun 21 00:25:51 2016 +0200
@@ -119,20 +119,19 @@ The eigenvalues returned by @code{eig} a
   if (arg_a.is_single_type ())
     {
       FloatEIG result;
-
       if (AEPcase)
         {
           if (arg_a.is_real_type ())
             {
               ftmp_a = arg_a.float_matrix_value ();
 
-              result = FloatEIG (ftmp_a, nargout > 1, balance);
+              result = FloatEIG (ftmp_a, nargout > 1, nargout > 2, balance);
             }
           else
             {
               fctmp_a = arg_a.float_complex_matrix_value ();
 
-              result = FloatEIG (fctmp_a, nargout > 1, balance);
+              result = FloatEIG (fctmp_a, nargout > 1, nargout > 2, balance);
             }
         }
       else
@@ -142,14 +141,14 @@ The eigenvalues returned by @code{eig} a
               ftmp_a = arg_a.float_matrix_value ();
               ftmp_b = arg_b.float_matrix_value ();
 
-              result = FloatEIG (ftmp_a, ftmp_b, nargout > 1);
+              result = FloatEIG (ftmp_a, ftmp_b, nargout > 1, nargout > 2);
             }
           else
             {
               fctmp_a = arg_a.float_complex_matrix_value ();
               fctmp_b = arg_b.float_complex_matrix_value ();
 
-              result = FloatEIG (fctmp_a, fctmp_b, nargout > 1);
+              result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2);
             }
         }
 
@@ -157,12 +156,18 @@ The eigenvalues returned by @code{eig} a
         {
           retval = ovl (result.eigenvalues ());
         }
-      else
+      else if (nargout == 2)
         {
           // Blame it on Matlab.
           FloatComplexDiagMatrix d (result.eigenvalues ());
 
-          retval = ovl (result.eigenvectors (), d);
+          retval = ovl (result.right_eigenvectors (), d);
+        }
+      else
+        {
+          FloatComplexDiagMatrix d (result.eigenvalues ());
+
+          retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
         }
     }
   else
@@ -175,13 +180,13 @@ The eigenvalues returned by @code{eig} a
             {
               tmp_a = arg_a.matrix_value ();
 
-              result = EIG (tmp_a, nargout > 1, balance);
+              result = EIG (tmp_a, nargout > 1, nargout > 2, balance);
             }
           else
             {
               ctmp_a = arg_a.complex_matrix_value ();
 
-              result = EIG (ctmp_a, nargout > 1, balance);
+              result = EIG (ctmp_a, nargout > 1, nargout > 2, balance);
             }
         }
       else
@@ -191,14 +196,14 @@ The eigenvalues returned by @code{eig} a
               tmp_a = arg_a.matrix_value ();
               tmp_b = arg_b.matrix_value ();
 
-              result = EIG (tmp_a, tmp_b, nargout > 1);
+              result = EIG (tmp_a, tmp_b, nargout > 1, nargout > 2);
             }
           else
             {
               ctmp_a = arg_a.complex_matrix_value ();
               ctmp_b = arg_b.complex_matrix_value ();
 
-              result = EIG (ctmp_a, ctmp_b, nargout > 1);
+              result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2);
             }
         }
 
@@ -206,12 +211,19 @@ The eigenvalues returned by @code{eig} a
         {
           retval = ovl (result.eigenvalues ());
         }
-      else
+      else if(nargout == 2)
         {
           // Blame it on Matlab.
           ComplexDiagMatrix d (result.eigenvalues ());
 
-          retval = ovl (result.eigenvectors (), d);
+          retval = ovl (result.right_eigenvectors (), d);
+        }
+      else
+        {
+          ComplexDiagMatrix d (result.eigenvalues ());
+
+          retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
+        
         }
     }
 
diff -r 78ba6d182780 -r 05420af0a41b libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc	Thu Jun 16 22:06:03 2016 +0200
+++ b/libinterp/corefcn/xpow.cc	Tue Jun 21 00:25:51 2016 +0200
@@ -123,7 +123,7 @@ xpow (double a, const Matrix& b)
       EIG b_eig (b);
 
       ComplexColumnVector lambda (b_eig.eigenvalues ());
-      ComplexMatrix Q (b_eig.eigenvectors ());
+      ComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -174,7 +174,7 @@ xpow (double a, const ComplexMatrix& b)
   try
     {
       ComplexColumnVector lambda (b_eig.eigenvalues ());
-      ComplexMatrix Q (b_eig.eigenvectors ());
+      ComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -263,7 +263,7 @@ precision, rcond = %g", rcond);
       try
         {
           ComplexColumnVector lambda (a_eig.eigenvalues ());
-          ComplexMatrix Q (a_eig.eigenvectors ());
+          ComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
             lambda(i) = std::pow (lambda(i), b);
@@ -339,7 +339,7 @@ xpow (const Matrix& a, const Complex& b)
   try
     {
       ComplexColumnVector lambda (a_eig.eigenvalues ());
-      ComplexMatrix Q (a_eig.eigenvectors ());
+      ComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         lambda(i) = std::pow (lambda(i), b);
@@ -387,7 +387,7 @@ xpow (const Complex& a, const Matrix& b)
   try
     {
       ComplexColumnVector lambda (b_eig.eigenvalues ());
-      ComplexMatrix Q (b_eig.eigenvectors ());
+      ComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -435,7 +435,7 @@ xpow (const Complex& a, const ComplexMat
   try
     {
       ComplexColumnVector lambda (b_eig.eigenvalues ());
-      ComplexMatrix Q (b_eig.eigenvectors ());
+      ComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -524,7 +524,7 @@ precision, rcond = %g", rcond);
       try
         {
           ComplexColumnVector lambda (a_eig.eigenvalues ());
-          ComplexMatrix Q (a_eig.eigenvectors ());
+          ComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
             lambda(i) = std::pow (lambda(i), b);
@@ -559,7 +559,7 @@ xpow (const ComplexMatrix& a, const Comp
   try
     {
       ComplexColumnVector lambda (a_eig.eigenvalues ());
-      ComplexMatrix Q (a_eig.eigenvectors ());
+      ComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         lambda(i) = std::pow (lambda(i), b);
@@ -1531,7 +1531,7 @@ xpow (float a, const FloatMatrix& b)
   try
     {
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
-      FloatComplexMatrix Q (b_eig.eigenvectors ());
+      FloatComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -1583,7 +1583,7 @@ xpow (float a, const FloatComplexMatrix&
   try
     {
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
-      FloatComplexMatrix Q (b_eig.eigenvectors ());
+      FloatComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -1672,7 +1672,7 @@ precision, rcond = %g", rcond);
       try
         {
           FloatComplexColumnVector lambda (a_eig.eigenvalues ());
-          FloatComplexMatrix Q (a_eig.eigenvectors ());
+          FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
             lambda(i) = std::pow (lambda(i), b);
@@ -1738,7 +1738,7 @@ xpow (const FloatMatrix& a, const FloatC
   try
     {
       FloatComplexColumnVector lambda (a_eig.eigenvalues ());
-      FloatComplexMatrix Q (a_eig.eigenvectors ());
+      FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         lambda(i) = std::pow (lambda(i), b);
@@ -1786,7 +1786,7 @@ xpow (const FloatComplex& a, const Float
   try
     {
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
-      FloatComplexMatrix Q (b_eig.eigenvectors ());
+      FloatComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -1834,7 +1834,7 @@ xpow (const FloatComplex& a, const Float
   try
     {
       FloatComplexColumnVector lambda (b_eig.eigenvalues ());
-      FloatComplexMatrix Q (b_eig.eigenvectors ());
+      FloatComplexMatrix Q (b_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         {
@@ -1923,7 +1923,7 @@ precision, rcond = %g", rcond);
       try
         {
           FloatComplexColumnVector lambda (a_eig.eigenvalues ());
-          FloatComplexMatrix Q (a_eig.eigenvectors ());
+          FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
           for (octave_idx_type i = 0; i < nr; i++)
             lambda(i) = std::pow (lambda(i), b);
@@ -1958,7 +1958,7 @@ xpow (const FloatComplexMatrix& a, const
   try
     {
       FloatComplexColumnVector lambda (a_eig.eigenvalues ());
-      FloatComplexMatrix Q (a_eig.eigenvectors ());
+      FloatComplexMatrix Q (a_eig.right_eigenvectors ());
 
       for (octave_idx_type i = 0; i < nr; i++)
         lambda(i) = std::pow (lambda(i), b);
diff -r 78ba6d182780 -r 05420af0a41b liboctave/numeric/EIG.cc
--- a/liboctave/numeric/EIG.cc	Thu Jun 16 22:06:03 2016 +0200
+++ b/liboctave/numeric/EIG.cc	Tue Jun 21 00:25:51 2016 +0200
@@ -151,14 +151,14 @@ extern "C"
 }
 
 octave_idx_type
-EIG::init (const Matrix& a, bool calc_ev, bool balance)
+EIG::init (const Matrix& a, bool calc_rev, bool calc_lev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
       ("EIG: matrix contains Inf or NaN values");
 
   if (a.is_symmetric ())
-    return symmetric_init (a, calc_ev);
+    return symmetric_init (a, calc_rev, calc_lev);
 
   octave_idx_type n = a.rows ();
 
@@ -176,15 +176,16 @@ EIG::init (const Matrix& a, bool calc_ev
   Array<double> wi (dim_vector (n, 1));
   double *pwi = wi.fortran_vec ();
 
-  octave_idx_type tnvr = calc_ev ? n : 0;
+  octave_idx_type tnvr = calc_rev ? n : 0;
   Matrix vr (tnvr, tnvr);
   double *pvr = vr.fortran_vec ();
 
+  octave_idx_type tnvl = calc_lev ? n : 0;
+  Matrix vl (tnvl, tnvl);
+  double *pvl = vl.fortran_vec ();
+  
   octave_idx_type lwork = -1;
   double dummy_work;
-
-  double *dummy = 0;
-  octave_idx_type idummy = 1;
   
   octave_idx_type ilo;
   octave_idx_type ihi;
@@ -203,11 +204,11 @@ EIG::init (const Matrix& a, bool calc_ev
   octave_idx_type dummy_iwork;
         
   F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pwr, pwi, dummy,
-                             idummy, pvr, n, ilo, ihi, pscale,
+                             n, tmp_data, n, pwr, pwi, pvl,
+                             n, pvr, n, ilo, ihi, pscale,
                              abnrm, prconde, prcondv, &dummy_work,
                              lwork, &dummy_iwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -223,11 +224,11 @@ EIG::init (const Matrix& a, bool calc_ev
   double *pwork = work.fortran_vec ();
   
   F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pwr, pwi, dummy,
-                             idummy, pvr, n, ilo, ihi, pscale,
+                             n, tmp_data, n, pwr, pwi, pvl,
+                             n, pvr, n, ilo, ihi, pscale,
                              abnrm, prconde, prcondv, pwork,
                              lwork, &dummy_iwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -242,8 +243,10 @@ EIG::init (const Matrix& a, bool calc_ev
     (*current_liboctave_error_handler) ("dgeevx failed to converge");
 
   lambda.resize (n);
-  octave_idx_type nvr = calc_ev ? n : 0;
+  octave_idx_type nvr = calc_rev ? n : 0;
   v.resize (nvr, nvr);
+  octave_idx_type nvl = calc_lev ? n : 0;
+  w.resize (nvl, nvl);
 
   for (octave_idx_type j = 0; j < n; j++)
     {
@@ -252,6 +255,10 @@ EIG::init (const Matrix& a, bool calc_ev
           lambda.elem (j) = Complex (wr.elem (j));
           for (octave_idx_type i = 0; i < nvr; i++)
             v.elem (i, j) = vr.elem (i, j);
+          
+          for (octave_idx_type i = 0; i < nvl; i++)
+            w.elem (i, j) = vl.elem (i, j);
+          
         }
       else
         {
@@ -268,6 +275,14 @@ EIG::init (const Matrix& a, bool calc_ev
               v.elem (i, j) = Complex (real_part, imag_part);
               v.elem (i, j+1) = Complex (real_part, -imag_part);
             }
+            
+          for (octave_idx_type i = 0; i < nvl; i++)
+            {
+              double real_part = vl.elem (i, j);
+              double imag_part = vl.elem (i, j+1);
+              w.elem (i, j) = Complex (real_part, imag_part);
+              w.elem (i, j+1) = Complex (real_part, -imag_part);
+            }
           j++;
         }
     }
@@ -276,7 +291,7 @@ EIG::init (const Matrix& a, bool calc_ev
 }
 
 octave_idx_type
-EIG::symmetric_init (const Matrix& a, bool calc_ev)
+EIG::symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
 
@@ -294,7 +309,7 @@ EIG::symmetric_init (const Matrix& a, bo
   octave_idx_type lwork = -1;
   double dummy_work;
 
-  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, &dummy_work, lwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -307,7 +322,7 @@ EIG::symmetric_init (const Matrix& a, bo
   Array<double> work (dim_vector (lwork, 1));
   double *pwork = work.fortran_vec ();
 
-  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, pwork, lwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -320,20 +335,21 @@ EIG::symmetric_init (const Matrix& a, bo
     (*current_liboctave_error_handler) ("dsyev failed to converge");
 
   lambda = ComplexColumnVector (wr);
-  v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
 
   return info;
 }
 
 octave_idx_type
-EIG::init (const ComplexMatrix& a, bool calc_ev, bool balance)
+EIG::init (const ComplexMatrix& a, bool calc_rev, bool calc_lev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
       ("EIG: matrix contains Inf or NaN values");
 
   if (a.is_hermitian ())
-    return hermitian_init (a, calc_ev);
+    return hermitian_init (a, calc_rev, calc_lev);
 
   octave_idx_type n = a.rows ();
 
@@ -345,13 +361,17 @@ EIG::init (const ComplexMatrix& a, bool 
   ComplexMatrix atmp = a;
   Complex *tmp_data = atmp.fortran_vec ();
 
-  ComplexColumnVector w (n);
-  Complex *pw = w.fortran_vec ();
+  ComplexColumnVector wr (n);
+  Complex *pw = wr.fortran_vec ();
 
-  octave_idx_type nvr = calc_ev ? n : 0;
-  ComplexMatrix vtmp (nvr, nvr);
-  Complex *pv = vtmp.fortran_vec ();
+  octave_idx_type nvr = calc_rev ? n : 0;
+  ComplexMatrix vrtmp (nvr, nvr);
+  Complex *pvr = vrtmp.fortran_vec ();
 
+  octave_idx_type nvl = calc_lev ? n : 0;
+  ComplexMatrix vltmp (nvl, nvl);
+  Complex *pvl = vltmp.fortran_vec ();
+  
   octave_idx_type lwork = -1;
   Complex dummy_work;
 
@@ -359,9 +379,6 @@ EIG::init (const ComplexMatrix& a, bool 
   Array<double> rwork (dim_vector (lrwork, 1));
   double *prwork = rwork.fortran_vec ();
 
-  Complex *dummy = 0;
-  octave_idx_type idummy = 1;
-
   octave_idx_type ilo;
   octave_idx_type ihi;
   
@@ -377,11 +394,11 @@ EIG::init (const ComplexMatrix& a, bool 
   double *prcondv = rcondv.fortran_vec ();
    
   F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pw, dummy, idummy,
-                             pv, n, ilo, ihi, pscale, abnrm,
+                             n, tmp_data, n, pw, pvl, n,
+                             pvr, n, ilo, ihi, pscale, abnrm,
                              prconde, prcondv,
                              &dummy_work, lwork, prwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -397,11 +414,11 @@ EIG::init (const ComplexMatrix& a, bool 
   Complex *pwork = work.fortran_vec ();
 
   F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pw, dummy, idummy,
-                             pv, n, ilo, ihi, pscale, abnrm,
+                             n, tmp_data, n, pw, pvl, n,
+                             pvr, n, ilo, ihi, pscale, abnrm,
                              prconde, prcondv,
                              pwork, lwork, prwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -415,14 +432,15 @@ EIG::init (const ComplexMatrix& a, bool 
   if (info > 0)
     (*current_liboctave_error_handler) ("zgeevx failed to converge");
 
-  lambda = w;
-  v = vtmp;
+  lambda = wr;
+  v = vrtmp;
+  w = vltmp;
 
   return info;
 }
 
 octave_idx_type
-EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev)
+EIG::hermitian_init (const ComplexMatrix& a, bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
 
@@ -444,7 +462,7 @@ EIG::hermitian_init (const ComplexMatrix
   Array<double> rwork (dim_vector (lrwork, 1));
   double *prwork = rwork.fortran_vec ();
 
-  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, &dummy_work, lwork,
                            prwork, info
@@ -458,7 +476,7 @@ EIG::hermitian_init (const ComplexMatrix
   Array<Complex> work (dim_vector (lwork, 1));
   Complex *pwork = work.fortran_vec ();
 
-  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, pwork, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -471,13 +489,14 @@ EIG::hermitian_init (const ComplexMatrix
     (*current_liboctave_error_handler) ("zheev failed to converge");
 
   lambda = ComplexColumnVector (wr);
-  v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
-
+  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  
   return info;
 }
 
 octave_idx_type
-EIG::init (const Matrix& a, const Matrix& b, bool calc_ev)
+EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -504,7 +523,7 @@ EIG::init (const Matrix& a, const Matrix
                              F77_CHAR_ARG_LEN (1)));
 
   if (a.is_symmetric () && b.is_symmetric () && info == 0)
-    return symmetric_init (a, b, calc_ev);
+    return symmetric_init (a, b, calc_rev, calc_lev);
 
   Matrix atmp = a;
   double *atmp_data = atmp.fortran_vec ();
@@ -521,21 +540,22 @@ EIG::init (const Matrix& a, const Matrix
   Array<double> beta (dim_vector (n, 1));
   double *pbeta = beta.fortran_vec ();
 
-  octave_idx_type tnvr = calc_ev ? n : 0;
+  octave_idx_type tnvr = calc_rev ? n : 0;
   Matrix vr (tnvr, tnvr);
   double *pvr = vr.fortran_vec ();
+  
+  octave_idx_type tnvl = calc_lev ? n : 0;
+  Matrix vl (tnvl, tnvl);
+  double *pvl = vl.fortran_vec ();
 
   octave_idx_type lwork = -1;
   double dummy_work;
 
-  double *dummy = 0;
-  octave_idx_type idummy = 1;
-
   F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            &dummy_work, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -548,10 +568,10 @@ EIG::init (const Matrix& a, const Matrix
   double *pwork = work.fortran_vec ();
 
   F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            pwork, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -563,8 +583,11 @@ EIG::init (const Matrix& a, const Matrix
     (*current_liboctave_error_handler) ("dggev failed to converge");
 
   lambda.resize (n);
-  octave_idx_type nvr = calc_ev ? n : 0;
+  octave_idx_type nvr = calc_rev ? n : 0;
   v.resize (nvr, nvr);
+  
+  octave_idx_type nvl = calc_lev ? n : 0;
+  w.resize (nvl, nvl);
 
   for (octave_idx_type j = 0; j < n; j++)
     {
@@ -573,6 +596,8 @@ EIG::init (const Matrix& a, const Matrix
           lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j));
           for (octave_idx_type i = 0; i < nvr; i++)
             v.elem (i, j) = vr.elem (i, j);
+          for (octave_idx_type i = 0; i < nvl; i++)
+            w.elem (i, j) = vl.elem (i, j);
         }
       else
         {
@@ -591,6 +616,13 @@ EIG::init (const Matrix& a, const Matrix
               v.elem (i, j) = Complex (real_part, imag_part);
               v.elem (i, j+1) = Complex (real_part, -imag_part);
             }
+          for (octave_idx_type i = 0; i < nvl; i++)
+            {
+              double real_part = vl.elem (i, j);
+              double imag_part = vl.elem (i, j+1);
+              w.elem (i, j) = Complex (real_part, imag_part);
+              w.elem (i, j+1) = Complex (real_part, -imag_part);
+            }
           j++;
         }
     }
@@ -599,7 +631,7 @@ EIG::init (const Matrix& a, const Matrix
 }
 
 octave_idx_type
-EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev)
+EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
   octave_idx_type nb = b.rows ();
@@ -624,7 +656,7 @@ EIG::symmetric_init (const Matrix& a, co
   octave_idx_type lwork = -1;
   double dummy_work;
 
-  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -639,7 +671,7 @@ EIG::symmetric_init (const Matrix& a, co
   Array<double> work (dim_vector (lwork, 1));
   double *pwork = work.fortran_vec ();
 
-  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -654,13 +686,14 @@ EIG::symmetric_init (const Matrix& a, co
     (*current_liboctave_error_handler) ("dsygv failed to converge");
 
   lambda = ComplexColumnVector (wr);
-  v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
-
+  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  
   return info;
 }
 
 octave_idx_type
-EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev)
+EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, bool calc_lev)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -687,7 +720,7 @@ EIG::init (const ComplexMatrix& a, const
                              F77_CHAR_ARG_LEN (1)));
 
   if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_ev);
+    return hermitian_init (a, b, calc_rev, calc_lev);
 
   ComplexMatrix atmp = a;
   Complex *atmp_data = atmp.fortran_vec ();
@@ -701,9 +734,13 @@ EIG::init (const ComplexMatrix& a, const
   ComplexColumnVector beta (n);
   Complex *pbeta = beta.fortran_vec ();
 
-  octave_idx_type nvr = calc_ev ? n : 0;
-  ComplexMatrix vtmp (nvr, nvr);
-  Complex *pv = vtmp.fortran_vec ();
+  octave_idx_type nvr = calc_rev ? n : 0;
+  ComplexMatrix vrtmp (nvr, nvr);
+  Complex *pvr = vrtmp.fortran_vec ();
+  
+  octave_idx_type nvl = calc_lev ? n : 0;
+  ComplexMatrix vltmp (nvl, nvl);
+  Complex *pvl = vltmp.fortran_vec ();
 
   octave_idx_type lwork = -1;
   Complex dummy_work;
@@ -712,14 +749,11 @@ EIG::init (const ComplexMatrix& a, const
   Array<double> rwork (dim_vector (lrwork, 1));
   double *prwork = rwork.fortran_vec ();
 
-  Complex *dummy = 0;
-  octave_idx_type idummy = 1;
-
   F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
-                           palpha, pbeta, dummy, idummy,
-                           pv, n, &dummy_work, lwork, prwork, info
+                           palpha, pbeta, pvl, n,
+                           pvr, n, &dummy_work, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -731,10 +765,10 @@ EIG::init (const ComplexMatrix& a, const
   Complex *pwork = work.fortran_vec ();
 
   F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
-                           palpha, pbeta, dummy, idummy,
-                           pv, n, pwork, lwork, prwork, info
+                           palpha, pbeta, pvl, n,
+                           pvr, n, pwork, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -749,14 +783,15 @@ EIG::init (const ComplexMatrix& a, const
   for (octave_idx_type j = 0; j < n; j++)
     lambda.elem (j) = alpha.elem (j) / beta.elem (j);
 
-  v = vtmp;
+  v = vrtmp;
+  w = vltmp;
 
   return info;
 }
 
 octave_idx_type
 EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b,
-                     bool calc_ev)
+                     bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
   octave_idx_type nb = b.rows ();
@@ -785,7 +820,7 @@ EIG::hermitian_init (const ComplexMatrix
   Array<double> rwork (dim_vector (lrwork, 1));
   double *prwork = rwork.fortran_vec ();
 
-  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -801,7 +836,7 @@ EIG::hermitian_init (const ComplexMatrix
   Array<Complex> work (dim_vector (lwork, 1));
   Complex *pwork = work.fortran_vec ();
 
-  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -816,7 +851,8 @@ EIG::hermitian_init (const ComplexMatrix
     (*current_liboctave_error_handler) ("zhegv failed to converge");
 
   lambda = ComplexColumnVector (wr);
-  v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
-
+  v = calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  w = calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ();
+  
   return info;
 }
diff -r 78ba6d182780 -r 05420af0a41b liboctave/numeric/EIG.h
--- a/liboctave/numeric/EIG.h	Thu Jun 16 22:06:03 2016 +0200
+++ b/liboctave/numeric/EIG.h	Tue Jun 21 00:25:51 2016 +0200
@@ -42,63 +42,63 @@ public:
 
   EIG (void) : lambda (), v () { }
 
-  EIG (const Matrix& a, bool calc_eigenvectors = true,
-       bool balance = true)
+  EIG (const Matrix& a, bool calc_rev = true,
+       bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors, balance);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const Matrix& a, octave_idx_type& info,
-       bool calc_eigenvectors = true, bool balance = true)
+       bool calc_rev = true, bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors, balance);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
-  EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true)
+  EIG (const Matrix& a, const Matrix& b,
+       bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev);
   }
 
   EIG (const Matrix& a, const Matrix& b, octave_idx_type& info,
-       bool calc_eigenvectors = true)
+       bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    info = init (a, b, calc_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev);
   }
 
-  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true,
-       bool balance = true)
+  EIG (const ComplexMatrix& a, bool calc_rev = true,
+       bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors, balance);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, octave_idx_type& info,
-       bool calc_eigenvectors = true, bool balance = true)
+       bool calc_rev = true, bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors, balance);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       bool calc_eigenvectors = true)
+       bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       octave_idx_type& info, bool calc_eigenvectors = true)
+       octave_idx_type& info, bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    info = init (a, b, calc_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev);
   }
 
-  EIG (const EIG& a)
-    : lambda (a.lambda), v (a.v) { }
+  EIG (const EIG& a) : lambda (a.lambda), v (a.v), w (a.w) { }
 
   EIG& operator = (const EIG& a)
   {
@@ -106,6 +106,7 @@ public:
       {
         lambda = a.lambda;
         v = a.v;
+        w = a.w;
       }
     return *this;
   }
@@ -113,8 +114,8 @@ public:
   ~EIG (void) { }
 
   ComplexColumnVector eigenvalues (void) const { return lambda; }
-
-  ComplexMatrix eigenvectors (void) const { return v; }
+  ComplexMatrix right_eigenvectors (void) const { return v; }
+  ComplexMatrix left_eigenvectors (void) const { return w; }
 
   friend std::ostream&  operator << (std::ostream& os, const EIG& a);
 
@@ -122,30 +123,31 @@ private:
 
   ComplexColumnVector lambda;
   ComplexMatrix v;
+  ComplexMatrix w;
 
-  octave_idx_type init (const Matrix& a, bool calc_eigenvectors,
+  octave_idx_type init (const Matrix& a, bool calc_rev, bool calc_lev,
                         bool balance);
 
   octave_idx_type init (const Matrix& a, const Matrix& b,
-                        bool calc_eigenvectors);
+                        bool calc_rev, bool calc_lev);
 
-  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors,
-                        bool balance);
+  octave_idx_type init (const ComplexMatrix& a, bool calc_rev,
+                        bool calc_lev, bool balance);
 
   octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b,
-                        bool calc_eigenvectors);
+                        bool calc_rev, bool calc_lev);
 
-  octave_idx_type symmetric_init (const Matrix& a, bool calc_eigenvectors);
+  octave_idx_type symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev);
 
   octave_idx_type symmetric_init (const Matrix& a, const Matrix& b,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
 
   octave_idx_type hermitian_init (const ComplexMatrix& a,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
 
   octave_idx_type hermitian_init (const ComplexMatrix& a,
                                   const ComplexMatrix& b,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
 
 };
 
diff -r 78ba6d182780 -r 05420af0a41b liboctave/numeric/fEIG.cc
--- a/liboctave/numeric/fEIG.cc	Thu Jun 16 22:06:03 2016 +0200
+++ b/liboctave/numeric/fEIG.cc	Tue Jun 21 00:25:51 2016 +0200
@@ -146,14 +146,14 @@ extern "C"
 }
 
 octave_idx_type
-FloatEIG::init (const FloatMatrix& a, bool calc_ev, bool balance)
+FloatEIG::init (const FloatMatrix& a, bool calc_rev, bool calc_lev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
       ("EIG: matrix contains Inf or NaN values");
 
   if (a.is_symmetric ())
-    return symmetric_init (a, calc_ev);
+    return symmetric_init (a, calc_rev, calc_lev);
 
   octave_idx_type n = a.rows ();
 
@@ -171,16 +171,17 @@ FloatEIG::init (const FloatMatrix& a, bo
   Array<float> wi (dim_vector (n, 1));
   float *pwi = wi.fortran_vec ();
 
-  volatile octave_idx_type nvr = calc_ev ? n : 0;
+  volatile octave_idx_type nvr = calc_rev ? n : 0;
   FloatMatrix vr (nvr, nvr);
   float *pvr = vr.fortran_vec ();
+  
+  volatile octave_idx_type nvl = calc_lev ? n : 0;
+  FloatMatrix vl (nvl, nvl);
+  float *pvl = vl.fortran_vec ();
 
   octave_idx_type lwork = -1;
   float dummy_work;
 
-  float *dummy = 0;
-  octave_idx_type idummy = 1;
-
   octave_idx_type ilo;
   octave_idx_type ihi;
   
@@ -199,10 +200,10 @@ FloatEIG::init (const FloatMatrix& a, bo
   
   F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pwr, pwi, dummy,
-                             idummy, pvr, n,
+                             n, tmp_data, n, pwr, pwi,
+                             pvl, n, pvr, n,
                              ilo, ihi, pscale, abnrm, prconde, prcondv,
                              &dummy_work, lwork, &dummy_iwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -219,10 +220,10 @@ FloatEIG::init (const FloatMatrix& a, bo
   
   F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pwr, pwi, dummy,
-                             idummy, pvr, n,
+                             n, tmp_data, n, pwr, pwi,
+                             pvl, n, pvr, n,
                              ilo, ihi, pscale, abnrm, prconde, prcondv,
                              pwork, lwork, &dummy_iwork, info
                              F77_CHAR_ARG_LEN (1)
@@ -238,6 +239,7 @@ FloatEIG::init (const FloatMatrix& a, bo
 
   lambda.resize (n);
   v.resize (nvr, nvr);
+  w.resize (nvl, nvl);
 
   for (octave_idx_type j = 0; j < n; j++)
     {
@@ -246,6 +248,9 @@ FloatEIG::init (const FloatMatrix& a, bo
           lambda.elem (j) = FloatComplex (wr.elem (j));
           for (octave_idx_type i = 0; i < nvr; i++)
             v.elem (i, j) = vr.elem (i, j);
+
+          for (octave_idx_type i = 0; i < nvl; i++)
+            w.elem (i, j) = vl.elem (i, j);
         }
       else
         {
@@ -262,6 +267,13 @@ FloatEIG::init (const FloatMatrix& a, bo
               v.elem (i, j) = FloatComplex (real_part, imag_part);
               v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
             }
+          for (octave_idx_type i = 0; i < nvl; i++)
+            {
+              float real_part = vl.elem (i, j);
+              float imag_part = vl.elem (i, j+1);
+              w.elem (i, j) = FloatComplex (real_part, imag_part);
+              w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
+            }
           j++;
         }
     }
@@ -270,7 +282,7 @@ FloatEIG::init (const FloatMatrix& a, bo
 }
 
 octave_idx_type
-FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev)
+FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
 
@@ -288,7 +300,7 @@ FloatEIG::symmetric_init (const FloatMat
   octave_idx_type lwork = -1;
   float dummy_work;
 
-  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, &dummy_work, lwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -301,7 +313,7 @@ FloatEIG::symmetric_init (const FloatMat
   Array<float> work (dim_vector (lwork, 1));
   float *pwork = work.fortran_vec ();
 
-  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, pwork, lwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -314,20 +326,21 @@ FloatEIG::symmetric_init (const FloatMat
     (*current_liboctave_error_handler) ("ssyev failed to converge");
 
   lambda = FloatComplexColumnVector (wr);
-  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
-
+  v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  
   return info;
 }
 
 octave_idx_type
-FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev, bool balance)
+FloatEIG::init (const FloatComplexMatrix& a, bool calc_rev, bool calc_lev, bool balance)
 {
   if (a.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
       ("EIG: matrix contains Inf or NaN values");
 
   if (a.is_hermitian ())
-    return hermitian_init (a, calc_ev);
+    return hermitian_init (a, calc_rev, calc_lev);
 
   octave_idx_type n = a.rows ();
 
@@ -339,12 +352,16 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplexMatrix atmp = a;
   FloatComplex *tmp_data = atmp.fortran_vec ();
 
-  FloatComplexColumnVector w (n);
-  FloatComplex *pw = w.fortran_vec ();
+  FloatComplexColumnVector wr (n);
+  FloatComplex *pw = wr.fortran_vec ();
 
-  octave_idx_type nvr = calc_ev ? n : 0;
-  FloatComplexMatrix vtmp (nvr, nvr);
-  FloatComplex *pv = vtmp.fortran_vec ();
+  octave_idx_type nvr = calc_rev ? n : 0;
+  FloatComplexMatrix vrtmp (nvr, nvr);
+  FloatComplex *pvr = vrtmp.fortran_vec ();
+  
+  octave_idx_type nvl = calc_lev ? n : 0;
+  FloatComplexMatrix vltmp (nvl, nvl);
+  FloatComplex *pvl = vltmp.fortran_vec ();
 
   octave_idx_type lwork = -1;
   FloatComplex dummy_work;
@@ -353,9 +370,6 @@ FloatEIG::init (const FloatComplexMatrix
   Array<float> rwork (dim_vector (lrwork, 1));
   float *prwork = rwork.fortran_vec ();
 
-  FloatComplex *dummy = 0;
-  octave_idx_type idummy = 1;
-
   octave_idx_type ilo;
   octave_idx_type ihi;
   
@@ -371,11 +385,11 @@ FloatEIG::init (const FloatComplexMatrix
   float *prcondv = rcondv.fortran_vec ();
   
   F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pw, dummy, idummy,
-                             pv, n, ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             n, tmp_data, n, pw, pvl, n, pvr, n,
+                             ilo, ihi, pscale, abnrm, prconde, prcondv,
                              &dummy_work, lwork, prwork, info
                              F77_CHAR_ARG_LEN (1)
                              F77_CHAR_ARG_LEN (1)
@@ -390,11 +404,11 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplex *pwork = work.fortran_vec ();
 
   F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                              F77_CONST_CHAR_ARG2 ("N", 1),
-                             F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                             F77_CONST_CHAR_ARG2 ("N", 1),
-                             n, tmp_data, n, pw, dummy, idummy,
-                             pv, n, ilo, ihi, pscale,abnrm, prconde, prcondv,
+                             n, tmp_data, n, pw, pvl, n, pvr, n,
+                             ilo, ihi, pscale,abnrm, prconde, prcondv,
                              pwork, lwork, prwork, info
                              F77_CHAR_ARG_LEN (1)
                              F77_CHAR_ARG_LEN (1)
@@ -407,14 +421,15 @@ FloatEIG::init (const FloatComplexMatrix
   if (info > 0)
     (*current_liboctave_error_handler) ("cgeevx failed to converge");
 
-  lambda = w;
-  v = vtmp;
+  lambda = wr;
+  v = vrtmp;
+  w = vltmp;
 
   return info;
 }
 
 octave_idx_type
-FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_ev)
+FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
 
@@ -436,7 +451,7 @@ FloatEIG::hermitian_init (const FloatCom
   Array<float> rwork (dim_vector (lrwork, 1));
   float *prwork = rwork.fortran_vec ();
 
-  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, &dummy_work, lwork,
                            prwork, info
@@ -450,7 +465,7 @@ FloatEIG::hermitian_init (const FloatCom
   Array<FloatComplex> work (dim_vector (lwork, 1));
   FloatComplex *pwork = work.fortran_vec ();
 
-  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, tmp_data, n, pwr, pwork, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -463,13 +478,14 @@ FloatEIG::hermitian_init (const FloatCom
     (*current_liboctave_error_handler) ("cheev failed to converge");
 
   lambda = FloatComplexColumnVector (wr);
-  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
 
   return info;
 }
 
 octave_idx_type
-FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev)
+FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev, bool calc_lev)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -496,7 +512,7 @@ FloatEIG::init (const FloatMatrix& a, co
                              F77_CHAR_ARG_LEN (1)));
 
   if (a.is_symmetric () && b.is_symmetric () && info == 0)
-    return symmetric_init (a, b, calc_ev);
+    return symmetric_init (a, b, calc_rev, calc_lev);
 
   FloatMatrix atmp = a;
   float *atmp_data = atmp.fortran_vec ();
@@ -513,21 +529,22 @@ FloatEIG::init (const FloatMatrix& a, co
   Array<float> beta (dim_vector (n, 1));
   float *pbeta = beta.fortran_vec ();
 
-  volatile octave_idx_type nvr = calc_ev ? n : 0;
+  volatile octave_idx_type nvr = calc_rev ? n : 0;
   FloatMatrix vr (nvr, nvr);
   float *pvr = vr.fortran_vec ();
+  
+  volatile octave_idx_type nvl = calc_lev ? n : 0;
+  FloatMatrix vl (nvl, nvl);
+  float *pvl = vl.fortran_vec ();
 
   octave_idx_type lwork = -1;
   float dummy_work;
 
-  float *dummy = 0;
-  octave_idx_type idummy = 1;
-
-  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            &dummy_work, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -539,11 +556,11 @@ FloatEIG::init (const FloatMatrix& a, co
   Array<float> work (dim_vector (lwork, 1));
   float *pwork = work.fortran_vec ();
 
-  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            pwork, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -556,6 +573,8 @@ FloatEIG::init (const FloatMatrix& a, co
 
   lambda.resize (n);
   v.resize (nvr, nvr);
+  w.resize (nvl, nvl);
+  
 
   for (octave_idx_type j = 0; j < n; j++)
     {
@@ -564,6 +583,9 @@ FloatEIG::init (const FloatMatrix& a, co
           lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
           for (octave_idx_type i = 0; i < nvr; i++)
             v.elem (i, j) = vr.elem (i, j);
+          
+          for (octave_idx_type i = 0; i < nvl; i++)
+            w.elem (i, j) = vl.elem (i, j);
         }
       else
         {
@@ -582,6 +604,13 @@ FloatEIG::init (const FloatMatrix& a, co
               v.elem (i, j) = FloatComplex (real_part, imag_part);
               v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
             }
+          for (octave_idx_type i = 0; i < nvl; i++)
+            {
+              float real_part = vl.elem (i, j);
+              float imag_part = vl.elem (i, j+1);
+              w.elem (i, j) = FloatComplex (real_part, imag_part);
+              w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
+            }
           j++;
         }
     }
@@ -591,7 +620,7 @@ FloatEIG::init (const FloatMatrix& a, co
 
 octave_idx_type
 FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b,
-                          bool calc_ev)
+                          bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
   octave_idx_type nb = b.rows ();
@@ -616,7 +645,7 @@ FloatEIG::symmetric_init (const FloatMat
   octave_idx_type lwork = -1;
   float dummy_work;
 
-  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -631,7 +660,7 @@ FloatEIG::symmetric_init (const FloatMat
   Array<float> work (dim_vector (lwork, 1));
   float *pwork = work.fortran_vec ();
 
-  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -646,14 +675,15 @@ FloatEIG::symmetric_init (const FloatMat
     (*current_liboctave_error_handler) ("ssygv failed to converge");
 
   lambda = FloatComplexColumnVector (wr);
-  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
 
   return info;
 }
 
 octave_idx_type
 FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-                bool calc_ev)
+                bool calc_rev, bool calc_lev)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -680,7 +710,7 @@ FloatEIG::init (const FloatComplexMatrix
                              F77_CHAR_ARG_LEN (1)));
 
   if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_ev);
+    return hermitian_init (a, b, calc_rev, calc_lev);
 
   FloatComplexMatrix atmp = a;
   FloatComplex *atmp_data = atmp.fortran_vec ();
@@ -694,9 +724,13 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplexColumnVector beta (n);
   FloatComplex *pbeta = beta.fortran_vec ();
 
-  octave_idx_type nvr = calc_ev ? n : 0;
-  FloatComplexMatrix vtmp (nvr, nvr);
-  FloatComplex *pv = vtmp.fortran_vec ();
+  octave_idx_type nvr = calc_rev ? n : 0;
+  FloatComplexMatrix vrtmp (nvr, nvr);
+  FloatComplex *pvr = vrtmp.fortran_vec ();
+  
+  octave_idx_type nvl = calc_lev ? n : 0;
+  FloatComplexMatrix vltmp (nvl, nvl);
+  FloatComplex *pvl = vltmp.fortran_vec ();
 
   octave_idx_type lwork = -1;
   FloatComplex dummy_work;
@@ -705,14 +739,11 @@ FloatEIG::init (const FloatComplexMatrix
   Array<float> rwork (dim_vector (lrwork, 1));
   float *prwork = rwork.fortran_vec ();
 
-  FloatComplex *dummy = 0;
-  octave_idx_type idummy = 1;
-
-  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
-                           palpha, pbeta, dummy, idummy,
-                           pv, n, &dummy_work, lwork, prwork, info
+                           palpha, pbeta, pvl, n,
+                           pvr, n, &dummy_work, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -723,11 +754,11 @@ FloatEIG::init (const FloatComplexMatrix
   Array<FloatComplex> work (dim_vector (lwork, 1));
   FloatComplex *pwork = work.fortran_vec ();
 
-  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
-                           palpha, pbeta, dummy, idummy,
-                           pv, n, pwork, lwork, prwork, info
+                           palpha, pbeta, pvl, n,
+                           pvr, n, pwork, lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -742,14 +773,15 @@ FloatEIG::init (const FloatComplexMatrix
   for (octave_idx_type j = 0; j < n; j++)
     lambda.elem (j) = alpha.elem (j) / beta.elem (j);
 
-  v = vtmp;
+  v = vrtmp;
+  w = vltmp;
 
   return info;
 }
 
 octave_idx_type
-FloatEIG::hermitian_init (const FloatComplexMatrix& a,
-                          const FloatComplexMatrix& b, bool calc_ev)
+FloatEIG::hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
+                          bool calc_rev, bool calc_lev)
 {
   octave_idx_type n = a.rows ();
   octave_idx_type nb = b.rows ();
@@ -778,7 +810,7 @@ FloatEIG::hermitian_init (const FloatCom
   Array<float> rwork (dim_vector (lrwork, 1));
   float *prwork = rwork.fortran_vec ();
 
-  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -794,7 +826,7 @@ FloatEIG::hermitian_init (const FloatCom
   Array<FloatComplex> work (dim_vector (lwork, 1));
   FloatComplex *pwork = work.fortran_vec ();
 
-  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 ("U", 1),
                            n, atmp_data, n,
                            btmp_data, n,
@@ -809,7 +841,8 @@ FloatEIG::hermitian_init (const FloatCom
     (*current_liboctave_error_handler) ("zhegv failed to converge");
 
   lambda = FloatComplexColumnVector (wr);
-  v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  v = calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+  w = calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
 
   return info;
 }
diff -r 78ba6d182780 -r 05420af0a41b liboctave/numeric/fEIG.h
--- a/liboctave/numeric/fEIG.h	Thu Jun 16 22:06:03 2016 +0200
+++ b/liboctave/numeric/fEIG.h	Tue Jun 21 00:25:51 2016 +0200
@@ -40,66 +40,65 @@ FloatEIG
 
 public:
 
-  FloatEIG (void)
-    : lambda (), v () { }
+  FloatEIG (void) : lambda (), v () { }
 
-  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true,
-            bool balance = true)
+  FloatEIG (const FloatMatrix& a, bool calc_rev = true,
+            bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors, balance);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true, bool balance = true)
+            bool calc_rev = true, bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors, balance);
-  }
-
-  FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
-  {
-    init (a, b, calc_eigenvectors);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
-            octave_idx_type& info, bool calc_eigenvectors = true)
+            bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    info = init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev);
   }
 
-  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true,
-            bool balance = true)
+  FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info,
+            bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    init (a, calc_eigenvectors, balance);
+    info = init (a, b, calc_rev, calc_lev);
+  }
+
+  FloatEIG (const FloatComplexMatrix& a, bool calc_rev = true,
+            bool calc_lev = true, bool balance = true)
+    : lambda (), v ()
+  {
+    init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true, bool balance = true)
+            bool calc_rev = true, bool calc_lev = true, bool balance = true)
     : lambda (), v ()
   {
-    info = init (a, calc_eigenvectors, balance);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            bool calc_eigenvectors = true)
+            bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            octave_idx_type& info, bool calc_eigenvectors = true)
+            octave_idx_type& info, bool calc_rev = true, bool calc_lev = true)
     : lambda (), v ()
   {
-    info = init (a, b, calc_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev);
   }
 
-  FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v) { }
+  FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v), w (a.w) { }
 
   FloatEIG& operator = (const FloatEIG& a)
   {
@@ -107,6 +106,7 @@ public:
       {
         lambda = a.lambda;
         v = a.v;
+        w = a.w;
       }
     return *this;
   }
@@ -114,8 +114,8 @@ public:
   ~FloatEIG (void) { }
 
   FloatComplexColumnVector eigenvalues (void) const { return lambda; }
-
-  FloatComplexMatrix eigenvectors (void) const { return v; }
+  FloatComplexMatrix right_eigenvectors (void) const { return v; }
+  FloatComplexMatrix left_eigenvectors (void) const { return w; }
 
   friend std::ostream&  operator << (std::ostream& os, const FloatEIG& a);
 
@@ -123,24 +123,32 @@ private:
 
   FloatComplexColumnVector lambda;
   FloatComplexMatrix v;
+  FloatComplexMatrix w;
 
-  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors,
-                        bool balance);
-  octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b,
-                        bool calc_eigenvectors);
-  octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors,
+  octave_idx_type init (const FloatMatrix& a, bool calc_rev, bool calc_lev,
                         bool balance);
-  octave_idx_type init (const FloatComplexMatrix& a,
-                        const FloatComplexMatrix& b, bool calc_eigenvectors);
+                        
+  octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b,
+                        bool calc_rev, bool calc_lev);
+                        
+  octave_idx_type init (const FloatComplexMatrix& a, bool calc_rev,
+                        bool calc_lev, bool balance);
+                        
+  octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
+                        bool calc_rev, bool calc_lev);
 
-  octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_eigenvectors);
+  octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev);
+  
   octave_idx_type symmetric_init (const FloatMatrix& a, const FloatMatrix& b,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
+                                  
   octave_idx_type hermitian_init (const FloatComplexMatrix& a,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
+                                  
   octave_idx_type hermitian_init (const FloatComplexMatrix& a,
                                   const FloatComplexMatrix& b,
-                                  bool calc_eigenvectors);
+                                  bool calc_rev, bool calc_lev);
+                                  
 };
 
 #endif
 
# HG changeset patch
# User Barbara Lócsi <locsi.barbara@gmail.com>
# Date 1467407083 -7200
#      Fri Jul 01 23:04:43 2016 +0200
# Node ID 1e6628b3dc0ea6a75f1e245ebdeaa56f2082306b
# Parent  cc55436c4cbc3db2e42e6e70f83051d91417cc9e
add chol/qz flag to eig

diff -r cc55436c4cbc -r 1e6628b3dc0e libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Fri Jul 01 22:04:34 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Fri Jul 01 23:04:43 2016 +0200
@@ -52,7 +52,7 @@ The eigenvalues returned by @code{eig} a
 {
   int nargin = args.length ();
 
-  if (nargin > 2 || nargin == 0)
+  if (nargin > 3 || nargin == 0)
     print_usage ();
 
   octave_value_list retval;
@@ -75,42 +75,60 @@ The eigenvalues returned by @code{eig} a
   if (! arg_a.is_double_type () && ! arg_a.is_single_type ())
     err_wrong_type_arg ("eig", arg_a);
 
+  if (nr_a != nc_a)
+    err_square_matrix_required ("eig", "A");
+  
   bool balance = true;
-  if (nargin == 2)
+  bool force_qz = false;
+  
+  // determine if it's AEP or GEP
+  bool AEPcase = nargin == 1 || args(1).is_string ();
+  
+  if (AEPcase)
     {
-      if(args(1).is_string ())
+      if (nargin == 2)
         {
           std::string a1s = args(1).string_value ();
           if(a1s != "balance" && a1s != "nobalance")
             error ("eig: unexpected second input: %s", a1s.c_str ());
           balance = a1s == "balance";
         }
-      else
+      if (nargin == 3)
         {
-          arg_b = args(1);
-          nr_b = arg_b.rows ();
-          nc_b = arg_b.columns ();
-
-          arg_is_empty = empty_arg ("eig", nr_b, nc_b);
-          if (arg_is_empty < 0)
-            return retval;
-          else if (arg_is_empty > 0)
-            return ovl (2, Matrix ());
-
-          if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
-            err_wrong_type_arg ("eig", arg_b);
+          if(!args(2).is_string ())
+            err_wrong_type_arg ("eig", args(2));
+          std::string a2s = args(2).string_value ();
+            error ("eig: unexpected third input: %s", a2s.c_str ());
         }
     }
+  else
+    {
+      arg_b = args(1);
+      nr_b = arg_b.rows ();
+      nc_b = arg_b.columns ();
 
-  if (nr_a != nc_a)
-    err_square_matrix_required ("eig", "A");
-  
-  // determine if it's AEP or GEP
-  bool AEPcase = nargin == 1 || args(1).is_string ();
-  
+      arg_is_empty = empty_arg ("eig", nr_b, nc_b);
+      if (arg_is_empty < 0)
+        return retval;
+      else if (arg_is_empty > 0)
+        return ovl (2, Matrix ());
+
+      if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
+        err_wrong_type_arg ("eig", arg_b);
+      if(nargin == 3)
+        {
+          if(!args(2).is_string ())
+            err_wrong_type_arg ("eig", args(2));
+          std::string a2s = args(2).string_value ();
+          if(a2s != "chol" && a2s != "qz")
+            error ("eig: unexpected third input: %s", a2s.c_str ());
+          force_qz = a2s == "qz";
+        }
+    }
+    
   if (!AEPcase && nr_b != nc_b)
     err_square_matrix_required ("eig", "B");
-
+  
   Matrix tmp_a, tmp_b;
   ComplexMatrix ctmp_a, ctmp_b;
   FloatMatrix ftmp_a, ftmp_b;
@@ -141,14 +159,14 @@ The eigenvalues returned by @code{eig} a
               ftmp_a = arg_a.float_matrix_value ();
               ftmp_b = arg_b.float_matrix_value ();
 
-              result = FloatEIG (ftmp_a, ftmp_b, nargout > 1, nargout > 2);
+              result = FloatEIG (ftmp_a, ftmp_b, nargout > 1, nargout > 2, force_qz);
             }
           else
             {
               fctmp_a = arg_a.float_complex_matrix_value ();
               fctmp_b = arg_b.float_complex_matrix_value ();
 
-              result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2);
+              result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2, force_qz);
             }
         }
 
@@ -196,14 +214,14 @@ The eigenvalues returned by @code{eig} a
               tmp_a = arg_a.matrix_value ();
               tmp_b = arg_b.matrix_value ();
 
-              result = EIG (tmp_a, tmp_b, nargout > 1, nargout > 2);
+              result = EIG (tmp_a, tmp_b, nargout > 1, nargout > 2, force_qz);
             }
           else
             {
               ctmp_a = arg_a.complex_matrix_value ();
               ctmp_b = arg_b.complex_matrix_value ();
 
-              result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2);
+              result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2, force_qz);
             }
         }
 
@@ -240,11 +258,21 @@ The eigenvalues returned by @code{eig} a
 %! assert (v, [-x, x; x, x], sqrt (eps));
 
 %!test
+%! [v, d, w] = eig ([1, 2; 2, 1]);
+%! x = 1 / sqrt (2);
+%! assert (w, [-x, x; x, x], sqrt (eps));
+
+%!test
 %! [v, d] = eig ([1, 2; 2, 1], 'balance');
 %! x = 1 / sqrt (2);
 %! assert (d, [-1, 0; 0, 3], sqrt (eps));
 %! assert (v, [-x, x; x, x], sqrt (eps));
 
+%!test
+%! [v, d, w] = eig ([1, 2; 2, 1], 'balance');
+%! x = 1 / sqrt (2);
+%! assert (w, [-x, x; x, x], sqrt (eps));
+
 %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single")))
 
 %!assert (eig (single ([1, 2; 2, 1]), 'balance'), single ([-1; 3]), sqrt (eps ("single")))
@@ -256,84 +284,276 @@ The eigenvalues returned by @code{eig} a
 %! assert (v, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
+%! [v, d, w] = eig (single ([1, 2; 2, 1]));
+%! x = single (1 / sqrt (2));
+%! assert (w, [-x, x; x, x], sqrt (eps ("single")));
+
+%!test
 %! [v, d] = eig (single ([1, 2; 2, 1]), 'balance');
 %! x = single (1 / sqrt (2));
 %! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")));
 %! assert (v, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
+%! [v, d, w] = eig (single ([1, 2; 2, 1]), 'balance');
+%! x = single (1 / sqrt (2));
+%! assert (w, [-x, x; x, x], sqrt (eps ("single")));
+
+%!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'qz');
+%! [v3, d3] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (v, v3);
+%! assert (d, d2);
+%! assert (d, d3);
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'qz');
+%! [v3, d3, w3] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (v, v3);
+%! assert (d, d2);
+%! assert (d, d3);
+%! assert (w, w2);
+%! assert (w, w3);
 
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'qz');
+%! [v3, d3] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (v, v3);
+%! assert (d, d2);
+%! assert (d, d3);
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'qz');
+%! [v3, d3, w3] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (v, v3);
+%! assert (d, d2);
+%! assert (d, d3);
+%! assert (w, w2);
+%! assert (w, w3);
 
 %!test
 %! A = [1, 2; 2, 1];  B = [3, -2; -2, 3];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+
+%!test
+%! A = [1, 2; 2, 1];  B = [3, -2; -2, 3];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
 
 %!test
 %! A = single ([1, 2; 2, 1]);  B = single ([3, -2; -2, 3]);
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);  B = single ([3, -2; -2, 3]);
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
 
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
 
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
 
 %!test
 %! A = [1+3i, 2+3i; 3-8i, 8+3i];  B = [8+i, 3+i; 4-9i, 3+i];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'qz');
+%! [v3, d3] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (v, v3);
+%! assert (d, d3);
+
+%!test
+%! A = [1+3i, 2+3i; 3-8i, 8+3i];  B = [8+i, 3+i; 4-9i, 3+i];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'qz');
+%! [v3, d3, w3] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+%! assert (v, v3);
+%! assert (d, d3);
+%! assert (w, w3);
 
 %!test
 %! A = single ([1+3i, 2+3i; 3-8i, 8+3i]);  B = single ([8+i, 3+i; 4-9i, 3+i]);
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'qz');
+%! [v3, d3] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (v, v3);
+%! assert (d, d3);
+
+%!test
+%! A = single ([1+3i, 2+3i; 3-8i, 8+3i]);  B = single ([8+i, 3+i; 4-9i, 3+i]);
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'qz');
+%! [v3, d3, w3] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+%! assert (v, v3);
+%! assert (d, d3);
+%! assert (w, w3);
 
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'qz');
+%! [v3, d3] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (v, v3);
+%! assert (d, d3);
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'qz');
+%! [v3, d3, w3] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+%! assert (v, v3);
+%! assert (d, d3);
+%! assert (w, w3);
 
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
 
 %!test
 %! A = single ([1, 1+i; 1-i, 1]);  B = single ([2, 0; 0, 2]);
 %! [v, d] = eig (A, B);
+%! [v2, d2] = eig (A, B, 'chol');
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
 
-%!error eig ()
-%!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1)
-%!error <EIG requires same size matrices> eig ([1, 2; 3, 4], 2)
-%!error <must be a square matrix> eig ([1, 2; 3, 4; 5, 6])
-%!error <wrong type argument> eig ("abcd")
-%!error <eig: unexpected second input: abcd> eig ([1 2 ; 2 3], "abcd")
-%!error <wrong type argument> eig (false, [1 2 ; 2 3])
+%!test
+%! A = single ([1, 1+i; 1-i, 1]);  B = single ([2, 0; 0, 2]);
+%! [v, d, w] = eig (A, B);
+%! [v2, d2, w2] = eig (A, B, 'chol');
+%! wt = w';
+%! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
+%! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
 %!test
 %! A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1];
 %! [V1, D1] = eig (A);
@@ -342,4 +562,29 @@ The eigenvalues returned by @code{eig} a
 %! assert (V1, V2);
 %! assert (D1, D2);
 %! assert (isequal (V2, V3), false);
+
+%!test
+%! A = diag([10^-16, 10^-15]);
+%! [V1,D1] = eig(A,A);
+%! [V2,D2] = eig(A,A,'qz');
+%! [V3,D3] = eig(A,A,'chol');
+%! assert (V1, V3);
+%! assert (D1, D3);
+%! assert (isequal (V1, V2), false);
+
+%!error eig ()
+%!error eig (false)
+%!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1)
+%!error <EIG requires same size matrices> eig ([1, 2; 3, 4], 2)
+%!error <must be a square matrix> eig ([1, 2; 3, 4; 5, 6])
+%!error <wrong type argument> eig ("abcd")
+%!error <eig: unexpected second input: abcd> eig ([1 2 ; 2 3], "abcd")
+%!error <eig: unexpected second input: chol> eig ([1 2 ; 2 3], "chol")
+%!error <eig: unexpected second input: qz> eig ([1 2 ; 2 3], "qz")
+%!error <wrong type argument> eig (false, [1 2 ; 2 3])
+%!error <eig: unexpected third input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd")
+%!error <eig: unexpected third input: qz> eig ([1 2 ; 2 3], 'balance', 'qz')
+%!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], false)
+%!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3])
+
 */
\ No newline at end of file
diff -r cc55436c4cbc -r 1e6628b3dc0e liboctave/numeric/EIG.cc
--- a/liboctave/numeric/EIG.cc	Fri Jul 01 22:04:34 2016 +0200
+++ b/liboctave/numeric/EIG.cc	Fri Jul 01 23:04:43 2016 +0200
@@ -494,7 +494,7 @@ EIG::hermitian_init (const ComplexMatrix
 }
 
 octave_idx_type
-EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev)
+EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -513,14 +513,17 @@ EIG::init (const Matrix& a, const Matrix
 
   Matrix tmp = b;
   double *tmp_data = tmp.fortran_vec ();
-
-  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
-                             n, tmp_data, n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
-
-  if (a.is_symmetric () && b.is_symmetric () && info == 0)
-    return symmetric_init (a, b, calc_rev, calc_lev);
+  
+  if(!force_qz)
+    {
+      F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, tmp_data, n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
+                                 
+      if (a.is_symmetric () && b.is_symmetric () && info == 0)
+        return symmetric_init (a, b, calc_rev, calc_lev);
+    }
 
   Matrix atmp = a;
   double *atmp_data = atmp.fortran_vec ();
@@ -548,7 +551,7 @@ EIG::init (const Matrix& a, const Matrix
   octave_idx_type lwork = -1;
   double dummy_work;
 
-  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
@@ -564,7 +567,7 @@ EIG::init (const Matrix& a, const Matrix
   Array<double> work (dim_vector (lwork, 1));
   double *pwork = work.fortran_vec ();
 
-  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
                            F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
                            n, atmp_data, n, btmp_data, n,
                            par, pai, pbeta,
@@ -690,7 +693,7 @@ EIG::symmetric_init (const Matrix& a, co
 }
 
 octave_idx_type
-EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, bool calc_lev)
+EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, bool calc_lev, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -709,14 +712,17 @@ EIG::init (const ComplexMatrix& a, const
 
   ComplexMatrix tmp = b;
   Complex*tmp_data = tmp.fortran_vec ();
+  
+  if(!force_qz)
+    {
+      F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, tmp_data, n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
 
-  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
-                             n, tmp_data, n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
-
-  if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_rev, calc_lev);
+      if (a.is_hermitian () && b.is_hermitian () && info == 0)
+        return hermitian_init (a, b, calc_rev, calc_lev);
+    }
 
   ComplexMatrix atmp = a;
   Complex *atmp_data = atmp.fortran_vec ();
diff -r cc55436c4cbc -r 1e6628b3dc0e liboctave/numeric/EIG.h
--- a/liboctave/numeric/EIG.h	Fri Jul 01 22:04:34 2016 +0200
+++ b/liboctave/numeric/EIG.h	Fri Jul 01 23:04:43 2016 +0200
@@ -40,62 +40,62 @@ EIG
 
 public:
 
-  EIG (void) : lambda (), v () { }
+  EIG (void) : lambda (), v (), w () { }
 
   EIG (const Matrix& a, bool calc_rev = true,
        bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const Matrix& a, octave_idx_type& info,
        bool calc_rev = true, bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     info = init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const Matrix& a, const Matrix& b,
-       bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_rev, calc_lev);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const Matrix& a, const Matrix& b, octave_idx_type& info,
-       bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_rev, calc_lev);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const ComplexMatrix& a, bool calc_rev = true,
        bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, octave_idx_type& info,
        bool calc_rev = true, bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     info = init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_rev, calc_lev);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       octave_idx_type& info, bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+       octave_idx_type& info, bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_rev, calc_lev);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const EIG& a) : lambda (a.lambda), v (a.v), w (a.w) { }
@@ -129,13 +129,13 @@ private:
                         bool balance);
 
   octave_idx_type init (const Matrix& a, const Matrix& b,
-                        bool calc_rev, bool calc_lev);
+                        bool calc_rev, bool calc_lev, bool force_qz);
 
   octave_idx_type init (const ComplexMatrix& a, bool calc_rev,
                         bool calc_lev, bool balance);
 
   octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b,
-                        bool calc_rev, bool calc_lev);
+                        bool calc_rev, bool calc_lev, bool force_qz);
 
   octave_idx_type symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev);
 
diff -r cc55436c4cbc -r 1e6628b3dc0e liboctave/numeric/fEIG.cc
--- a/liboctave/numeric/fEIG.cc	Fri Jul 01 22:04:34 2016 +0200
+++ b/liboctave/numeric/fEIG.cc	Fri Jul 01 23:04:43 2016 +0200
@@ -483,7 +483,7 @@ FloatEIG::hermitian_init (const FloatCom
 }
 
 octave_idx_type
-FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev, bool calc_lev)
+FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev, bool calc_lev, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -502,14 +502,16 @@ FloatEIG::init (const FloatMatrix& a, co
 
   FloatMatrix tmp = b;
   float *tmp_data = tmp.fortran_vec ();
+  if(!force_qz)
+    {
+      F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, tmp_data, n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
 
-  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
-                             n, tmp_data, n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
-
-  if (a.is_symmetric () && b.is_symmetric () && info == 0)
-    return symmetric_init (a, b, calc_rev, calc_lev);
+      if (a.is_symmetric () && b.is_symmetric () && info == 0)
+        return symmetric_init (a, b, calc_rev, calc_lev);
+    }
 
   FloatMatrix atmp = a;
   float *atmp_data = atmp.fortran_vec ();
@@ -680,7 +682,7 @@ FloatEIG::symmetric_init (const FloatMat
 
 octave_idx_type
 FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-                bool calc_rev, bool calc_lev)
+                bool calc_rev, bool calc_lev, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -700,13 +702,16 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplexMatrix tmp = b;
   FloatComplex *tmp_data = tmp.fortran_vec ();
 
-  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
-                             n, tmp_data, n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
+  if(!force_qz)
+    {
+      F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, tmp_data, n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
 
-  if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_rev, calc_lev);
+      if (a.is_hermitian () && b.is_hermitian () && info == 0)
+        return hermitian_init (a, b, calc_rev, calc_lev);
+    }
 
   FloatComplexMatrix atmp = a;
   FloatComplex *atmp_data = atmp.fortran_vec ();
diff -r cc55436c4cbc -r 1e6628b3dc0e liboctave/numeric/fEIG.h
--- a/liboctave/numeric/fEIG.h	Fri Jul 01 22:04:34 2016 +0200
+++ b/liboctave/numeric/fEIG.h	Fri Jul 01 23:04:43 2016 +0200
@@ -40,62 +40,62 @@ FloatEIG
 
 public:
 
-  FloatEIG (void) : lambda (), v () { }
+  FloatEIG (void) : lambda (), v (), w () { }
 
   FloatEIG (const FloatMatrix& a, bool calc_rev = true,
             bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, octave_idx_type& info,
             bool calc_rev = true, bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
-            bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_rev, calc_lev);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info,
-            bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_rev, calc_lev);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatComplexMatrix& a, bool calc_rev = true,
             bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info,
             bool calc_rev = true, bool calc_lev = true, bool balance = true)
-    : lambda (), v ()
+    : lambda (), v (), w ()
   {
     info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_rev, calc_lev);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            octave_idx_type& info, bool calc_rev = true, bool calc_lev = true)
-    : lambda (), v ()
+            octave_idx_type& info, bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_rev, calc_lev);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatEIG& a) : lambda (a.lambda), v (a.v), w (a.w) { }
@@ -129,13 +129,13 @@ private:
                         bool balance);
                         
   octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b,
-                        bool calc_rev, bool calc_lev);
+                        bool calc_rev, bool calc_lev, bool force_qz);
                         
   octave_idx_type init (const FloatComplexMatrix& a, bool calc_rev,
                         bool calc_lev, bool balance);
                         
   octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-                        bool calc_rev, bool calc_lev);
+                        bool calc_rev, bool calc_lev, bool force_qz);
 
   octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev);
   
# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1468241214 -7200
#      Mon Jul 11 14:46:54 2016 +0200
# Node ID 60a971769d83ec535541a079628cad4ef605f12c
# Parent  1e6628b3dc0ea6a75f1e245ebdeaa56f2082306b
eig.cc: add vector/matrix flag to eig

diff -r 1e6628b3dc0e -r 60a971769d83 libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Fri Jul 01 23:04:43 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Mon Jul 11 14:46:54 2016 +0200
@@ -52,7 +52,7 @@ The eigenvalues returned by @code{eig} a
 {
   int nargin = args.length ();
 
-  if (nargin > 3 || nargin == 0)
+  if (nargin > 4 || nargin == 0)
     print_usage ();
 
   octave_value_list retval;
@@ -80,6 +80,10 @@ The eigenvalues returned by @code{eig} a
   
   bool balance = true;
   bool force_qz = false;
+  bool v_flag = false;
+  bool m_flag = false;
+  bool b_flag = false;
+  bool nb_flag = false;
   
   // determine if it's AEP or GEP
   bool AEPcase = nargin == 1 || args(1).is_string ();
@@ -89,17 +93,38 @@ The eigenvalues returned by @code{eig} a
       if (nargin == 2)
         {
           std::string a1s = args(1).string_value ();
-          if(a1s != "balance" && a1s != "nobalance")
+          if (a1s != "balance" && a1s != "nobalance" && a1s!="vector" && a1s!="matrix")
             error ("eig: unexpected second input: %s", a1s.c_str ());
-          balance = a1s == "balance";
+          b_flag = a1s == "balance";
+          nb_flag = a1s == "nobalance";
+          v_flag = a1s == "vector";
+          m_flag = a1s == "matrix";
         }
       if (nargin == 3)
         {
-          if(!args(2).is_string ())
+          std::string a1s = args(1).string_value ();
+          if (a1s != "balance" && a1s != "nobalance" && a1s!="vector" && a1s!="matrix")
+            error ("eig: unexpected second input: %s", a1s.c_str ());
+          if (!args(2).is_string ())
+            err_wrong_type_arg ("eig", args(2));
+          
+          std::string a2s = args(2).string_value ();
+          if (a2s != "balance" && a2s != "nobalance" && a2s!="vector" && a2s!="matrix")
+            error ("eig: unexpected third input: %s", a2s.c_str ());
+          b_flag = a1s == "balance" || a2s == "balance";
+          nb_flag = a1s == "nobalance" || a2s == "nobalance";
+          v_flag = a1s == "vector" || a2s == "vector";
+          m_flag = a1s == "matrix" || a2s == "matrix";
+        }
+      if (nargin == 4)
+        {
+          if (!args(2).is_string ())
             err_wrong_type_arg ("eig", args(2));
           std::string a2s = args(2).string_value ();
-            error ("eig: unexpected third input: %s", a2s.c_str ());
+          if (a2s != "balance" && a2s != "nobalance" && a2s!="vector" && a2s!="matrix")
+            error ("eig: unexpected fourth input: %s", a2s.c_str ());
         }
+
     }
   else
     {
@@ -115,16 +140,44 @@ The eigenvalues returned by @code{eig} a
 
       if (! arg_b.is_single_type () && ! arg_b.is_double_type ())
         err_wrong_type_arg ("eig", arg_b);
-      if(nargin == 3)
+      if (nargin == 3)
         {
-          if(!args(2).is_string ())
+          if (!args(2).is_string ())
             err_wrong_type_arg ("eig", args(2));
           std::string a2s = args(2).string_value ();
-          if(a2s != "chol" && a2s != "qz")
+          if (a2s != "chol" && a2s != "qz" && a2s!="vector" && a2s!="matrix")
             error ("eig: unexpected third input: %s", a2s.c_str ());
           force_qz = a2s == "qz";
+          v_flag = a2s == "vector";
+          m_flag = a2s == "matrix";
+        }
+      if (nargin == 4)
+        {
+          if (!args(2).is_string ())
+            err_wrong_type_arg ("eig", args(2));
+          std::string a2s = args(2).string_value ();
+          if (a2s != "chol" && a2s != "qz" && a2s!="vector" && a2s!="matrix")
+            error ("eig: unexpected third input: %s", a2s.c_str ());
+          
+          if (!args(3).is_string ())
+            err_wrong_type_arg ("eig", args(3));
+          std::string a3s = args(3).string_value ();
+          if (a3s != "chol" && a3s != "qz" && a3s!="vector" && a3s!="matrix")
+            error ("eig: unexpected third input: %s", a3s.c_str ());
+          
+          force_qz = a2s == "qz" || a3s == "qz";
+          v_flag = a2s == "vector" || a3s == "vector";
+          m_flag = a2s == "matrix" || a3s == "matrix";
         }
     }
+  
+  if (b_flag && nb_flag)
+    error ("eig: the combination of flags is invalid.");
+  if (v_flag && m_flag)
+    error ("eig: the combination of flags is invalid.");
+  if(nb_flag)
+    balance = false;
+  
     
   if (!AEPcase && nr_b != nc_b)
     err_square_matrix_required ("eig", "B");
@@ -172,20 +225,39 @@ The eigenvalues returned by @code{eig} a
 
       if (nargout == 0 || nargout == 1)
         {
-          retval = ovl (result.eigenvalues ());
+          if (m_flag)
+            {
+              FloatComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (d);
+            }
+          else
+            {
+              retval = ovl (result.eigenvalues ());
+            }
         }
       else if (nargout == 2)
         {
-          // Blame it on Matlab.
-          FloatComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.right_eigenvectors (), d);
+          if (v_flag)
+            {
+              retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
+            }
+          else
+            {
+              FloatComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (result.right_eigenvectors (), d);
+            }
         }
       else
         {
-          FloatComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
+          if (v_flag)
+            {
+             retval = ovl (result.right_eigenvectors (), result.eigenvalues (), result.left_eigenvectors ());
+            }
+          else
+            {
+               FloatComplexDiagMatrix d (result.eigenvalues ());
+               retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
+            }
         }
     }
   else
@@ -227,21 +299,40 @@ The eigenvalues returned by @code{eig} a
 
       if (nargout == 0 || nargout == 1)
         {
-          retval = ovl (result.eigenvalues ());
+          if (m_flag)
+            {
+              ComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (d);
+            }
+          else
+            {
+              retval = ovl (result.eigenvalues ());
+            }    
         }
       else if(nargout == 2)
         {
-          // Blame it on Matlab.
-          ComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.right_eigenvectors (), d);
+          if (v_flag)
+            {
+              retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
+            }
+          else
+            {
+              ComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (result.right_eigenvectors (), d);
+            }
         }
       else
         {
-          ComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
-        
+          
+          if (v_flag)
+            {
+             retval = ovl (result.right_eigenvectors (), result.eigenvalues (), result.left_eigenvectors ());
+            }
+          else
+            {
+               ComplexDiagMatrix d (result.eigenvalues ());
+               retval = ovl (result.right_eigenvectors (), d, result.left_eigenvectors ());
+            }
         }
     }
 
@@ -572,6 +663,439 @@ The eigenvalues returned by @code{eig} a
 %! assert (D1, D3);
 %! assert (isequal (V1, V2), false);
 
+%!test
+%! A = [1, 2; 2, 1];
+%! d = eig (A);
+%! assert (columns(d), 1);
+%! d2 = eig (A, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = [1, 2; 2, 1];
+%! [v, d, w] = eig (A, 'balance', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, 'vector', 'balance');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 2, 1];
+%! [v, d, w] = eig (A, 'balance', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, 'matrix', 'balance');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 2, 1];
+%! [v, d, w] = eig (A, 'nobalance', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, 'matrix', 'nobalance');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 2, 1];
+%! [v, d, w] = eig (A, 'nobalance', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, 'vector', 'nobalance');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! d = eig (A);
+%! assert (columns(d), 1);
+%! d2 = eig (A, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! [v, d, w] = eig (A, 'balance', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, 'vector', 'balance');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! [v, d, w] = eig (A, 'balance', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, 'matrix', 'balance');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! [v, d, w] = eig (A, 'nobalance', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, 'matrix', 'nobalance');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! [v, d, w] = eig (A, 'nobalance', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, 'vector', 'nobalance');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! d = eig (A, B);
+%! assert (columns(d), 1);
+%! d2 = eig (A, B, 'vector');
+%! assert (columns(d2), 1);
+%! [v, d3] = eig (A, B);
+%! assert (isdiag(d3), true);
+%! d4 = eig (A, B, 'matrix');
+%! assert (isdiag(d4), true);
+%! assert (d, d2);
+%! assert (d3, d4);
+%! assert (d, diag(d3));
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! [v, d, w] = eig (A, B, 'qz', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! assert (isdiag(d), true);
+%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! assert (isdiag(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
+%! [v, d, w] = eig (A, B, 'chol', 'vector');
+%! assert (columns(d), 1);
+%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
+%! assert (columns(d2), 1);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
 %!error eig ()
 %!error eig (false)
 %!error eig ([1, 2; 3, 4], [4, 3; 2, 1], 1)
@@ -584,6 +1108,12 @@ The eigenvalues returned by @code{eig} a
 %!error <wrong type argument> eig (false, [1 2 ; 2 3])
 %!error <eig: unexpected third input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd")
 %!error <eig: unexpected third input: qz> eig ([1 2 ; 2 3], 'balance', 'qz')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'balance', 'nobalance')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'nobalance', 'balance')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'matrix', 'vector')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'vector', 'matrix')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'matrix', 'vector')
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'vector', 'matrix')
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], false)
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3])

# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1468493706 -7200
#      Thu Jul 14 12:55:06 2016 +0200
# Node ID c15c6f5a9d630c126f8af06394aa600738d00671
# Parent  4cb607691bb250bff3145148a102d45504a814d9
eig.cc: fix error messages, add tests

diff -r 4cb607691bb2 -r c15c6f5a9d63 libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Mon Jul 11 14:46:54 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Thu Jul 14 12:55:06 2016 +0200
@@ -122,7 +122,10 @@ The eigenvalues returned by @code{eig} a
             err_wrong_type_arg ("eig", args(2));
           std::string a2s = args(2).string_value ();
           if (a2s != "balance" && a2s != "nobalance" && a2s!="vector" && a2s!="matrix")
-            error ("eig: unexpected fourth input: %s", a2s.c_str ());
+            error ("eig: unexpected third input: %s", a2s.c_str ());
+          std::string a3s = args(3).string_value ();
+          if (a3s != "balance" && a3s != "nobalance" && a3s!="vector" && a3s!="matrix")
+            error ("eig: unexpected fourth input: %s", a3s.c_str ());
         }
 
     }
@@ -163,7 +166,7 @@ The eigenvalues returned by @code{eig} a
             err_wrong_type_arg ("eig", args(3));
           std::string a3s = args(3).string_value ();
           if (a3s != "chol" && a3s != "qz" && a3s!="vector" && a3s!="matrix")
-            error ("eig: unexpected third input: %s", a3s.c_str ());
+            error ("eig: unexpected fourth input: %s", a3s.c_str ());
           
           force_qz = a2s == "qz" || a3s == "qz";
           v_flag = a2s == "vector" || a3s == "vector";
@@ -1108,6 +1111,8 @@ The eigenvalues returned by @code{eig} a
 %!error <wrong type argument> eig (false, [1 2 ; 2 3])
 %!error <eig: unexpected third input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd")
 %!error <eig: unexpected third input: qz> eig ([1 2 ; 2 3], 'balance', 'qz')
+%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'vector', 'abcd')
+%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], 'balance', 'matrix', 'abcd')
 %!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'balance', 'nobalance')
 %!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'nobalance', 'balance')
 %!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'matrix', 'vector')
@@ -1117,4 +1122,4 @@ The eigenvalues returned by @code{eig} a
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], false)
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3])
 
-*/
\ No newline at end of file
+*/

# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1468830958 -7200
#      Mon Jul 18 10:35:58 2016 +0200
# Node ID e4462c09aa92ed5529316d368474c617f7fb9d7a
# Parent  0545e0dacd333f91e6ac859f22c4c6fa5a390808
eig: modify code to match coding guidelines

diff -r 0545e0dacd33 -r e4462c09aa92 libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Fri Jul 15 03:16:15 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Mon Jul 18 10:35:58 2016 +0200
@@ -178,11 +178,11 @@ The eigenvalues returned by @code{eig} a
     error ("eig: the combination of flags is invalid.");
   if (v_flag && m_flag)
     error ("eig: the combination of flags is invalid.");
-  if(nb_flag)
+  if (nb_flag)
     balance = false;
   
     
-  if (!AEPcase && nr_b != nc_b)
+  if (! AEPcase && nr_b != nc_b)
     err_square_matrix_required ("eig", "B");
   
   Matrix tmp_a, tmp_b;
@@ -312,7 +312,7 @@ The eigenvalues returned by @code{eig} a
               retval = ovl (result.eigenvalues ());
             }    
         }
-      else if(nargout == 2)
+      else if (nargout == 2)
         {
           if (v_flag)
             {
@@ -357,19 +357,19 @@ The eigenvalues returned by @code{eig} a
 %! assert (w, [-x, x; x, x], sqrt (eps));
 
 %!test
-%! [v, d] = eig ([1, 2; 2, 1], 'balance');
+%! [v, d] = eig ([1, 2; 2, 1], "balance");
 %! x = 1 / sqrt (2);
 %! assert (d, [-1, 0; 0, 3], sqrt (eps));
 %! assert (v, [-x, x; x, x], sqrt (eps));
 
 %!test
-%! [v, d, w] = eig ([1, 2; 2, 1], 'balance');
+%! [v, d, w] = eig ([1, 2; 2, 1], "balance");
 %! x = 1 / sqrt (2);
 %! assert (w, [-x, x; x, x], sqrt (eps));
 
 %!assert (eig (single ([1, 2; 2, 1])), single ([-1; 3]), sqrt (eps ("single")))
 
-%!assert (eig (single ([1, 2; 2, 1]), 'balance'), single ([-1; 3]), sqrt (eps ("single")))
+%!assert (eig (single ([1, 2; 2, 1]), "balance"), single ([-1; 3]), sqrt (eps ("single")))
 
 %!test
 %! [v, d] = eig (single ([1, 2; 2, 1]));
@@ -383,21 +383,21 @@ The eigenvalues returned by @code{eig} a
 %! assert (w, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
-%! [v, d] = eig (single ([1, 2; 2, 1]), 'balance');
+%! [v, d] = eig (single ([1, 2; 2, 1]), "balance");
 %! x = single (1 / sqrt (2));
 %! assert (d, single ([-1, 0; 0, 3]), sqrt (eps ("single")));
 %! assert (v, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
-%! [v, d, w] = eig (single ([1, 2; 2, 1]), 'balance');
+%! [v, d, w] = eig (single ([1, 2; 2, 1]), "balance");
 %! x = single (1 / sqrt (2));
 %! assert (w, [-x, x; x, x], sqrt (eps ("single")));
 
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'qz');
-%! [v3, d3] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "qz");
+%! [v3, d3] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -408,8 +408,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'qz');
-%! [v3, d3, w3] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "qz");
+%! [v3, d3, w3] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -423,8 +423,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'qz');
-%! [v3, d3] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "qz");
+%! [v3, d3] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
 %! assert (v, v2);
@@ -435,8 +435,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'qz');
-%! [v3, d3, w3] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "qz");
+%! [v3, d3, w3] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
@@ -450,7 +450,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 2, 1];  B = [3, -2; -2, 3];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -459,7 +459,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 2, 1];  B = [3, -2; -2, 3];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -470,7 +470,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; 2, 1]);  B = single ([3, -2; -2, 3]);
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
 %! assert (v, v2);
@@ -479,7 +479,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; 2, 1]);  B = single ([3, -2; -2, 3]);
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
@@ -490,7 +490,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -499,7 +499,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -510,7 +510,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
 %! assert (v, v2);
@@ -519,7 +519,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
@@ -530,8 +530,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1+3i, 2+3i; 3-8i, 8+3i];  B = [8+i, 3+i; 4-9i, 3+i];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'qz');
-%! [v3, d3] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "qz");
+%! [v3, d3] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -542,8 +542,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1+3i, 2+3i; 3-8i, 8+3i];  B = [8+i, 3+i; 4-9i, 3+i];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'qz');
-%! [v3, d3, w3] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "qz");
+%! [v3, d3, w3] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -557,8 +557,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1+3i, 2+3i; 3-8i, 8+3i]);  B = single ([8+i, 3+i; 4-9i, 3+i]);
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'qz');
-%! [v3, d3] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "qz");
+%! [v3, d3] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
 %! assert (v, v2);
@@ -569,8 +569,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1+3i, 2+3i; 3-8i, 8+3i]);  B = single ([8+i, 3+i; 4-9i, 3+i]);
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'qz');
-%! [v3, d3, w3] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "qz");
+%! [v3, d3, w3] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
@@ -584,8 +584,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'qz');
-%! [v3, d3] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "qz");
+%! [v3, d3] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -596,8 +596,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'qz');
-%! [v3, d3, w3] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "qz");
+%! [v3, d3, w3] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -611,7 +611,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
 %! assert (v, v2);
@@ -620,7 +620,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps));
@@ -631,7 +631,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 1+i; 1-i, 1]);  B = single ([2, 0; 0, 2]);
 %! [v, d] = eig (A, B);
-%! [v2, d2] = eig (A, B, 'chol');
+%! [v2, d2] = eig (A, B, "chol");
 %! assert (A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps ("single")));
 %! assert (A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps ("single")));
 %! assert (v, v2);
@@ -640,7 +640,7 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 1+i; 1-i, 1]);  B = single ([2, 0; 0, 2]);
 %! [v, d, w] = eig (A, B);
-%! [v2, d2, w2] = eig (A, B, 'chol');
+%! [v2, d2, w2] = eig (A, B, "chol");
 %! wt = w';
 %! assert (wt(1, :)* A, d(1, 1) * wt(1, :) * B, sqrt (eps ("single")));
 %! assert (wt(2, :)* A, d(2, 2) * wt(2, :) * B, sqrt (eps ("single")));
@@ -651,8 +651,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [3 -2 -0.9 0; -2 4 1 -0; -0 0 -1 0; -0.5 -0.5 0.1 1];
 %! [V1, D1] = eig (A);
-%! [V2, D2] = eig (A, 'balance');
-%! [V3, D3] = eig (A, 'nobalance');
+%! [V2, D2] = eig (A, "balance");
+%! [V3, D3] = eig (A, "nobalance");
 %! assert (V1, V2);
 %! assert (D1, D2);
 %! assert (isequal (V2, V3), false);
@@ -660,8 +660,8 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = diag([10^-16, 10^-15]);
 %! [V1,D1] = eig(A,A);
-%! [V2,D2] = eig(A,A,'qz');
-%! [V3,D3] = eig(A,A,'chol');
+%! [V2,D2] = eig(A,A,"qz");
+%! [V3,D3] = eig(A,A,"chol");
 %! assert (V1, V3);
 %! assert (D1, D3);
 %! assert (isequal (V1, V2), false);
@@ -669,12 +669,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 2, 1];
 %! d = eig (A);
-%! assert (columns(d), 1);
-%! d2 = eig (A, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, 'matrix');
+%! d4 = eig (A, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -682,19 +682,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 2, 1];
-%! [v, d, w] = eig (A, 'balance', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, 'vector', 'balance');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, "balance", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "balance");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = [1, 2; 2, 1];
-%! [v, d, w] = eig (A, 'balance', 'matrix');
+%! [v, d, w] = eig (A, "balance", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, 'matrix', 'balance');
+%! [v2, d2, w2] = eig (A, "matrix", "balance");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -702,9 +702,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 2, 1];
-%! [v, d, w] = eig (A, 'nobalance', 'matrix');
+%! [v, d, w] = eig (A, "nobalance", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, 'matrix', 'nobalance');
+%! [v2, d2, w2] = eig (A, "matrix", "nobalance");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -712,10 +712,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 2, 1];
-%! [v, d, w] = eig (A, 'nobalance', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, 'vector', 'nobalance');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, "nobalance", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "nobalance");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -723,12 +723,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; 2, 1]);
 %! d = eig (A);
-%! assert (columns(d), 1);
-%! d2 = eig (A, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, 'matrix');
+%! d4 = eig (A, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -736,19 +736,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; 2, 1]);
-%! [v, d, w] = eig (A, 'balance', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, 'vector', 'balance');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, "balance", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "balance");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = single ([1, 2; 2, 1]);
-%! [v, d, w] = eig (A, 'balance', 'matrix');
+%! [v, d, w] = eig (A, "balance", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, 'matrix', 'balance');
+%! [v2, d2, w2] = eig (A, "matrix", "balance");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -756,9 +756,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; 2, 1]);
-%! [v, d, w] = eig (A, 'nobalance', 'matrix');
+%! [v, d, w] = eig (A, "nobalance", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, 'matrix', 'nobalance');
+%! [v2, d2, w2] = eig (A, "matrix", "nobalance");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -766,10 +766,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; 2, 1]);
-%! [v, d, w] = eig (A, 'nobalance', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, 'vector', 'nobalance');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, "nobalance", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "nobalance");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -777,12 +777,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -790,19 +790,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -810,9 +810,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -820,10 +820,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; -1, 1];  B = [3, 3; 1, 2];
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -831,12 +831,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -844,19 +844,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -864,9 +864,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -874,10 +874,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1, 2; -1, 1]);  B = single ([3, 3; 1, 2]);
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -885,12 +885,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -898,19 +898,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -918,9 +918,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -928,10 +928,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1+3i, 2+i; 2-i, 1+3i];  B = [5+9i, 2+i; 2-i, 5+9i];
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -939,12 +939,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -952,19 +952,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -972,9 +972,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -982,10 +982,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = single ([1+3i, 2+i; 2-i, 1+3i]);  B = single ([5+9i, 2+i; 2-i, 5+9i]);
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -993,12 +993,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -1006,10 +1006,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -1017,9 +1017,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -1027,9 +1027,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -1037,10 +1037,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 2; 3, 8];  B = [8, 3; 4, 3];
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -1048,12 +1048,12 @@ The eigenvalues returned by @code{eig} a
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
 %! d = eig (A, B);
-%! assert (columns(d), 1);
-%! d2 = eig (A, B, 'vector');
-%! assert (columns(d2), 1);
+%! assert (isvector(d), true);
+%! d2 = eig (A, B, "vector");
+%! assert (isvector(d2), true);
 %! [v, d3] = eig (A, B);
 %! assert (isdiag(d3), true);
-%! d4 = eig (A, B, 'matrix');
+%! d4 = eig (A, B, "matrix");
 %! assert (isdiag(d4), true);
 %! assert (d, d2);
 %! assert (d3, d4);
@@ -1061,19 +1061,19 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
-%! [v, d, w] = eig (A, B, 'qz', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'qz');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "qz", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "qz");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
 
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
-%! [v, d, w] = eig (A, B, 'qz', 'matrix');
+%! [v, d, w] = eig (A, B, "qz", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'qz');
+%! [v2, d2, w2] = eig (A, B, "matrix", "qz");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -1081,9 +1081,9 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
-%! [v, d, w] = eig (A, B, 'chol', 'matrix');
+%! [v, d, w] = eig (A, B, "chol", "matrix");
 %! assert (isdiag(d), true);
-%! [v2, d2, w2] = eig (A, B, 'matrix', 'chol');
+%! [v2, d2, w2] = eig (A, B, "matrix", "chol");
 %! assert (isdiag(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
@@ -1091,10 +1091,10 @@ The eigenvalues returned by @code{eig} a
 
 %!test
 %! A = [1, 1+i; 1-i, 1];  B = [2, 0; 0, 2];
-%! [v, d, w] = eig (A, B, 'chol', 'vector');
-%! assert (columns(d), 1);
-%! [v2, d2, w2] = eig (A, B, 'vector', 'chol');
-%! assert (columns(d2), 1);
+%! [v, d, w] = eig (A, B, "chol", "vector");
+%! assert (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
 %! assert (v, v2);
 %! assert (d, d2);
 %! assert (w, w2);
@@ -1110,15 +1110,15 @@ The eigenvalues returned by @code{eig} a
 %!error <eig: unexpected second input: qz> eig ([1 2 ; 2 3], "qz")
 %!error <wrong type argument> eig (false, [1 2 ; 2 3])
 %!error <eig: unexpected third input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], "abcd")
-%!error <eig: unexpected third input: qz> eig ([1 2 ; 2 3], 'balance', 'qz')
-%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'vector', 'abcd')
-%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], 'balance', 'matrix', 'abcd')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'balance', 'nobalance')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'nobalance', 'balance')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'matrix', 'vector')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], 'vector', 'matrix')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'matrix', 'vector')
-%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], 'vector', 'matrix')
+%!error <eig: unexpected third input: qz> eig ([1 2 ; 2 3], "balance", "qz")
+%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "abcd")
+%!error <eig: unexpected fourth input: abcd> eig ([1 2 ; 2 3], "balance", "matrix", "abcd")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], "balance", "nobalance")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], "nobalance", "balance")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], "matrix", "vector")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], "vector", "matrix")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], "matrix", "vector")
+%!error <eig: the combination of flags is invalid.> eig ([1 2 ; 2 3], [1 2 ; 2 3], "vector", "matrix")
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], false)
 %!error <wrong type argument> eig ([1 2 ; 2 3], [1 2 ; 2 3], [1 2 ; 2 3])
 
diff -r 0545e0dacd33 -r e4462c09aa92 liboctave/numeric/EIG.cc
--- a/liboctave/numeric/EIG.cc	Fri Jul 15 03:16:15 2016 +0200
+++ b/liboctave/numeric/EIG.cc	Mon Jul 18 10:35:58 2016 +0200
@@ -514,7 +514,7 @@ EIG::init (const Matrix& a, const Matrix
   Matrix tmp = b;
   double *tmp_data = tmp.fortran_vec ();
   
-  if(!force_qz)
+  if (! force_qz)
     {
       F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
                                  n, tmp_data, n,
@@ -713,7 +713,7 @@ EIG::init (const ComplexMatrix& a, const
   ComplexMatrix tmp = b;
   Complex*tmp_data = tmp.fortran_vec ();
   
-  if(!force_qz)
+  if (! force_qz)
     {
       F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
                                  n, tmp_data, n,
diff -r 0545e0dacd33 -r e4462c09aa92 liboctave/numeric/fEIG.cc
--- a/liboctave/numeric/fEIG.cc	Fri Jul 15 03:16:15 2016 +0200
+++ b/liboctave/numeric/fEIG.cc	Mon Jul 18 10:35:58 2016 +0200
@@ -502,7 +502,7 @@ FloatEIG::init (const FloatMatrix& a, co
 
   FloatMatrix tmp = b;
   float *tmp_data = tmp.fortran_vec ();
-  if(!force_qz)
+  if (! force_qz)
     {
       F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
                                  n, tmp_data, n,
@@ -702,7 +702,7 @@ FloatEIG::init (const FloatComplexMatrix
   FloatComplexMatrix tmp = b;
   FloatComplex *tmp_data = tmp.fortran_vec ();
 
-  if(!force_qz)
+  if (! force_qz)
     {
       F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
                                  n, tmp_data, n,
