# HG changeset patch
# User Barbara Locsi <locsi.barbara@gmail.com>
# Date 1469034897 -7200
#      Wed Jul 20 19:14:57 2016 +0200
# Node ID fa85bd7065f3e86d08d655be392ba0ffcbf3a68a
# Parent  ba8a9d2934c7c20e427a72f5101514aa53a87249
Generalised eigenvalue problem

diff -r ba8a9d2934c7 -r fa85bd7065f3 libinterp/corefcn/__qp__.cc
--- a/libinterp/corefcn/__qp__.cc	Mon Jul 18 13:52:57 2016 +0200
+++ b/libinterp/corefcn/__qp__.cc	Wed Jul 20 19:14:57 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 ba8a9d2934c7 -r fa85bd7065f3 libinterp/corefcn/eig.cc
--- a/libinterp/corefcn/eig.cc	Mon Jul 18 13:52:57 2016 +0200
+++ b/libinterp/corefcn/eig.cc	Wed Jul 20 19:14:57 2016 +0200
@@ -52,7 +52,7 @@ The eigenvalues returned by @code{eig} a
 {
   int nargin = args.length ();
 
-  if (nargin > 2 || nargin == 0)
+  if (nargin > 4 || nargin == 0)
     print_usage ();
 
   octave_value_list retval;
@@ -75,7 +75,61 @@ 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 (nargin == 2)
+  if (nr_a != nc_a)
+    err_square_matrix_required ("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 ();
+  
+  if (AEPcase)
+    {
+      if (nargin == 2)
+        {
+          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 ());
+          b_flag = a1s == "balance";
+          nb_flag = a1s == "nobalance";
+          v_flag = a1s == "vector";
+          m_flag = a1s == "matrix";
+        }
+      if (nargin == 3)
+        {
+          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 ();
+          if (a2s != "balance" && a2s != "nobalance" && a2s!="vector" && a2s!="matrix")
+            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 ());
+        }
+
+    }
+  else
     {
       arg_b = args(1);
       nr_b = arg_b.rows ();
@@ -89,14 +143,48 @@ 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 (!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 ());
+          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 fourth input: %s", a3s.c_str ());
+          
+          force_qz = a2s == "qz" || a3s == "qz";
+          v_flag = a2s == "vector" || a3s == "vector";
+          m_flag = a2s == "matrix" || a3s == "matrix";
+        }
     }
-
-  if (nr_a != nc_a)
-    err_square_matrix_required ("eig", "A");
-
-  if (nargin == 2 && nr_b != nc_b)
+  
+  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");
-
+  
   Matrix tmp_a, tmp_b;
   ComplexMatrix ctmp_a, ctmp_b;
   FloatMatrix ftmp_a, ftmp_b;
@@ -105,99 +193,149 @@ The eigenvalues returned by @code{eig} a
   if (arg_a.is_single_type ())
     {
       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, nargout > 2, balance);
             }
           else
             {
               fctmp_a = arg_a.float_complex_matrix_value ();
 
-              result = FloatEIG (fctmp_a, nargout > 1);
+              result = FloatEIG (fctmp_a, nargout > 1, nargout > 2, balance);
             }
         }
-      else if (nargin == 2)
+      else
         {
           if (arg_a.is_real_type () && arg_b.is_real_type ())
             {
               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, 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);
+              result = FloatEIG (fctmp_a, fctmp_b, nargout > 1, nargout > 2, force_qz);
             }
         }
 
       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)
+        {
+          if (v_flag)
+            {
+              retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
+            }
+          else
+            {
+              FloatComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (result.right_eigenvectors (), d);
+            }
         }
       else
         {
-          // Blame it on Matlab.
-          FloatComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.eigenvectors (), d);
+          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
     {
       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, nargout > 2, balance);
             }
           else
             {
               ctmp_a = arg_a.complex_matrix_value ();
 
-              result = EIG (ctmp_a, nargout > 1);
+              result = EIG (ctmp_a, nargout > 1, nargout > 2, balance);
             }
         }
-      else if (nargin == 2)
+      else
         {
           if (arg_a.is_real_type () && arg_b.is_real_type ())
             {
               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, 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);
+              result = EIG (ctmp_a, ctmp_b, nargout > 1, nargout > 2, force_qz);
             }
         }
 
       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)
