Mercurial > octave

changeset 25736:138f91fb2883

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ordeig.m: Use Octave coding conventions in function. * linalg.txi: Move docstring location to be near ordschur for relevance. * ordeig.m: Use common GPL text header at start of file. Change docstring to start with an initial one sentence summary of the function. Add semicolons to all error() invocations. Eliminate (! isquasitri (B)) validation check as B has not yet been determined to even be a valid square, numeric matrix yet. Use in-place '+=' operator for efficiency. Add semicolon to assert statements that are within %!test blocks. Add tests for input validation.
author Rik <rik@octave.org>
date Fri, 03 Aug 2018 09:00:49 -0700
parents d45b0027d325
children 799dcddaf158
files doc/interpreter/linalg.txi scripts/linear-algebra/ordeig.m
diffstat 2 files changed, 47 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/doc/interpreter/linalg.txi	Fri Aug 03 17:21:00 2018 +0200
+++ b/doc/interpreter/linalg.txi	Fri Aug 03 09:00:49 2018 -0700
@@ -122,8 +122,6 @@
 
 @DOCSTRING(null)
 
-@DOCSTRING(ordeig)
-
 @DOCSTRING(orth)
 
 @DOCSTRING(mgorth)
@@ -185,6 +183,8 @@
 
 @DOCSTRING(ordschur)
 
+@DOCSTRING(ordeig)
+
 @DOCSTRING(subspace)
 
 @DOCSTRING(svd)
--- a/scripts/linear-algebra/ordeig.m	Fri Aug 03 17:21:00 2018 +0200
+++ b/scripts/linear-algebra/ordeig.m	Fri Aug 03 09:00:49 2018 -0700
@@ -3,9 +3,9 @@
 ##
 ## This file is part of Octave.
 ##
-## Octave is free software; you can redistribute it and/or modify it
+## 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
+## 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
@@ -15,20 +15,21 @@
 ##
 ## 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/>.
+## <https://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
 ## @deftypefn  {} {@var{lambda} =} ordeig (@var{A})
 ## @deftypefnx {} {@var{lambda} =} ordeig (@var{A}, @var{B})
-## Return the eigenvalues of quasi-triangular matrices in their order
-## of appearance on the matrix @var{A}.  The matrix @var{A} is usually
-## the result of a Schur factorization.  If @var{B} is given, then the
-## generalized eigenvalues of the pair @var{A}, @var{B} are computed,
-## in their order of appearance on the matrix
-## @code{@var{A}-@var{lambda}*@var{B}}. The pair @var{A}, @var{B} is
-## usually the result of a QZ decomposition.
+## Return the eigenvalues of quasi-triangular matrices in their order of
+## appearance in the matrix @var{A}.
 ##
-## @seealso{eig, schur, qz}
+## The quasi-triangular matrix @var{A} is usually the result of a Schur
+## factorization.  If called with a second input @var{B} then the generalized
+## eigenvalues of the pair @var{A}, @var{B} are returned in the order of
+## appearance of the matrix @code{@var{A}-@var{lambda}*@var{B}}.  The pair
+## @var{A}, @var{B} is usually the result of a QZ decomposition.
+##
+## @seealso{ordschur, eig, schur, qz}
 ## @end deftypefn
 
 function lambda = ordeig (A, B)
@@ -47,31 +48,31 @@
     B = eye (n);
     if (isreal (A))
       if (! isquasitri (A))
-        error ("ordeig: A must be quasi-triangular (i.e. upper block triangular with 1x1 or 2x2 blocks on the diagonal)")
+        error ("ordeig: A must be quasi-triangular (i.e., upper block triangular with 1x1 or 2x2 blocks on the diagonal)");
       endif
     else
       if (! istriu (A))
-        error ("ordeig: A must be upper-triangular when it is complex")
+        error ("ordeig: A must be upper-triangular when it is complex");
       endif
     endif
-  elseif (! isquasitri (B))
+  else
     if (! isnumeric (B) || ! issquare (B))
