Mercurial > octave
changeset 22595:acfb81e6992a
maint: merge stable to default
| author | Carlo de Falco <carlo.defalco@polimi.it> |
|---|---|
| date | Thu, 06 Oct 2016 07:36:59 +0200 |
| parents | 354b7a6e642c (diff) b8d525710075 (current diff) |
| children | f812283c4367 |
| files | scripts/ode/private/integrate_const.m scripts/ode/private/integrate_n_steps.m scripts/ode/private/ode_struct_value_check.m |
| diffstat | 13 files changed, 1061 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Thu Oct 06 07:13:24 2016 +0200 +++ b/NEWS Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/configure.ac Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/doc/interpreter/linalg.txi Thu Oct 06 07:36:59 2016 +0200 @@ -108,6 +108,8 @@ @DOCSTRING(givens) +@DOCSTRING(gsvd) + @DOCSTRING(planerot) @DOCSTRING(inv)
--- a/libinterp/corefcn/data.cc Thu Oct 06 07:13:24 2016 +0200 +++ b/libinterp/corefcn/data.cc Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/libinterp/corefcn/module.mk Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/liboctave/numeric/lo-lapack-proto.h Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/liboctave/numeric/module.mk Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/liboctave/operators/mx-defs.h Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/liboctave/operators/mx-ext.h Thu Oct 06 07:36:59 2016 +0200 @@ -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 Thu Oct 06 07:13:24 2016 +0200 +++ b/scripts/help/__unimplemented__.m Thu Oct 06 07:36:59 2016 +0200 @@ -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",