+        {
+          if (v_flag)
+            {
+              retval = ovl (result.right_eigenvectors (), result.eigenvalues ());
+            }
+          else
+            {
+              ComplexDiagMatrix d (result.eigenvalues ());
+              retval = ovl (result.right_eigenvectors (), d);
+            }
         }
       else
         {
-          // Blame it on Matlab.
-          ComplexDiagMatrix d (result.eigenvalues ());
-
-          retval = ovl (result.eigenvectors (), d);
+          
+          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 ());
+            }
         }
     }
 
@@ -213,8 +351,26 @@ 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, 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")))
+
 %!test
 %! [v, d] = eig (single ([1, 2; 2, 1]));
 %! x = single (1 / sqrt (2));
@@ -222,76 +378,748 @@ 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);
+
+%!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);
+%! [V2, D2] = eig (A, "balance");
+%! [V3, D3] = eig (A, "nobalance");
+%! 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);
+
+%!test
+%! A = [1, 2; 2, 1];
+%! d = eig (A);
+%! assert (isvector(d), true);
+%! d2 = eig (A, "vector");
+%! assert (isvector(d2), true);
+%! [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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "nobalance");
+%! assert (isvector(d2), true);
+%! assert (v, v2);
+%! assert (d, d2);
+%! assert (w, w2);
+
+%!test
+%! A = single ([1, 2; 2, 1]);
+%! d = eig (A);
+%! assert (isvector(d), true);
+%! d2 = eig (A, "vector");
+%! assert (isvector(d2), true);
+%! [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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, "vector", "nobalance");
+%! assert (isvector(d2), true);
+%! 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 (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");
+%! 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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! 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]);
+%! d = eig (A, B);
+%! 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");
+%! 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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! 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];
+%! d = eig (A, B);
+%! 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");
+%! 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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! 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]);
+%! d = eig (A, B);
+%! 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");
+%! 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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
+%! 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 (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");
+%! 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 (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; 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! 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];
+%! d = eig (A, B);
+%! 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");
+%! 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 (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");
+%! 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 (isvector(d), true);
+%! [v2, d2, w2] = eig (A, B, "vector", "chol");
+%! assert (isvector(d2), true);
+%! 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)
 %!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 <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 <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 ba8a9d2934c7 -r fa85bd7065f3 libinterp/corefcn/xpow.cc
--- a/libinterp/corefcn/xpow.cc	Mon Jul 18 13:52:57 2016 +0200
+++ b/libinterp/corefcn/xpow.cc	Wed Jul 20 19:14:57 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);
@@ -1529,7 +1529,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++)
         {
@@ -1581,7 +1581,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++)
         {
@@ -1670,7 +1670,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);
@@ -1736,7 +1736,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);
@@ -1784,7 +1784,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++)
         {
@@ -1832,7 +1832,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++)
         {
@@ -1921,7 +1921,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);
@@ -1956,7 +1956,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 ba8a9d2934c7 -r fa85bd7065f3 liboctave/numeric/EIG.cc
--- a/liboctave/numeric/EIG.cc	Mon Jul 18 13:52:57 2016 +0200
+++ b/liboctave/numeric/EIG.cc	Wed Jul 20 19:14:57 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 F77_INT&, F77_DBLE*,
-                           const F77_INT&, F77_DBLE*, F77_DBLE*,
-                           F77_DBLE*, const F77_INT&, F77_DBLE*,
-                           const F77_INT&, F77_DBLE*,
-                           const F77_INT&, F77_INT&
-                           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 F77_INT&, F77_DBLE*,
+                             const F77_INT&, F77_DBLE*, F77_DBLE*,
+                             F77_DBLE*, const F77_INT&, F77_DBLE*,
+                             const F77_INT&, F77_INT&,
+                             F77_INT&, F77_DBLE*, F77_DBLE&,
+                             F77_DBLE*, F77_DBLE*, F77_DBLE*,
+                             const F77_INT&, F77_INT*,
+                             F77_INT&
+                             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 F77_INT&, F77_DBLE_CMPLX*,
-                           const F77_INT&, F77_DBLE_CMPLX*,
-                           F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*,
-                           const F77_INT&, F77_DBLE_CMPLX*,
-                           const F77_INT&, F77_DBLE*, F77_INT&
-                           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 F77_INT&, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_DBLE_CMPLX*,
+                             F77_DBLE_CMPLX*, const F77_INT&, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_INT&,
+                             F77_INT&, F77_DBLE*, F77_DBLE&,
+                             F77_DBLE*, F77_DBLE*, F77_DBLE_CMPLX*,
+                             const F77_INT&, F77_DBLE*,
+                             F77_INT&
+                             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,
@@ -135,14 +149,14 @@ extern "C"
 }
 
 octave_idx_type