-      error ("ordeig: B must be a square matrix")
-    endif
-    if (length (B) != n)
-      error ("ordeig: A and B must be of the same size")
+      error ("ordeig: B must be a square matrix");
+    elseif (length (B) != n)
+      error ("ordeig: A and B must be the same size");
     endif
     if (isreal (A) && isreal (B))
       if (! isquasitri (A) || ! isquasitri (B))
-        error ("ordeig: A and B must be quasi-triangular (i.e. upper block triangular with 1x1 or 2x2 blocks on the diagonal)")
+        error ("ordeig: A and B must be quasi-triangular (i.e., upper block triangular with 1x1 or 2x2 blocks on the diagonal)");
       endif
     else
       if (! istriu (A) || ! istriu (B))
-        error ("ordeig: A and B must be both upper-triangular if any of the two is complex")
+        error ("ordeig: A and B must both be upper-triangular if either is complex");
       endif
     endif
   endif
 
+  ## Start of algorithm
   lambda = zeros (n, 1);
 
   i = 1;
@@ -85,15 +86,15 @@
           (A(i,i+1) - B(i,i+1)) * (A(i+1,i) - B(i+1,i));
       if (b > 0)
         lambda(i) = 2*c / (-b - sqrt (b^2 - 4*a*c));
-        i = i + 1;
+        i += 1;
         lambda(i) = (-b - sqrt (b^2 - 4*a*c)) / 2 / a;
       else
         lambda(i) = (-b + sqrt (b^2 - 4*a*c)) / 2 / a;
-        i = i + 1;
+        i += 1;
         lambda(i) = 2*c / (-b + sqrt (b^2 - 4*a*c));
       endif
     endif
-    i = i + 1;
+    i += 1;
   endwhile
 
 endfunction
@@ -104,20 +105,36 @@
   retval = all (tril (A, -2)(:) == 0) && all (v(1:end-1) + v(2:end) < 2);
 endfunction
 
+
 %!test
 %! A = toeplitz ([0, 1, 0, 0], [0, -1, 0, 0]);
 %! T = schur (A);
 %! lambda = ordeig (T);
-%! assert (lambda, [1.61803i; -1.61803i; 0.61803i; -0.61803i], 1e-4)
+%! assert (lambda, [1.61803i; -1.61803i; 0.61803i; -0.61803i], 1e-4);
 
 %!test
 %! A = toeplitz ([0, 1, 0, 0], [0, -1, 0, 0]);
 %! B = toeplitz ([0, 0, 0, 1], [0, -1, 0, 2]);
 %! [AA, BB] = qz (A, B);
-%! assert (isreal (AA) && isreal (BB))
+%! assert (isreal (AA) && isreal (BB));
 %! lambda = ordeig (AA, BB);
-%! assert (lambda, [0.5+0.86603i; 0.5-0.86603i; i; -i], 1e-4)
+%! assert (lambda, [0.5+0.86603i; 0.5-0.86603i; i; -i], 1e-4);
 %! [AA, BB] = qz (complex (A), complex (B));
-%! assert (iscomplex (AA) && iscomplex (BB))
+%! assert (iscomplex (AA) && iscomplex (BB));
 %! lambda = ordeig (AA, BB);
-%! assert (lambda, diag (AA) ./ diag (BB))
+%! assert (lambda, diag (AA) ./ diag (BB));
+
+## Test input validation
+%!error ordeig ()
+%!error ordeig (1,2,3)
+%!error <A must be a square matrix> ordeig ('a')
+%!error <A must be a square matrix> ordeig ([1, 2, 3])
+%!error <A must be quasi-triangular> ordeig (magic (3))
+%!error <A must be upper-triangular> ordeig ([1, 0; i, 1])
+%!error <B must be a square matrix> ordeig (1, 'a')
+%!error <B must be a square matrix> ordeig (1, [1, 2])
+%!error <A and B must be the same size> ordeig (1, ones (2,2))
+%!error <A and B must be quasi-triangular>
+%! ordeig (triu (magic (3)), magic (3))
+%!error <A and B must both be upper-triangular>
+%! ordeig ([1, 1; 0, 1], [1, 0; i, 1])