--- a/doc/interpreter/linalg.txi	Sun Jul 17 20:40:06 2016 -0400
+++ b/doc/interpreter/linalg.txi	Mon Jul 18 07:15:07 2016 -0700
@@ -100,6 +100,8 @@
 
 @DOCSTRING(cond)
 
+@DOCSTRING(condeig)
+
 @DOCSTRING(det)
 
 @DOCSTRING(eig)

--- a/scripts/linear-algebra/cond.m	Sun Jul 17 20:40:06 2016 -0400
+++ b/scripts/linear-algebra/cond.m	Mon Jul 18 07:15:07 2016 -0700
@@ -19,7 +19,8 @@
 ## -*- texinfo -*-
 ## @deftypefn  {} {} cond (@var{A})
 ## @deftypefnx {} {} cond (@var{A}, @var{p})
-## Compute the @var{p}-norm condition number of a matrix.
+## Compute the @var{p}-norm condition number of a matrix with respect to
+## inversion.
 ##
 ## @code{cond (@var{A})} is defined as
 ## @tex
@@ -40,7 +41,7 @@
 ## indicates small changes (such as underflow or round-off error) will produce
 ## large changes in the resulting output.  In such cases the solution results
 ## from numerical computing are not likely to be accurate.
-## @seealso{condest, rcond, norm, svd}
+## @seealso{condest, rcond, condeig, norm, svd}
 ## @end deftypefn
 
 ## Author: jwe

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/condeig.m	Mon Jul 18 07:15:07 2016 -0700
@@ -0,0 +1,135 @@
+## Copyright (C) 2006, 2007 Arno Onken <asnelt@asnelt.org>
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {@var{c} =} condeig (@var{a})
+## @deftypefnx {} {[@var{v}, @var{lambda}, @var{c}] =} condeig (@var{a})
+## Compute condition numbers of a matrix with respect to eigenvalues.
+##
+## The condition numbers are the reciprocals of the cosines of the angles
+## between the left and right eigenvectors; Large values indicate that the
+## matrix has multiple distinct eigenvalues.
+##
+## The input @var{a} must be a square numeric matrix.
+##
+## The outputs are:
+##
+## @itemize @bullet
+## @item
+## @var{c} is a vector of condition numbers for the eigenvalues of
+## @var{a}.
+##
+## @item
+## @var{v} is the matrix of right eigenvectors of @var{a}.  The result is
+## equivalent to calling @code{[@var{v}, @var{lambda}] = eig (@var{a)}}.
+##
+## @item
+## @var{lambda} is the diagonal matrix of eigenvalues of @var{a}.  The
+## result is equivalent to calling
+## @code{[@var{v}, @var{lambda}] = eig (@var{a})}.
+## @end itemize
+##
+## Example
+##
+## @example
+## @group
+## a = [1, 2; 3, 4];
+## c = condeig (a)
+## @result{} [1.0150; 1.0150]
+## @end group
+## @end example
+## @seealso{eig, cond, balance}
+## @end deftypefn
+
+function [v, lambda, c] = condeig (a)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  if (! (isnumeric (a) && issquare (a)))
+    error ("condeig: A must be a square numeric matrix");
+  endif
+
+  if (issparse (a) && nargout <= 1)
+    ## Try to use svds to calculate the condition number as it will typically
+    ## be much faster than calling eig as only the smallest and largest
+    ## eigenvalue are calculated.
+
+    ## FIXME: This calculates one condition number for the entire matrix.
+    ## In the full case, separate condition numbers are calculated for every
+    ## eigenvalue.
+    try
+      s0 = svds (a, 1, 0);    # min eigenvalue
+      v = svds (a, 1) / s0;   # max eigenvalue
+    catch
+      ## Caught an error as there is a singular value exactly at zero!!
+      v = Inf;
+    end_try_catch
+    return;
+  endif
+
+  ## Right eigenvectors
+  [v, lambda] = eig (a);
+
+  if (isempty (a))
+    c = [];
+  else
+    ## Corresponding left eigenvectors
+    vl = inv (v)';
+    ## Normalize vectors
+    vl ./= repmat (sqrt (sum (abs (vl .^ 2))), rows (vl), 1);
+
+    ## Condition numbers
+    ## Definition: cos (angle) = (norm (v1) * norm (v2)) / dot (v1, v2)
+    ## Eigenvectors have been normalized so `norm (v1) * norm (v2)' = 1
+    c = abs (1 ./ dot (vl, v)');
+  endif
+
+  if (nargout <= 1)
+    v = c;
+  endif
+
+endfunction
+
+
+%!test
+%! a = [1, 2; 3, 4];
+%! c = condeig (a);
+%! expected_c = [1.0150; 1.0150];
+%! assert (c, expected_c, 0.001);
+
+%!test
+%! a = [1, 3; 5, 8];
+%! [v, lambda, c] = condeig (a);
+%! [expected_v, expected_lambda] = eig (a);
+%! expected_c = [1.0182; 1.0182];
+%! assert (v, expected_v, 0.001);
+%! assert (lambda, expected_lambda, 0.001);
+%! assert (c, expected_c, 0.001);
+
+# Test empty input
+%!assert (condeig ([]), [])
+
+## Test input validation
+%!error condeig ()
+%!error condeig (1,2)
+%!error <A must be a square numeric matrix> condeig ({1})
+%!error <A must be a square numeric matrix> condeig (ones (3,2))
+%!error <A must be a square numeric matrix> condeig (ones (2,2,2))
+

--- a/scripts/linear-algebra/module.mk	Sun Jul 17 20:40:06 2016 -0400
+++ b/scripts/linear-algebra/module.mk	Mon Jul 18 07:15:07 2016 -0700
@@ -5,6 +5,7 @@
   scripts/linear-algebra/commutation_matrix.m \
   scripts/linear-algebra/cond.m \
   scripts/linear-algebra/condest.m \
+  scripts/linear-algebra/condeig.m \
   scripts/linear-algebra/cross.m \
   scripts/linear-algebra/duplication_matrix.m \
   scripts/linear-algebra/expm.m \