Mercurial > octave
changeset 22573:d29f3679e954 stable
gsvd: remove function from the upcoming 4.2 release (bug #48807)
* libinterp/corefcn/gsvd.cc: remove Octave gsvd function because it is
unfinished work and still has Matlab incompatibility issues whose later
fix would become backwards compatibility issues.
* liboctave/numeric/gsvd.cc, liboctave/numeric/gsvd.h,
liboctave/numeric/lo-lapack-proto.h, liboctave/operators/mx-defs.h,
liboctave/operators/mx-ext.h: remove liboctave gsvd class for same
reason.
* libinterp/corefcn/module.mk, liboctave/numeric/module.mk: remove
files from build system.
* doc/interpreter/linalg.txi: remove reference to gsvd on manual.
* NEWS: remove mention of gsvd
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Sat, 01 Oct 2016 17:31:13 +0100 |
| parents | d65ed83dfd13 |
| children | 730b385e5646 |
| files | NEWS doc/interpreter/linalg.txi libinterp/corefcn/gsvd.cc libinterp/corefcn/module.mk liboctave/numeric/gsvd.cc liboctave/numeric/gsvd.h liboctave/numeric/lo-lapack-proto.h liboctave/numeric/module.mk liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 10 files changed, 0 insertions(+), 1071 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sat Oct 01 17:02:51 2016 +0100 +++ b/NEWS Sat Oct 01 17:31:13 2016 +0100 @@ -164,7 +164,6 @@ deg2rad dialog evalc - gsvd hash im2double isocaps
--- a/doc/interpreter/linalg.txi Sat Oct 01 17:02:51 2016 +0100 +++ b/doc/interpreter/linalg.txi Sat Oct 01 17:31:13 2016 +0100 @@ -108,8 +108,6 @@ @DOCSTRING(givens) -@DOCSTRING(gsvd) - @DOCSTRING(planerot) @DOCSTRING(inv)
--- a/libinterp/corefcn/gsvd.cc Sat Oct 01 17:02:51 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,465 +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}) -@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 Sat Oct 01 17:02:51 2016 +0100 +++ b/libinterp/corefcn/module.mk Sat Oct 01 17:31:13 2016 +0100 @@ -173,7 +173,6 @@ 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/liboctave/numeric/gsvd.cc Sat Oct 01 17:02:51 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,371 +0,0 @@ -// 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>; - } -}
--- a/liboctave/numeric/gsvd.h Sat Oct 01 17:02:51 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ -// 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 Sat Oct 01 17:02:51 2016 +0100 +++ b/liboctave/numeric/lo-lapack-proto.h Sat Oct 01 17:31:13 2016 +0100 @@ -807,126 +807,6 @@ 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 Sat Oct 01 17:02:51 2016 +0100 +++ b/liboctave/numeric/module.mk Sat Oct 01 17:31:13 2016 +0100 @@ -18,7 +18,6 @@ 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 \ @@ -70,7 +69,6 @@ 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 Sat Oct 01 17:02:51 2016 +0100 +++ b/liboctave/operators/mx-defs.h Sat Oct 01 17:31:13 2016 +0100 @@ -66,8 +66,6 @@ class EIG; -template <typename T> class gsvd; - template <typename T> class hess; template <typename T> class schur;
--- a/liboctave/operators/mx-ext.h Sat Oct 01 17:02:51 2016 +0100 +++ b/liboctave/operators/mx-ext.h Sat Oct 01 17:31:13 2016 +0100 @@ -57,11 +57,6 @@ #include "EIG.h" -// Result of a Generalized Singular Value Decomposition. - -#include "gsvd.h" - - // Result of an LU decomposition. #include "lu.h"