Mercurial > octave
changeset 22612:c05377052b50
maint: Merge stable to default.
| author | Rik <rik@octave.org> |
|---|---|
| date | Tue, 11 Oct 2016 08:56:03 -0700 |
| parents | dc872d5d74c5 (diff) 95aff68c443d (current diff) |
| children | edd04ce99891 |
| files | configure.ac |
| diffstat | 15 files changed, 1204 insertions(+), 130 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Mon Oct 10 17:34:38 2016 -0400 +++ b/NEWS Tue Oct 11 08:56:03 2016 -0700 @@ -1,3 +1,11 @@ +Summary of important user-visible changes for version 4.4: +--------------------------------------------------------- + + ** Other new functions added in 4.4: + + gsvd + + Summary of important user-visible changes for version 4.2: ---------------------------------------------------------
--- a/configure.ac Mon Oct 10 17:34:38 2016 -0400 +++ b/configure.ac Tue Oct 11 08:56:03 2016 -0700 @@ -19,14 +19,14 @@ ### <http://www.gnu.org/licenses/>. AC_PREREQ([2.63]) -AC_INIT([GNU Octave], [4.2.0-rc2], [http://octave.org/bugs.html], [octave]) +AC_INIT([GNU Octave], [4.3.0+], [http://octave.org/bugs.html], [octave]) dnl Note that the version number is duplicated here and in AC_INIT dnl because AC_INIT requires it to be static, not computed from dnl shell variables. OCTAVE_MAJOR_VERSION=4 -OCTAVE_MINOR_VERSION=2 -OCTAVE_PATCH_VERSION=0-rc2 +OCTAVE_MINOR_VERSION=3 +OCTAVE_PATCH_VERSION=0+ dnl PACKAGE_VERSION is set by the AC_INIT VERSION arg OCTAVE_VERSION="$PACKAGE_VERSION"
--- a/doc/interpreter/linalg.txi Mon Oct 10 17:34:38 2016 -0400 +++ b/doc/interpreter/linalg.txi Tue Oct 11 08:56:03 2016 -0700 @@ -108,6 +108,8 @@ @DOCSTRING(givens) +@DOCSTRING(gsvd) + @DOCSTRING(planerot) @DOCSTRING(inv)
--- a/libinterp/corefcn/data.cc Mon Oct 10 17:34:38 2016 -0400 +++ b/libinterp/corefcn/data.cc Tue Oct 11 08:56:03 2016 -0700 @@ -3549,6 +3549,49 @@ return retval; } +/* +%!error <undefined> 1+Infj +%!error <undefined> 1+Infi + +%!test <31974> +%! assert (Inf + Inf*i, complex (Inf, Inf)) +%! +%! assert (1 + Inf*i, complex (1, Inf)) +%! assert (1 + Inf*j, complex (1, Inf)) +%! +%! ## whitespace should not affect parsing +%! assert (1+Inf*i, complex (1, Inf)) +%! assert (1+Inf*j, complex (1, Inf)) +%! +%! assert (NaN*j, complex (0, NaN)) +%! +%! assert (Inf * 4j, complex (0, Inf)) + +%!test <31974> +%! x = Inf; +%! assert (x * j, complex (0, Inf)) +%! j = complex (0, 1); +%! assert (Inf * j, complex (0, Inf)) + +%!test <31974> +%! exp = complex (zeros (2, 2), Inf (2, 2)); +%! assert (Inf (2, 2) * j, exp) +%! assert (Inf (2, 2) .* j, exp) +%! assert (Inf * (ones (2, 2) * j), exp) +%! assert (Inf (2, 2) .* (ones (2, 2) * j), exp) + +%!test <31974> +%! assert ([Inf; 0] * [i, 0], complex ([NaN NaN; 0 0], [Inf NaN; 0 0])) +%! assert ([Inf, 0] * [i; 0], complex (NaN, Inf)) +%! assert ([Inf, 0] .* [i, 0], complex ([0 0], [Inf 0])) + +%!test <31974> +%! m = @(x, y) x * y; +%! d = @(x, y) x / y; +%! assert (m (Inf, i), complex (0, +Inf)) +%! assert (d (Inf, i), complex (0, -Inf)) +*/ + DEFUN (isreal, args, , doc: /* -*- texinfo -*- @deftypefn {} {} isreal (@var{x})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/gsvd.cc Tue Oct 11 08:56:03 2016 -0700 @@ -0,0 +1,402 @@ +// Copyright (C) 2016 Barbara Lócsi +// Copyright (C) 2006, 2010 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 1996, 1997 John W. Eaton +// +// This program 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. +// +// This program 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 +// this program; if not, see <http://www.gnu.org/licenses/>. + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include "dMatrix.h" +#include "CMatrix.h" +#include "dDiagMatrix.h" +#include "gsvd.h" + +#include "defun.h" +#include "defun-int.h" +#include "error.h" +#include "errwarn.h" +#include "utils.h" +#include "ovl.h" +#include "ov.h" + + +template <typename T> +static typename octave::math::gsvd<T>::Type +gsvd_type (int nargout) +{ + return ((nargout == 0 || nargout == 1) + ? octave::math::gsvd<T>::Type::sigma_only + : (nargout > 5) ? octave::math::gsvd<T>::Type::std + : octave::math::gsvd<T>::Type::economy); +} + +// Named like this to avoid conflicts with the gsvd class. +template <typename T> +static octave_value_list +function_gsvd (const T& A, const T& B, const octave_idx_type nargout) +{ + octave::math::gsvd<T> result (A, B, gsvd_type<T> (nargout)); + + octave_value_list retval (nargout); + if (nargout < 2) + { + DiagMatrix sigA = result.singular_values_A (); + DiagMatrix sigB = result.singular_values_B (); + for (int i = sigA.rows () - 1; i >= 0; i--) + sigA.dgxelem(i) /= sigB.dgxelem(i); + retval(0) = sigA.diag (); + } + else + { + retval(0) = result.left_singular_matrix_A (); + retval(1) = result.left_singular_matrix_B (); + if (nargout > 2) + retval(2) = result.right_singular_matrix (); + if (nargout > 3) + retval(3) = result.singular_values_A (); + if (nargout > 4) + retval(4) = result.singular_values_B (); + if (nargout > 5) + retval(5) = result.R_matrix (); + } + return retval; +} + +DEFUN (gsvd, args, nargout, + doc: /* -*- texinfo -*- +@deftypefn {} {@var{S} =} gsvd (@var{A}, @var{B}) +@deftypefnx {} {[@var{U}, @var{V}, @var{X}, @var{C}, @var{S}] =} gsvd (@var{A}, @var{B}) +@deftypefnx {} {[@var{U}, @var{V}, @var{X}, @var{C}, @var{S}] =} gsvd (@var{A}, @var{B}, 0) +Compute the generalized singular value decomposition of (@var{A}, @var{B}): + +@tex +$$ A = U C X^\dagger $$ +$$ B = V S X^\dagger $$ +$$ C^\dagger C + S^\dagger S = eye (columns (A)) $$ +@end tex +@ifnottex + +@example +@group +A = U*C*X' +B = V*S*X' +C'*C + S'*S = eye (columns (A)) +@end group +@end example + +@end ifnottex + +The function @code{gsvd} normally returns just the vector of generalized singular values +@tex +$$ \sqrt{{{diag (C^\dagger C)} \over {diag (S^\dagger S)}}} $$ +@end tex +@ifnottex +@code{sqrt (diag (C'*C) ./ diag (S'*S))}. +@end ifnottex +If asked for five return values, it also computes +@tex +$U$, $V$, $X$, and $C$. +@end tex +@ifnottex +U, V, X, and C. +@end ifnottex + +If the optional third input is present, @code{gsvd} constructs the +"economy-sized" decomposition where the number of columns of @var{U}, @var{V} +and the number of rows of @var{C}, @var{S} is less than or equal to the number +of columns of @var{A}. This option is not yet implemented. + +Programming Note: the code is a wrapper to the corresponding @sc{lapack} dggsvd +and zggsvd routines. + +@seealso{svd} +@end deftypefn */) +{ + int nargin = args.length (); + + if (nargin < 2 || nargin > 3) + print_usage (); + else if (nargin == 3) + warning ("gsvd: economy-sized decomposition is not yet implemented, returning full decomposition"); + + octave_value_list retval; + + octave_value argA = args(0); + octave_value argB = args(1); + + octave_idx_type nr = argA.rows (); + octave_idx_type nc = argA.columns (); + + octave_idx_type np = argB.columns (); + + // FIXME: This "special" case should be handled in the gsvd class, not here + if (nr == 0 || nc == 0) + { + retval = octave_value_list (nargout); + if (nargout < 2) // S = gsvd (A, B) + retval(0) = Matrix (0, 1); + else // [U, V, X, C, S, R] = gsvd (A, B) + { + retval(0) = identity_matrix (nc, nc); + retval(1) = identity_matrix (nc, nc); + if (nargout > 2) + retval(2) = identity_matrix (nr, nr); + if (nargout > 3) + retval(3) = Matrix (nr, nc); + if (nargout > 4) + retval(4) = identity_matrix (nr, nr); + if (nargout > 5) + retval(5) = identity_matrix (nr, nr); + } + } + else + { + if (nc != np) + print_usage (); + + // FIXME: Remove when interface to gsvd single class has been written + if (argA.is_single_type () && argB.is_single_type ()) + warning ("gsvd: no implementation for single matrices, converting to double"); + + if (argA.is_real_type () && argB.is_real_type ()) + { + Matrix tmpA = argA.xmatrix_value ("gsvd: A must be a real or complex matrix"); + Matrix tmpB = argB.xmatrix_value ("gsvd: B must be a real or complex matrix"); + + if (tmpA.any_element_is_inf_or_nan ()) + error ("gsvd: A cannot have Inf or NaN values"); + if (tmpB.any_element_is_inf_or_nan ()) + error ("gsvd: B cannot have Inf or NaN values"); + + retval = function_gsvd (tmpA, tmpB, nargout); + } + else if (argA.is_complex_type () || argB.is_complex_type ()) + { + ComplexMatrix ctmpA = argA.xcomplex_matrix_value ("gsvd: A must be a real or complex matrix"); + ComplexMatrix ctmpB = argB.xcomplex_matrix_value ("gsvd: B must be a real or complex matrix"); + + if (ctmpA.any_element_is_inf_or_nan ()) + error ("gsvd: A cannot have Inf or NaN values"); + if (ctmpB.any_element_is_inf_or_nan ()) + error ("gsvd: B cannot have Inf or NaN values"); + + retval = function_gsvd (ctmpA, ctmpB, nargout); + } + else + error ("gsvd: A and B must be real or complex matrices"); + } + + return retval; +} + +/* + +## Basic test of decomposition +%!test <48807> +%! A = reshape (1:15,5,3); +%! B = magic (3); +%! [U,V,X,C,S] = gsvd (A,B); +%! assert (U*C*X', A, 50*eps); +%! assert (V*S*X', B, 50*eps); +%! S0 = gsvd (A, B); +%! S1 = svd (A / B); +%! assert (S0, S1, 10*eps); + +## a few tests for gsvd.m +%!shared A, A0, B, B0, U, V, C, S, X, R, D1, D2 +%! A0 = randn (5, 3); +%! B0 = diag ([1 2 4]); +%! A = A0; +%! B = B0; + +## A (5x3) and B (3x3) are full rank +%!test <48807> +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros (5, 3); D1(1:3, 1:3) = C; +%! D2 = S; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 5x3 full rank, B: 3x3 rank deficient +%!test <48807> +%! B(2, 2) = 0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros (5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; +%! D2 = [zeros(2, 1) S; zeros(1, 3)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 5x3 rank deficient, B: 3x3 full rank +%!test <48807> +%! B = B0; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 3); D1(1:3, 1:3) = C; +%! D2 = S; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A and B are both rank deficient +%!test <48807> +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 2); D1(1:2, 1:2) = C; +%! D2 = [S; zeros(1, 2)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6); + +## A (now 3x5) and B (now 5x5) are full rank +%!test <48807> +%! A = A0.'; +%! B0 = diag ([1 2 4 8 16]); +%! B = B0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = [C zeros(3,2)]; +%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 3x5 full rank, B: 5x5 rank deficient +%!test <48807> +%! B(2, 2) = 0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; +%! D2 = zeros(5, 5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 3x5 rank deficient, B: 5x5 full rank +%!test <48807> +%! B = B0; +%! A(3, :) = 2*A(1, :) - A(2, :); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros (3, 5); D1(1:3, 1:3) = C; +%! D2 = zeros (5, 5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A and B are both rank deficient +%!test <48807> +%! A = A0.'; B = B0.'; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! [U, V, X, C, S, R]=gsvd (A, B); +%! D1 = zeros(3, 4); D1(1:3, 1:3) = C; +%! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6); + +## A: 5x3 complex full rank, B: 3x3 complex full rank +%!test <48807> +%! A0 = A0 + j*randn (5, 3); +%! B0 = diag ([1 2 4]) + j*diag ([4 -2 -1]); +%! A = A0; +%! B = B0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 3); D1(1:3, 1:3) = C; +%! D2 = S; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 5x3 complex full rank, B: 3x3 complex rank deficient +%!test <48807> +%! B(2, 2) = 0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; +%! D2 = [zeros(2, 1) S; zeros(1, 3)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 5x3 complex rank deficient, B: 3x3 complex full rank +%!test <48807> +%! B = B0; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 3); D1(1:3, 1:3) = C; +%! D2 = S; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A (5x3) and B (3x3) are both complex rank deficient +%!test <48807> +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(5, 2); D1(1:2, 1:2) = C; +%! D2 = [S; zeros(1, 2)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6); + +## A (now 3x5) complex and B (now 5x5) complex are full rank +## now, A is 3x5 +%!test <48807> +%! A = A0.'; +%! B0 = diag ([1 2 4 8 16]) + j*diag ([-5 4 -3 2 -1]); +%! B = B0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = [C zeros(3,2)]; +%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 3x5 complex full rank, B: 5x5 complex rank deficient +%!test <48807> +%! B(2, 2) = 0; +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; +%! D2 = zeros(5,5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A: 3x5 complex rank deficient, B: 5x5 complex full rank +%!test <48807> +%! B = B0; +%! A(3, :) = 2*A(1, :) - A(2, :); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(3, 5); D1(1:3, 1:3) = C; +%! D2 = zeros(5,5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); + +## A and B are both complex rank deficient +%!test <48807> +%! A = A0.'; +%! B = B0.'; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! [U, V, X, C, S, R] = gsvd (A, B); +%! D1 = zeros(3, 4); D1(1:3, 1:3) = C; +%! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; +%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); +%! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6); +%! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6); + +*/ +
--- a/libinterp/corefcn/module.mk Mon Oct 10 17:34:38 2016 -0400 +++ b/libinterp/corefcn/module.mk Tue Oct 11 08:56:03 2016 -0700 @@ -173,6 +173,7 @@ libinterp/corefcn/gl2ps-print.cc \ libinterp/corefcn/graphics.cc \ libinterp/corefcn/gripes.cc \ + libinterp/corefcn/gsvd.cc \ libinterp/corefcn/hash.cc \ libinterp/corefcn/help.cc \ libinterp/corefcn/hess.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gsvd.cc Tue Oct 11 08:56:03 2016 -0700 @@ -0,0 +1,371 @@ +// Copyright (C) 2016 Barbara Lócsi +// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 1996, 1997 John W. Eaton +// +// This program 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. +// +// This program 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 +// this program; if not, see <http://www.gnu.org/licenses/>. + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include <vector> + +#include "gsvd.h" + +#include "lo-error.h" +#include "lo-lapack-proto.h" +#include "dMatrix.h" +#include "fMatrix.h" +#include "CMatrix.h" +#include "fCMatrix.h" +#include "dDiagMatrix.h" +#include "fDiagMatrix.h" + +namespace octave +{ + namespace math + { + template <> + void + gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, + octave_idx_type m, octave_idx_type n, + octave_idx_type p, + octave_idx_type& k, octave_idx_type& l, + double *tmp_dataA, octave_idx_type m1, + double *tmp_dataB, octave_idx_type p1, Matrix& alpha, + Matrix& beta, double *u, octave_idx_type nrow_u, + double *v, octave_idx_type nrow_v, double *q, + octave_idx_type nrow_q, + Matrix& work, octave_idx_type* iwork, + octave_idx_type& info) + { + F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, tmp_dataA, m1, + tmp_dataB, p1, alpha.fortran_vec (), + beta.fortran_vec (), u, nrow_u, + v, nrow_v, q, nrow_q, work.fortran_vec (), + iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + + template <> + void + gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, + octave_idx_type m, octave_idx_type n, + octave_idx_type p, + octave_idx_type& k, octave_idx_type& l, + float *tmp_dataA, octave_idx_type m1, + float *tmp_dataB, octave_idx_type p1, + FloatMatrix& alpha, FloatMatrix& beta, float *u, + octave_idx_type nrow_u, float *v, + octave_idx_type nrow_v, float *q, + octave_idx_type nrow_q, FloatMatrix& work, + octave_idx_type* iwork, octave_idx_type& info) + { + F77_XFCN (sggsvd, SGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, tmp_dataA, m1, + tmp_dataB, p1, alpha.fortran_vec (), + beta.fortran_vec (), u, nrow_u, + v, nrow_v, q, nrow_q, work.fortran_vec (), + iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + + template <> + void + gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, + octave_idx_type m, octave_idx_type n, + octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, Complex *tmp_dataA, + octave_idx_type m1, Complex *tmp_dataB, + octave_idx_type p1, Matrix& alpha, + Matrix& beta, Complex *u, + octave_idx_type nrow_u, + Complex *v, octave_idx_type nrow_v, Complex *q, + octave_idx_type nrow_q, ComplexMatrix& work, + octave_idx_type* iwork, octave_idx_type& info) + { + Matrix rwork(2*n, 1); + F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, F77_DBLE_CMPLX_ARG (tmp_dataA), + m1, F77_DBLE_CMPLX_ARG (tmp_dataB), p1, + alpha.fortran_vec (), beta.fortran_vec (), + F77_DBLE_CMPLX_ARG (u), nrow_u, + F77_DBLE_CMPLX_ARG (v), nrow_v, + F77_DBLE_CMPLX_ARG (q), nrow_q, + F77_DBLE_CMPLX_ARG (work.fortran_vec ()), + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + + template <> + void + gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, + octave_idx_type m, octave_idx_type n, + octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, + FloatComplex *tmp_dataA, + octave_idx_type m1, + FloatComplex *tmp_dataB, + octave_idx_type p1, FloatMatrix& alpha, + FloatMatrix& beta, FloatComplex *u, + octave_idx_type nrow_u, FloatComplex *v, + octave_idx_type nrow_v, FloatComplex *q, + octave_idx_type nrow_q, + FloatComplexMatrix& work, + octave_idx_type* iwork, + octave_idx_type& info) + { + FloatMatrix rwork(2*n, 1); + F77_XFCN (cggsvd, CGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, F77_CMPLX_ARG (tmp_dataA), m1, + F77_CMPLX_ARG (tmp_dataB), p1, + alpha.fortran_vec (), beta.fortran_vec (), + F77_CMPLX_ARG (u), nrow_u, + F77_CMPLX_ARG (v), nrow_v, + F77_CMPLX_ARG (q), nrow_q, + F77_CMPLX_ARG (work.fortran_vec ()), + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + + template <typename T> + T + gsvd<T>::left_singular_matrix_A (void) const + { + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: U not computed because type == gsvd::sigma_only"); + return T (); + } + else + return left_smA; + } + + template <typename T> + T + gsvd<T>::left_singular_matrix_B (void) const + { + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: V not computed because type == gsvd::sigma_only"); + return T (); + } + else + return left_smB; + } + + template <typename T> + T + gsvd<T>::right_singular_matrix (void) const + { + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: X not computed because type == gsvd::sigma_only"); + return T (); + } + else + return right_sm; + } + + template <typename T> + T + gsvd<T>::R_matrix (void) const + { + if (type != gsvd::Type::std) + { + (*current_liboctave_error_handler) + ("gsvd: R not computed because type != gsvd::std"); + return T (); + } + else + return R; + } + + template <typename T> + gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type) + { + octave_idx_type info; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + octave_idx_type p = b.rows (); + + T atmp = a; + P *tmp_dataA = atmp.fortran_vec (); + + T btmp = b; + P *tmp_dataB = btmp.fortran_vec (); + + char jobu = 'U'; + char jobv = 'V'; + char jobq = 'Q'; + + octave_idx_type nrow_u = m; + octave_idx_type nrow_v = p; + octave_idx_type nrow_q = n; + + octave_idx_type k, l; + + switch (gsvd_type) + { + case gsvd<T>::Type::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = 'N'; + jobv = 'N'; + jobq = 'N'; + nrow_u = nrow_v = nrow_q = 1; + break; + + default: + break; + } + + type = gsvd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_smA.resize (nrow_u, m); + + P *u = left_smA.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + left_smB.resize (nrow_v, p); + + P *v = left_smB.fortran_vec (); + + if (! (jobq == 'N' || jobq == 'O')) + right_sm.resize (nrow_q, n); + + P *q = right_sm.fortran_vec (); + + octave_idx_type lwork = 3*n; + lwork = lwork > m ? lwork : m; + lwork = (lwork > p ? lwork : p) + n; + + T work (lwork, 1); + real_matrix alpha (n, 1); + real_matrix beta (n, 1); + + std::vector<octave_idx_type> iwork (n); + + gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l, + tmp_dataA, m, tmp_dataB, p, alpha, beta, u, + nrow_u, v, nrow_v, q, nrow_q, work, iwork.data (), info); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd"); + + if (info < 0) + (*current_liboctave_error_handler) ("*ggsvd.f: argument %d illegal", + -info); + else + { + if (info > 0) + (*current_liboctave_error_handler) + ("*ggsvd.f: Jacobi-type procedure failed to converge."); + else + { + octave_idx_type i, j; + + if (gsvd<T>::Type::std == gsvd_type) + { + R.resize(k+l, k+l); + int astart = n-k-l; + if (m - k - l >= 0) + { + astart = n-k-l; + // R is stored in A(1:K+L,N-K-L+1:N) + for (i = 0; i < k+l; i++) + for (j = 0; j < k+l; j++) + R.xelem (i, j) = atmp.xelem (i, astart + j); + } + else + { + // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), + // ( 0 R22 R23 ) + + for (i = 0; i < m; i++) + for (j = 0; j < k+l; j++) + R.xelem (i, j) = atmp.xelem (i, astart + j); + // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N) + for (i = k+l-1; i >=m; i--) + { + for (j = 0; j < m; j++) + R.xelem(i, j) = 0.0; + for (j = m; j < k+l; j++) + R.xelem (i, j) = btmp.xelem (i - k, astart + j); + } + } + } + + if (m-k-l >= 0) + { + // Fills in C and S + sigmaA.resize (l, l); + sigmaB.resize (l, l); + for (i = 0; i < l; i++) + { + sigmaA.dgxelem(i) = alpha.elem(k+i); + sigmaB.dgxelem(i) = beta.elem(k+i); + } + } + else + { + // Fills in C and S + sigmaA.resize (m-k, m-k); + sigmaB.resize (m-k, m-k); + for (i = 0; i < m-k; i++) + { + sigmaA.dgxelem(i) = alpha.elem(k+i); + sigmaB.dgxelem(i) = beta.elem(k+i); + } + } + } + } + } + + // Instantiations we need. + template class gsvd<Matrix>; + template class gsvd<FloatMatrix>; + template class gsvd<ComplexMatrix>; + template class gsvd<FloatComplexMatrix>; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gsvd.h Tue Oct 11 08:56:03 2016 -0700 @@ -0,0 +1,102 @@ +// Copyright (C) 2016 Barbara Lócsi +// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 1996, 1997 John W. Eaton +// +// This program 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. +// +// This program 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 +// this program; if not, see <http://www.gnu.org/licenses/>. + +#if !defined (octave_gsvd_h) +#define octave_gsvd_h 1 + +#include "octave-config.h" + +namespace octave +{ + namespace math + { + template <typename T> + class + gsvd + { + public: + + enum class Type + { + std, + economy, + sigma_only + }; + + gsvd (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () + { } + + gsvd (const T& a, const T& b, + gsvd::Type gsvd_type = gsvd<T>::Type::economy); + + gsvd (const gsvd& a) + : type (a.type), + sigmaA (a.sigmaA), sigmaB (a.sigmaB), + left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm), + R(a.R) { } + + gsvd& operator = (const gsvd& a) + { + if (this != &a) + { + type = a.type; + sigmaA = a.sigmaA; + sigmaB = a.sigmaB; + left_smA = a.left_smA; + left_smB = a.left_smB; + right_sm = a.right_sm; + R = a.R; + } + + return *this; + } + + ~gsvd (void) { } + + typename T::real_diag_matrix_type + singular_values_A (void) const { return sigmaA; } + + typename T::real_diag_matrix_type + singular_values_B (void) const { return sigmaB; } + + T left_singular_matrix_A (void) const; + T left_singular_matrix_B (void) const; + + T right_singular_matrix (void) const; + T R_matrix (void) const; + + private: + typedef typename T::value_type P; + typedef typename T::real_matrix_type real_matrix; + + gsvd::Type type; + typename T::real_diag_matrix_type sigmaA, sigmaB; + T left_smA, left_smB; + T right_sm, R; + + void ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m, + octave_idx_type n, octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, P *tmp_dataA, octave_idx_type m1, + P *tmp_dataB, octave_idx_type p1, real_matrix& alpha, + real_matrix& beta, P *u, octave_idx_type nrow_u, P *v, + octave_idx_type nrow_v, P *q, octave_idx_type nrow_q, T& work, + octave_idx_type* iwork, octave_idx_type& info); + }; + } +} + +#endif
--- a/liboctave/numeric/lo-lapack-proto.h Mon Oct 10 17:34:38 2016 -0400 +++ b/liboctave/numeric/lo-lapack-proto.h Tue Oct 11 08:56:03 2016 -0700 @@ -807,6 +807,126 @@ F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + // GGSVD + + F77_RET_T + F77_FUNC (dggsvd, DGGSVD) + (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_CONST_CHAR_ARG_DECL, // JOBV + F77_CONST_CHAR_ARG_DECL, // JOBQ + const F77_INT&, // M + const F77_INT&, // N + const F77_INT&, // P + F77_INT &, // K + F77_INT &, // L + F77_DBLE*, // A(LDA,N) + const F77_INT&, // LDA + F77_DBLE*, // B(LDB,N) + const F77_INT&, // LDB + F77_DBLE*, // ALPHA(N) + F77_DBLE*, // BETA(N) + F77_DBLE*, // U(LDU,M) + const F77_INT&, // LDU + F77_DBLE*, // V(LDV,P) + const F77_INT&, // LDV + F77_DBLE*, // Q(LDQ,N) + const F77_INT&, // LDQ + F77_DBLE*, // WORK + F77_INT*, // IWORK(N) + F77_INT& // INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (sggsvd, SGGSVD) + (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_CONST_CHAR_ARG_DECL, // JOBV + F77_CONST_CHAR_ARG_DECL, // JOBQ + const F77_INT&, // M + const F77_INT&, // N + const F77_INT&, // P + F77_INT &, // K + F77_INT &, // L + F77_REAL*, // A + const F77_INT&, // LDA + F77_REAL*, // B + const F77_INT&, // LDB + F77_REAL*, // ALPHA + F77_REAL*, // BETA + F77_REAL*, // U + const F77_INT&, // LDU + F77_REAL*, // V + const F77_INT&, // LDV + F77_REAL*, // Q + const F77_INT&, // LDQ + F77_REAL*, // WORK + F77_INT*, // IWORK + F77_INT& // INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (zggsvd, ZGGSVD) + (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_CONST_CHAR_ARG_DECL, // JOBV + F77_CONST_CHAR_ARG_DECL, // JOBQ + const F77_INT&, // M + const F77_INT&, // N + const F77_INT&, // P + F77_INT &, // K + F77_INT &, // L + F77_DBLE_CMPLX*, // A(LDA,N) + const F77_INT&, // LDA + F77_DBLE_CMPLX*, // B(LDB,N) + const F77_INT&, // LDB + F77_DBLE*, // ALPHA(N) + F77_DBLE*, // BETA(N) + F77_DBLE_CMPLX*, // U(LDU,M) + const F77_INT&, // LDU + F77_DBLE_CMPLX*, // V(LDV,P) + const F77_INT&, // LDV + F77_DBLE_CMPLX*, // Q(LDQ,N) + const F77_INT&, // LDQ + F77_DBLE_CMPLX*, // WORK + F77_DBLE*, // RWORK + F77_INT*, // IWORK(N) + F77_INT& // INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (cggsvd, CGGSVD) + (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_CONST_CHAR_ARG_DECL, // JOBV + F77_CONST_CHAR_ARG_DECL, // JOBQ + const F77_INT&, // M + const F77_INT&, // N + const F77_INT&, // P + F77_INT &, // K + F77_INT &, // L + F77_CMPLX*, // A + const F77_INT&, // LDA + F77_CMPLX*, // B + const F77_INT&, // LDB + F77_REAL*, // ALPHA + F77_REAL*, // BETA + F77_CMPLX*, // U + const F77_INT&, // LDU + F77_CMPLX*, // V + const F77_INT&, // LDV + F77_CMPLX*, // Q + const F77_INT&, // LDQ + F77_CMPLX*, // WORK + F77_REAL*, // RWORK + F77_INT*, // IWORK + F77_INT& // INFO + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); + // GTSV F77_RET_T
--- a/liboctave/numeric/module.mk Mon Oct 10 17:34:38 2016 -0400 +++ b/liboctave/numeric/module.mk Tue Oct 11 08:56:03 2016 -0700 @@ -18,6 +18,7 @@ liboctave/numeric/DASSL.h \ liboctave/numeric/DET.h \ liboctave/numeric/EIG.h \ + liboctave/numeric/gsvd.h \ liboctave/numeric/LSODE.h \ liboctave/numeric/ODE.h \ liboctave/numeric/ODEFunc.h \ @@ -69,6 +70,7 @@ liboctave/numeric/DASRT.cc \ liboctave/numeric/DASSL.cc \ liboctave/numeric/EIG.cc \ + liboctave/numeric/gsvd.cc \ liboctave/numeric/LSODE.cc \ liboctave/numeric/ODES.cc \ liboctave/numeric/Quad.cc \
--- a/liboctave/operators/mx-defs.h Mon Oct 10 17:34:38 2016 -0400 +++ b/liboctave/operators/mx-defs.h Tue Oct 11 08:56:03 2016 -0700 @@ -66,6 +66,8 @@ class EIG; +template <typename T> class gsvd; + template <typename T> class hess; template <typename T> class schur;
--- a/liboctave/operators/mx-ext.h Mon Oct 10 17:34:38 2016 -0400 +++ b/liboctave/operators/mx-ext.h Tue Oct 11 08:56:03 2016 -0700 @@ -57,6 +57,11 @@ #include "EIG.h" +// Result of a Generalized Singular Value Decomposition. + +#include "gsvd.h" + + // Result of an LU decomposition. #include "lu.h"
--- a/scripts/help/__unimplemented__.m Mon Oct 10 17:34:38 2016 -0400 +++ b/scripts/help/__unimplemented__.m Tue Oct 11 08:56:03 2016 -0700 @@ -49,11 +49,6 @@ txt = ["exifread is deprecated. " ... "The functionality is available in the imfinfo function."]; - case "gsvd" - txt = ["gsvd is not currently part of core Octave. ", ... - "See the linear-algebra package at ", ... - "@url{http://octave.sourceforge.net/linear-algebra/}."]; - case "funm" txt = ["funm is not currently part of core Octave. ", ... "See the linear-algebra package at ", ... @@ -657,7 +652,6 @@ "graph", "graymon", "griddedInterpolant", - "gsvd", "guide", "h5create", "h5disp",
--- a/scripts/plot/util/private/__gnuplot_draw_axes__.m Mon Oct 10 17:34:38 2016 -0400 +++ b/scripts/plot/util/private/__gnuplot_draw_axes__.m Tue Oct 11 08:56:03 2016 -0700 @@ -500,7 +500,10 @@ while (! isempty (kids)) - obj = get (kids(end)); + h_obj = kids(end); + kids = kids(1:(end-1)); + + obj = get (h_obj); if (isfield (obj, "xdata")) obj.xdata = double (obj.xdata); @@ -515,13 +518,12 @@ if (isfield (obj, "units")) units = obj.units; unwind_protect - set (kids(end), "units", "data"); - obj = get (kids(end)); + set (h_obj, "units", "data"); + obj = get (h_obj); unwind_protect_cleanup - set (kids(end), "units", units); + set (h_obj, "units", units); end_unwind_protect endif - kids = kids(1:(end-1)); if (strcmp (obj.visible, "off")) continue; @@ -540,10 +542,6 @@ obj.zdata(obj.zdata<=0) = NaN; endif - if (strcmp (get (obj.parent, "type"), "hggroup")) - obj.displayname = get (obj.parent, "displayname"); - endif - switch (obj.type) case "image" img_data = mapcdata (obj.cdata, obj.cdatamapping, clim, cmap_sz); @@ -634,7 +632,6 @@ titlespec{data_idx} = ['title "' tmp '"']; endif usingclause{data_idx} = sprintf ("record=%d", numel (obj.xdata)); - errbars = ""; if (nd == 3) xdat = obj.xdata(:); ydat = obj.ydata(:); @@ -657,7 +654,7 @@ endif [style, sidx] = do_linestyle_command (obj, obj.color, data_idx, - plot_stream, errbars); + plot_stream); if isempty (style{1}) style{1} = "points"; @@ -691,6 +688,43 @@ style{3}, sidx(3)); endif + if (strcmp (get (obj.parent, "type"), "hggroup")) + hg = get (obj.parent, "children"); + if (hg(1) == h_obj) + # Place phantom errorbar data for legend + data_idx += 1; + is_image_data(data_idx) = is_image_data(data_idx - 1); + parametric(data_idx) = parametric(data_idx - 1); + have_cdata(data_idx) = have_cdata(data_idx - 1); + have_3d_patch(data_idx) = have_3d_patch(data_idx - 1); + obj.displayname = get (obj.parent, "displayname"); + if (isempty (get (obj.parent, "displayname"))) + titlespec{data_idx} = "title \"\""; + else + tmp = undo_string_escapes ( + __maybe_munge_text__ (enhanced, obj, "displayname") + ); + titlespec{data_idx} = ['title "' tmp '"']; + endif + data{data_idx} = nan (4,1); + usingclause{data_idx} = sprintf ("record=1 using ($1):($2):($3):($4)"); + switch (get (obj.parent, "format")) + case {"box" "boxy" "boxxy"} + errbars = "boxxy"; + case "xyerr" + errbars = "xyerrorbars"; + case "yerr" + errbars = "yerrorbars"; + case "xerr" + errbars = "xerrorbars"; + otherwise + errbars = "xerrorbars"; + endswitch + withclause{data_idx} = sprintf ("with %s linestyle %d", + errbars, sidx(1)); + endif + endif + case "patch" [nr, nc] = size (obj.xdata); @@ -1850,7 +1884,7 @@ endfunction function [style, ltidx] = do_linestyle_command (obj, linecolor, idx, - plot_stream, errbars = "") + plot_stream) idx = idx + 8; style = {}; ltidx = []; @@ -1882,9 +1916,6 @@ if (! isempty (lt)) fprintf (plot_stream, " %s", lt); endif - if (! isempty (errbars)) - found_style = true; - endif if (isfield (obj, "linewidth")) fprintf (plot_stream, " linewidth %f", obj.linewidth); @@ -1898,111 +1929,105 @@ endif sidx = 1; - if (isempty (errbars)) - if (isempty (lt)) - style{sidx} = ""; - else - style{sidx} = "lines"; - endif - ltidx(sidx) = idx; + if (isempty (lt)) + style{sidx} = ""; + else + style{sidx} = "lines"; + endif + ltidx(sidx) = idx; - facesame = true; - if (! isequal (pt, pt2) && isfield (obj, "markerfacecolor") - && ! strcmp (obj.markerfacecolor, "none")) - if (strcmp (obj.markerfacecolor, "auto") - || (isnumeric (obj.markerfacecolor) - && isequal (color, obj.markerfacecolor))) - if (! isempty (pt2)) - fprintf (plot_stream, " pointtype %s", pt2); - style{sidx} = [style{sidx} "points"]; - endif - if (isfield (obj, "markersize")) - fprintf (plot_stream, " pointsize %f", obj.markersize / 3); - endif + facesame = true; + if (! isequal (pt, pt2) && isfield (obj, "markerfacecolor") + && ! strcmp (obj.markerfacecolor, "none")) + if (strcmp (obj.markerfacecolor, "auto") + || (isnumeric (obj.markerfacecolor) + && isequal (color, obj.markerfacecolor))) + if (! isempty (pt2)) + fprintf (plot_stream, " pointtype %s", pt2); + style{sidx} = [style{sidx} "points"]; + endif + if (isfield (obj, "markersize")) + fprintf (plot_stream, " pointsize %f", obj.markersize / 3); + endif + else + facesame = false; + if (! found_style) + fputs (plot_stream, " default"); + endif + fputs (plot_stream, ";\n"); + if (! isempty (style{sidx})) + sidx += 1; + idx += 1; else - facesame = false; - if (! found_style) - fputs (plot_stream, " default"); - endif fputs (plot_stream, ";\n"); - if (! isempty (style{sidx})) - sidx += 1; - idx += 1; - else - fputs (plot_stream, ";\n"); - endif - fprintf (plot_stream, "set style line %d default;\n", idx); - fprintf (plot_stream, "set style line %d", idx); - if (isnumeric (obj.markerfacecolor)) - fprintf (plot_stream, " linecolor rgb \"#%02x%02x%02x\"", - round (255*obj.markerfacecolor)); - else - fprintf (plot_stream, " palette"); - endif - if (! isempty (pt2)) - style{sidx} = "points"; - ltidx(sidx) = idx; - fprintf (plot_stream, " pointtype %s", pt2); - endif - if (isfield (obj, "markersize")) - fprintf (plot_stream, " pointsize %f", obj.markersize / 3); - endif + endif + fprintf (plot_stream, "set style line %d default;\n", idx); + fprintf (plot_stream, "set style line %d", idx); + if (isnumeric (obj.markerfacecolor)) + fprintf (plot_stream, " linecolor rgb \"#%02x%02x%02x\"", + round (255*obj.markerfacecolor)); + else + fprintf (plot_stream, " palette"); + endif + if (! isempty (pt2)) + style{sidx} = "points"; + ltidx(sidx) = idx; + fprintf (plot_stream, " pointtype %s", pt2); + endif + if (isfield (obj, "markersize")) + fprintf (plot_stream, " pointsize %f", obj.markersize / 3); endif endif - if (! isempty(pt) && isfield (obj, "markeredgecolor") - && ! strcmp (obj.markeredgecolor, "none")) - if (facesame && ! isempty (pt) - && (strcmp (obj.markeredgecolor, "auto") - || (isnumeric (obj.markeredgecolor) - && isequal (color, obj.markeredgecolor)))) - if (sidx == 1 && ((length (style{sidx}) == 5 - && strncmp (style{sidx}, "lines", 5)) || isempty (style{sidx}))) - if (! isempty (pt)) - style{sidx} = [style{sidx} "points"]; - fprintf (plot_stream, " pointtype %s", pt); - endif - if (isfield (obj, "markersize")) - fprintf (plot_stream, " pointsize %f", obj.markersize / 3); - endif - endif - else - if (! found_style) - fputs (plot_stream, " default"); - endif - fputs (plot_stream, ";\n"); - if (! isempty (style{sidx})) - sidx += 1; - idx += 1; - else - fputs (plot_stream, ";\n"); - endif - fprintf (plot_stream, "set style line %d default;\n", idx); - fprintf (plot_stream, "set style line %d", idx); - if (isnumeric (obj.markeredgecolor) || strcmp (obj.markeredgecolor, "auto")) - if (isnumeric (obj.markeredgecolor)) - edgecolor = obj.markeredgecolor; - else - edgecolor = obj.color; - end - fprintf (plot_stream, " linecolor rgb \"#%02x%02x%02x\"", - round (255*edgecolor)); - else - fprintf (plot_stream, " palette"); - endif + endif + if (! isempty(pt) && isfield (obj, "markeredgecolor") + && ! strcmp (obj.markeredgecolor, "none")) + if (facesame && ! isempty (pt) + && (strcmp (obj.markeredgecolor, "auto") + || (isnumeric (obj.markeredgecolor) + && isequal (color, obj.markeredgecolor)))) + if (sidx == 1 && ((length (style{sidx}) == 5 + && strncmp (style{sidx}, "lines", 5)) || isempty (style{sidx}))) if (! isempty (pt)) - style{sidx} = "points"; - ltidx(sidx) = idx; + style{sidx} = [style{sidx} "points"]; fprintf (plot_stream, " pointtype %s", pt); endif if (isfield (obj, "markersize")) fprintf (plot_stream, " pointsize %f", obj.markersize / 3); endif endif + else + if (! found_style) + fputs (plot_stream, " default"); + endif + fputs (plot_stream, ";\n"); + if (! isempty (style{sidx})) + sidx += 1; + idx += 1; + else + fputs (plot_stream, ";\n"); + endif + fprintf (plot_stream, "set style line %d default;\n", idx); + fprintf (plot_stream, "set style line %d", idx); + if (isnumeric (obj.markeredgecolor) || strcmp (obj.markeredgecolor, "auto")) + if (isnumeric (obj.markeredgecolor)) + edgecolor = obj.markeredgecolor; + else + edgecolor = obj.color; + end + fprintf (plot_stream, " linecolor rgb \"#%02x%02x%02x\"", + round (255*edgecolor)); + else + fprintf (plot_stream, " palette"); + endif + if (! isempty (pt)) + style{sidx} = "points"; + ltidx(sidx) = idx; + fprintf (plot_stream, " pointtype %s", pt); + endif + if (isfield (obj, "markersize")) + fprintf (plot_stream, " pointsize %f", obj.markersize / 3); + endif endif - else - style{1} = errbars; - ltidx(1) = idx; - fputs (plot_stream, " pointtype 0"); endif if (! found_style && isempty (style{1})) @@ -2106,7 +2131,7 @@ ## FIXME: Should be 6 pt start, using "*" instead pt = pt2 = "3"; case "none" - pt = pt2 = ""; + pt = pt2 = "-1"; otherwise pt = pt2 = ""; endswitch
--- a/scripts/polynomial/residue.m Mon Oct 10 17:34:38 2016 -0400 +++ b/scripts/polynomial/residue.m Tue Oct 11 08:56:03 2016 -0700 @@ -333,21 +333,10 @@ endfor ## Add the direct term. - if (numel (k)) pnum += conv (pden, k); endif - ## Check for leading zeros and trim the polynomial coefficients. - if (isa (r, "single") || isa (p, "single") || isa (k, "single")) - small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single"); - else - small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; - endif - - pnum(abs (pnum) < small) = 0; - pden(abs (pden) < small) = 0; - pnum = polyreduce (pnum); pden = polyreduce (pden); @@ -383,7 +372,7 @@ %! assert (isempty (k)); %! assert (e, [1; 2; 1; 2]); %! [br, ar] = residue (r, p, k); -%! assert (br, b, 1e-12); +%! assert (br, [0,b], 1e-12); %! assert (ar, a, 1e-12); %!test @@ -431,6 +420,14 @@ %! b = [1, z1]; %! a = [1, -(p1 + p2), p1*p2, 0, 0]; %! [br, ar] = residue (r, p, k, e); -%! assert (br, b, 1e-8); +%! assert (br, [0,0,b], 1e-7); %! assert (ar, a, 1e-8); +%!test <49291> +%! rf = [1e3, 2e3, 1e3, 2e3]; +%! cf = [316.2e-9, 50e-9, 31.6e-9, 5e-9]; +%! [num, den] = residue (1./cf,-1./(rf.*cf),0); +%! assert (numel (num), 4); +%! assert (numel (den), 5); +%! assert (den(1), 1); +