# 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);
 
