Mercurial > octave
changeset 22576:cdbf931c2024
maint: merge stable to default.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Sat, 01 Oct 2016 20:35:36 +0100 |
| parents | 17f6ed063c69 (current diff) 726f32d9fd92 (diff) |
| children | 98eeed41f372 |
| files | libinterp/corefcn/__gsvd__.cc scripts/linear-algebra/gsvd.m |
| diffstat | 8 files changed, 481 insertions(+), 458 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__gsvd__.cc Fri Sep 30 17:45:22 2016 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,350 +0,0 @@ -// 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}, @var{r}] =} __gsvd__ (@var{a}, @var{b}) -Undocumented internal function. -@end deftypefn */) -{ - if (args.length () != 2) - print_usage (); - - 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 (); - - // 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 (); - - if (argA.is_real_type () && argB.is_real_type ()) - { - Matrix tmpA = argA.matrix_value (); - Matrix tmpB = argB.matrix_value (); - - // FIXME: This code is still using error_state - if (! error_state) - { - if (tmpA.any_element_is_inf_or_nan ()) - error ("gsvd: B 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.complex_matrix_value (); - ComplexMatrix ctmpB = argB.complex_matrix_value (); - - if (! error_state) - { - 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; -} - -/* -## FIXME: All tests are commented out for the 4.2.0 release. -## The m-file gsvd.m needs to be replaced with C++ code that achieves Matlab -## compatible outputs, and the BIST tests need to be updated to reflect the new -## outputs. - -## 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 -%! [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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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 -%! 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); -*/ -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/gsvd.cc Sat Oct 01 20:35:36 2016 +0100 @@ -0,0 +1,465 @@ +// 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}, @var{r}] =} gsvd (@var{a}, @var{b}) +@cindex generalized singular value decomposition +Compute the generalized singular value decomposition of (@var{a}, @var{b}): +@tex +$$ + U^H A X = [I 0; 0 C] [0 R] + V^H B X = [0 S; 0 0] [0 R] + C*C + S*S = eye (columns (A)) + I and 0 are padding matrices of suitable size + R is upper triangular +$$ +@end tex +@ifinfo + +@example +@group +u' * a * x = [I 0; 0 c] * [0 r] +v' * b * x = [0 s; 0 0] * [0 r] +c * c + s * s = eye (columns (a)) +I and 0 are padding matrices of suitable size +r is upper triangular +@end group +@end example + +@end ifinfo + +The function @code{gsvd} normally returns the vector of generalized singular +values +@tex +diag (C) ./ diag (S). +@end tex +@ifinfo +diag (r) ./ diag (s). +@end ifinfo +If asked for five return values, it computes +@tex +$U$, $V$, and $X$. +@end tex +@ifinfo +U, V, and X. +@end ifinfo +With a sixth output argument, it also returns +@tex +R, +@end tex +@ifinfo +r, +@end ifinfo +The common upper triangular right term. Other authors, like +@nospell{S. Van Huffel}, define this transformation as the simultaneous +diagonalization of the input matrices, this can be achieved by multiplying +@tex +X +@end tex +@ifinfo +x +@end ifinfo +by the inverse of +@tex +[I 0; 0 R]. +@end tex +@ifinfo +[I 0; 0 r]. +@end ifinfo + +For example, + +@example +gsvd (hilb (3), [1 2 3; 3 2 1]) + +@result{} + 0.1055705 + 0.0031759 +@end example + +@noindent +and + +@example +[u, v, c, s, x, r] = gsvd (hilb (3), [1 2 3; 3 2 1]) +@result{} + +u = + + -0.965609 0.240893 0.097825 + -0.241402 -0.690927 -0.681429 + -0.096561 -0.681609 0.725317 + +v = + + -0.41974 0.90765 + -0.90765 -0.41974 + +x = + + 0.408248 0.902199 0.139179 + -0.816497 0.429063 -0.386314 + 0.408248 -0.044073 -0.911806 + +c = + + 0.10499 0.00000 + 0.00000 0.00318 + +s = + 0.99447 0.00000 + 0.00000 0.99999 + +r = + -0.14093 -1.24345 0.43737 + 0.00000 -3.90043 2.57818 + 0.00000 0.00000 -2.52599 + +@end example + +The code is a wrapper to the corresponding @sc{lapack} dggsvd and zggsvd +routines. + +@end deftypefn */) +{ + if (args.length () != 2) + print_usage (); + + 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 (); + + // 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 (); + + if (argA.is_real_type () && argB.is_real_type ()) + { + Matrix tmpA = argA.matrix_value (); + Matrix tmpB = argB.matrix_value (); + + // FIXME: This code is still using error_state + if (! error_state) + { + if (tmpA.any_element_is_inf_or_nan ()) + error ("gsvd: B 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.complex_matrix_value (); + ComplexMatrix ctmpB = argB.complex_matrix_value (); + + if (! error_state) + { + 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 + { + // Actually, can't tell which arg is at fault + err_wrong_type_arg ("gsvd", argA); + //err_wrong_type_arg ("gsvd", argB); + } + } + + return retval; +} + +/* +## 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 +%! [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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 +%! 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 Fri Sep 30 17:45:22 2016 -0700 +++ b/libinterp/corefcn/module.mk Sat Oct 01 20:35:36 2016 +0100 @@ -118,7 +118,6 @@ libinterp/corefcn/__contourc__.cc \ libinterp/corefcn/__dispatch__.cc \ libinterp/corefcn/__dsearchn__.cc \ - libinterp/corefcn/__gsvd__.cc \ libinterp/corefcn/__ichol__.cc \ libinterp/corefcn/__ilu__.cc \ libinterp/corefcn/__lin_interpn__.cc \ @@ -174,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 \
--- a/scripts/linear-algebra/gsvd.m Fri Sep 30 17:45:22 2016 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ -## Copyright (C) 2016 Rik Wehbring -## -## 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{S} =} gsvd (@var{A}, @var{B}) -## @deftypefnx {} {[@var{U}, @var{V}, @var{X}, @var{C}, @var{S}] =} gsvd (@var{A}, @var{B}) -## 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 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$, and $X$. -## @end tex -## @ifnottex -## U, V, and X. -## @end ifnottex -## -## The code is a wrapper to the corresponding @sc{lapack} dggsvd and zggsvd -## routines. -## -## @seealso{svd} -## @end deftypefn - -## FIXME: This m-file is a wrapper around __gsvd__ in libinterp/corefcn. -## It was put in place strictly for the 4.2.0 release in order to achieve -## Matlab-compatible output for the gsvd function. Eventually the C++ code -## needs to be updated to reflect what is being done in this m-file and then -## this m-file should be deleted. - -function [U, V, X, C, S] = gsvd (A, B, econ = false) - - if (nargin < 2 || nargin > 3) - print_usage (); - endif - - if (nargin == 3) - warning ('gsvd: "economy" option is not yet implemented'); - endif - - if (nargout <= 1) - U = __gsvd__ (A, B); - else - [U, V, X, C, S, R] = __gsvd__ (A, B); - X = (R / X)'; - [m, n] = size (A); - if (m > n) - C = [C; zeros(m-n, n)]; - elseif (m < n) - C = [C, zeros(m, n-m)]; - S0 = S; - S = eye (n); - S(1:m, 1:m) = S0; - endif - endif - -endfunction - - -## FIXME: All BIST tests are in the C++ file __gsvd__.cc. They are currently -## not run because they have not been updated to reflect the actual expected -## output.
--- a/scripts/linear-algebra/module.mk Fri Sep 30 17:45:22 2016 -0700 +++ b/scripts/linear-algebra/module.mk Sat Oct 01 20:35:36 2016 +0100 @@ -9,7 +9,6 @@ scripts/linear-algebra/cross.m \ scripts/linear-algebra/duplication_matrix.m \ scripts/linear-algebra/expm.m \ - scripts/linear-algebra/gsvd.m \ scripts/linear-algebra/housh.m \ scripts/linear-algebra/isbanded.m \ scripts/linear-algebra/isdefinite.m \
--- a/scripts/plot/util/print.m Fri Sep 30 17:45:22 2016 -0700 +++ b/scripts/plot/util/print.m Sat Oct 01 20:35:36 2016 +0100 @@ -160,7 +160,9 @@ ## ## @item fig ## XFig. For the Gnuplot graphics toolkit, the additional options -## @option{-textspecial} or @option{-textnormal} can be used to control whether the special flag should be set for the text in the figure. (default is @option{-textnormal}) +## @option{-textspecial} or @option{-textnormal} can be used to control +## whether the special flag should be set for the text in the figure. +## (default is @option{-textnormal}) ## ## @item gif ## GIF image (only available for the Gnuplot graphics toolkit)
--- a/scripts/plot/util/private/__gnuplot_print__.m Fri Sep 30 17:45:22 2016 -0700 +++ b/scripts/plot/util/private/__gnuplot_print__.m Sat Oct 01 20:35:36 2016 +0100 @@ -115,7 +115,8 @@ local_drawnow ([opts.devopt " " gp_opts], opts.name, opts); case {"cairolatex", "epscairo", "epscairolatex", ... "epscairolatexstandalone", "pdfcairo", "pdfcairolatex", ... - "pdfcairolatexstandalone", "pngcairo"} + "pdfcairolatexstandalone", "pngcairo", ... + "pdflatexstandalone", "pdflatex"} term = opts.devopt; if (strfind (term, "standalone")) gp_opts = sprintf ("standalone %s", gp_opts); @@ -124,9 +125,10 @@ if (strfind (term, "epscairolatex")) gp_opts = sprintf ("eps %s", gp_opts); term = strrep (term, "epscairolatex", "cairolatex"); - elseif (strfind (term, "pdfcairolatex")) + elseif (strfind (term, "pdfcairolatex") || strfind (term, "pdflatex")) gp_opts = sprintf ("pdf %s", gp_opts); term = strrep (term, "pdfcairolatex", "cairolatex"); + term = strrep (term, "pdflatex", "cairolatex"); endif if (__gnuplot_has_terminal__ (term)) local_drawnow ([term " " gp_opts], opts.name, opts); @@ -273,7 +275,8 @@ endif case {"cairolatex", "epscairo", "epscairolatex", ... "epscairolatexstandalone", "pdfcairo", "pdfcairolatex", ... - "pdfcairolatexstandalone", "pngcairo"} + "pdfcairolatexstandalone", "pngcairo", ... + "pdflatexstandalone", "pdflatex"} if (! isempty (opts.font) && ! isempty (opts.fontsize)) f = sprintf ('font "%s,%d"', opts.font, opts.fontsize); elseif (! isempty (opts.font))
--- a/scripts/plot/util/private/__print_parse_opts__.m Fri Sep 30 17:45:22 2016 -0700 +++ b/scripts/plot/util/private/__print_parse_opts__.m Sat Oct 01 20:35:36 2016 +0100 @@ -228,7 +228,9 @@ "pstex", "tiff", "tiffn" "tikz", "pcxmono", ... "pcx24b", "pcx256", "pcx16", "pgm", "pgmraw", ... "ppm", "ppmraw", "pdflatex", "texdraw", ... - "pdfcairo", "pngcairo", "pstricks", ... + "epscairo", "pdfcairo", "pngcairo", "cairolatex", ... + "pdfcairolatex", "pdfcairolatexstandalone", ... + "epscairolatex", "epscairolatexstandalone", "pstricks", ... "epswrite", "eps2write", "pswrite", "ps2write", "pdfwrite"}; suffixes = {"ai", "cdr", "fig", "png", "jpg", ... @@ -240,7 +242,9 @@ "ps", "tiff", "tiff", "tikz", "pcx", ... "pcx", "pcx", "pcx", "pgm", "pgm", ... "ppm", "ppm", "tex", "tex", ... - "pdf", "png", "tex", ... + "eps", "pdf", "png", "tex", ... + "tex", "tex", ... + "tex", "tex", "tex", ... "eps", "eps", "ps", "ps", "pdf"}; if (isfigure (arg_st.figure))