-EIG::init (const Matrix& a, bool calc_ev)
+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 ();
 
@@ -160,46 +174,77 @@ 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;
-
-  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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             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)
+                             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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             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)
+                             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;
+  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++)
     {
@@ -208,6 +253,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
         {
@@ -224,6 +273,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++;
         }
     }
@@ -232,7 +289,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 ();
 
@@ -250,7 +307,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)
@@ -263,7 +320,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)
@@ -276,20 +333,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)
+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 ();
 
@@ -301,13 +359,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;
 
@@ -315,44 +377,68 @@ 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;
-
-  F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (dummy), idummy,
-                           F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (&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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n,
+                             F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm,
+                             prconde, prcondv,
+                             F77_DBLE_CMPLX_ARG (&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, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (dummy), idummy,
-                           F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, F77_DBLE_CMPLX_ARG (tmp_data), n, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), n,
+                             F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, pscale, abnrm,
+                             prconde, prcondv,
+                             F77_DBLE_CMPLX_ARG (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;
+  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 ();
 
@@ -374,7 +460,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, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, F77_DBLE_CMPLX_ARG (&dummy_work), lwork,
                            prwork, info
@@ -388,7 +474,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, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -401,13 +487,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, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -426,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_ev);
+  
+  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 ();
@@ -450,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_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,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            &dummy_work, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -476,11 +567,11 @@ 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_CONST_CHAR_ARG2 (calc_ev ? "V" : "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,
-                           dummy, idummy, pvr, n,
+                           pvl, n, pvr, n,
                            pwork, lwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
@@ -492,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++)
     {
@@ -502,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
         {
@@ -520,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++;
         }
     }
@@ -528,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 ();
@@ -553,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,
@@ -568,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,
@@ -583,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, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -608,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, F77_DBLE_CMPLX_ARG (tmp_data), n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
 
-  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
-                             n, F77_DBLE_CMPLX_ARG (tmp_data), n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
-
-  if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_ev);
+      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 ();
@@ -629,9 +736,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;
@@ -640,14 +751,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, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n,
-                           F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (dummy), idummy,
-                           F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info
+                           F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (pvl), n,
+                           F77_DBLE_CMPLX_ARG (pvr), n, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -659,10 +767,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),
-                           n, F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n,
-                           F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (dummy), idummy,
-                           F77_DBLE_CMPLX_ARG (pv), n, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
+                           F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                           n,  F77_DBLE_CMPLX_ARG (atmp_data), n, F77_DBLE_CMPLX_ARG (btmp_data), n,
+                           F77_DBLE_CMPLX_ARG (palpha), F77_DBLE_CMPLX_ARG (pbeta), F77_DBLE_CMPLX_ARG (pvl), n,
+                           F77_DBLE_CMPLX_ARG (pvr), n, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -677,14 +785,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 ();
@@ -713,7 +822,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, F77_DBLE_CMPLX_ARG (atmp_data), n,
                            F77_DBLE_CMPLX_ARG (btmp_data), n,
@@ -729,7 +838,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, F77_DBLE_CMPLX_ARG (atmp_data), n,
                            F77_DBLE_CMPLX_ARG (btmp_data), n,
@@ -744,7 +853,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 ba8a9d2934c7 -r fa85bd7065f3 liboctave/numeric/EIG.h
--- a/liboctave/numeric/EIG.h	Mon Jul 18 13:52:57 2016 +0200
+++ b/liboctave/numeric/EIG.h	Wed Jul 20 19:14:57 2016 +0200
@@ -40,62 +40,65 @@ EIG
 
 public:
 
-  EIG (void) : lambda (), v () { }
+  EIG (void) : lambda (), v (), w () { }
 
-  EIG (const Matrix& a, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  EIG (const Matrix& a, bool calc_rev = true,
+       bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_rev, calc_lev, balance);
   }
 
-  EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  EIG (const Matrix& a, octave_idx_type& info,
+       bool calc_rev = true, bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
-  EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  EIG (const Matrix& a, const Matrix& b,
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const Matrix& a, const Matrix& b, octave_idx_type& info,
-       bool calc_eigenvectors = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
-  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  EIG (const ComplexMatrix& a, bool calc_rev = true,
+       bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, octave_idx_type& info,
-       bool calc_eigenvectors = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       bool calc_eigenvectors = true)
-    : lambda (), v ()
+       bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   EIG (const ComplexMatrix& a, const ComplexMatrix& b,
-       octave_idx_type& info, bool calc_eigenvectors = 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_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
-  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)
   {
@@ -103,6 +106,7 @@ public:
       {
         lambda = a.lambda;
         v = a.v;
+        w = a.w;
       }
     return *this;
   }
@@ -110,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);
 
@@ -119,28 +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, bool force_qz);
 
-  octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors);
+  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, bool force_qz);
 
-  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 ba8a9d2934c7 -r fa85bd7065f3 liboctave/numeric/fEIG.cc
--- a/liboctave/numeric/fEIG.cc	Mon Jul 18 13:52:57 2016 +0200
+++ b/liboctave/numeric/fEIG.cc	Wed Jul 20 19:14:57 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 F77_INT&, F77_REAL*,
-                           const F77_INT&, F77_REAL*, F77_REAL*, F77_REAL*,
-                           const F77_INT&, F77_REAL*,
-                           const F77_INT&, F77_REAL*,
-                           const F77_INT&, F77_INT&
-                           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 F77_INT&, F77_REAL*,
+                             const F77_INT&, F77_REAL*, F77_REAL*, F77_REAL*,
+                             const F77_INT&, F77_REAL*,
+                             const F77_INT&, F77_INT&,
+                             F77_INT&, F77_REAL*, F77_REAL&, F77_REAL*,
+                             F77_REAL*, F77_REAL*, const F77_INT&,
+                             F77_INT*, F77_INT&
+                             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 F77_INT&, F77_CMPLX*,
-                           const F77_INT&, F77_CMPLX*, F77_CMPLX*,
-                           const F77_INT&, F77_CMPLX*,
-                           const F77_INT&, F77_CMPLX*,
-                           const F77_INT&, F77_REAL*, F77_INT&
-                           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 F77_INT&, F77_CMPLX*,
+                             const F77_INT&, F77_CMPLX*, F77_CMPLX*,
+                             const F77_INT&, F77_CMPLX*,
+                             const F77_INT&, F77_INT&,
+                             F77_INT&, F77_REAL*, F77_REAL&, F77_REAL*,
+                             F77_REAL*, F77_CMPLX*, const F77_INT&,
+                             F77_REAL*, F77_INT&
+                             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,
@@ -132,14 +144,14 @@ extern "C"
 }
 
 octave_idx_type
-FloatEIG::init (const FloatMatrix& a, bool calc_ev)
+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 ();
 
@@ -157,45 +169,75 @@ 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;
-
-  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_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             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)
+                             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_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             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)
+                             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);
+  w.resize (nvl, nvl);
 
   for (octave_idx_type j = 0; j < n; j++)
     {
@@ -204,6 +246,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
         {
@@ -220,6 +265,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++;
         }
     }
@@ -228,7 +280,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 ();
 
@@ -246,7 +298,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)
@@ -259,7 +311,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)
@@ -272,20 +324,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)
+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 ();
 
@@ -297,12 +350,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;
@@ -311,44 +368,66 @@ 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 (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
-                           n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (dummy), idummy,
-                           F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (&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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
+                             ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             F77_CMPLX_ARG (&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, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (dummy), idummy,
-                           F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (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 (calc_lev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
+                             F77_CONST_CHAR_ARG2 ("N", 1),
+                             n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
+                             ilo, ihi, pscale, abnrm, prconde, prcondv,
+                             F77_CMPLX_ARG (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;
+  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 ();
 
@@ -370,7 +449,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, F77_CMPLX_ARG (tmp_data), n, pwr, F77_CMPLX_ARG (&dummy_work), lwork,
                            prwork, info
@@ -384,7 +463,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, F77_CMPLX_ARG (tmp_data), n, pwr, F77_CMPLX_ARG (pwork), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
@@ -397,13 +476,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, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -422,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_ev);
+      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 ();
@@ -446,21 +528,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)));
@@ -472,11 +555,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)));
@@ -489,6 +572,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++)
     {
@@ -497,6 +582,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
         {
@@ -515,6 +603,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++;
         }
     }
@@ -524,7 +619,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 ();
@@ -549,7 +644,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,
@@ -564,7 +659,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,
@@ -579,14 +674,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, bool force_qz)
 {
   if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
     (*current_liboctave_error_handler)
@@ -606,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, F77_CMPLX_ARG (tmp_data), n,
-                             info
-                             F77_CHAR_ARG_LEN (1)));
+  if (! force_qz)
+    {
+      F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+                                 n, F77_CMPLX_ARG (tmp_data), n,
+                                 info
+                                 F77_CHAR_ARG_LEN (1)));
 
-  if (a.is_hermitian () && b.is_hermitian () && info == 0)
-    return hermitian_init (a, b, calc_ev);
+      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 ();
@@ -626,9 +725,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;
@@ -637,14 +740,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, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n,
-                           F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (dummy), idummy,
-                           F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (&dummy_work), lwork, prwork, info
+                           F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (pvl), n,
+                           F77_CMPLX_ARG (pvr), n, F77_CMPLX_ARG (&dummy_work), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -655,11 +755,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, F77_CMPLX_ARG (atmp_data), n, F77_CMPLX_ARG (btmp_data), n,
-                           F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (dummy), idummy,
-                           F77_CMPLX_ARG (pv), n, F77_CMPLX_ARG (pwork), lwork, prwork, info
+                           F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta), F77_CMPLX_ARG (pvl), n,
+                           F77_CMPLX_ARG (pvr), n, F77_CMPLX_ARG (pwork), lwork, prwork, info
                            F77_CHAR_ARG_LEN (1)
                            F77_CHAR_ARG_LEN (1)));
 
@@ -674,14 +774,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 ();
@@ -710,7 +811,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, F77_CMPLX_ARG (atmp_data), n,
                            F77_CMPLX_ARG (btmp_data), n,
@@ -726,7 +827,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, F77_CMPLX_ARG (atmp_data), n,
                            F77_CMPLX_ARG (btmp_data), n,
@@ -741,7 +842,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 ba8a9d2934c7 -r fa85bd7065f3 liboctave/numeric/fEIG.h
--- a/liboctave/numeric/fEIG.h	Mon Jul 18 13:52:57 2016 +0200
+++ b/liboctave/numeric/fEIG.h	Wed Jul 20 19:14:57 2016 +0200
@@ -40,64 +40,65 @@ FloatEIG
 
 public:
 
-  FloatEIG (void)
-    : lambda (), v () { }
+  FloatEIG (void) : lambda (), v (), w () { }
 
-  FloatEIG (const FloatMatrix& a, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  FloatEIG (const FloatMatrix& a, bool calc_rev = true,
+            bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    info = init (a, b, calc_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
-  FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true)
-    : lambda (), v ()
+  FloatEIG (const FloatComplexMatrix& a, bool calc_rev = true,
+            bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    init (a, calc_eigenvectors);
+    init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool balance = true)
+    : lambda (), v (), w ()
   {
-    info = init (a, calc_eigenvectors);
+    info = init (a, calc_rev, calc_lev, balance);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            bool calc_eigenvectors = true)
-    : lambda (), v ()
+            bool calc_rev = true, bool calc_lev = true, bool force_qz = false)
+    : lambda (), v (), w ()
   {
-    init (a, b, calc_eigenvectors);
+    init (a, b, calc_rev, calc_lev, force_qz);
   }
 
   FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
-            octave_idx_type& info, bool calc_eigenvectors = 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_eigenvectors);
+    info = init (a, b, calc_rev, calc_lev, force_qz);
   }
 
-  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)
   {
@@ -105,6 +106,7 @@ public:
       {
         lambda = a.lambda;
         v = a.v;
+        w = a.w;
       }
     return *this;
   }
@@ -112,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);
 
@@ -121,22 +123,32 @@ private:
 
   FloatComplexColumnVector lambda;
   FloatComplexMatrix v;
+  FloatComplexMatrix w;
 
-  octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const FloatMatrix& a, bool calc_rev, bool calc_lev,
+                        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,
-                        const FloatComplexMatrix& b, bool calc_eigenvectors);
+                        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 force_qz);
 
-  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
