# HG changeset patch # User Barbara Locsi # Date 1470289831 -7200 # Node ID 63b41167ef1efc33b2221df28d88f8d81ed4dcb3 # Parent 66dd260512a441f9ed70c5ab13fc296f4b8153bb gsvd: new function imported from Octave-Forge linear-algebra package. * libinterp/corefcn/gsvd.cc: New function to the interpreter. Imported from the linear-algebra package. * CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, dbleGSVD.h: new classes imported from the linear-algebra package to compute gsvd of Matrix and ComplexMatrix. * liboctave/operators/mx-defs.h, liboctave/operators/mx-ext.h: use new classes. * libinterp/corefcn/module.mk, liboctave/numeric/module.mk: Add to the * scripts/help/__unimplemented__.m: Remove "gsvd" from list. * doc/interpreter/linalg.txi: Add to manual. build system. * NEWS: Add function to list of new functions for 4.2. diff -r 66dd260512a4 -r 63b41167ef1e NEWS --- a/NEWS Tue Aug 09 14:41:52 2016 -0400 +++ b/NEWS Thu Aug 04 07:50:31 2016 +0200 @@ -113,6 +113,7 @@ deg2rad dialog evalc + gsvd hash im2double localfunctions diff -r 66dd260512a4 -r 63b41167ef1e doc/interpreter/linalg.txi --- a/doc/interpreter/linalg.txi Tue Aug 09 14:41:52 2016 -0400 +++ b/doc/interpreter/linalg.txi Thu Aug 04 07:50:31 2016 +0200 @@ -108,6 +108,8 @@ @DOCSTRING(givens) +@DOCSTRING(gsvd) + @DOCSTRING(planerot) @DOCSTRING(inv) diff -r 66dd260512a4 -r 63b41167ef1e libinterp/corefcn/gsvd.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/gsvd.cc Thu Aug 04 07:50:31 2016 +0200 @@ -0,0 +1,500 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006, 2010 Pascal Dupuis +// +// 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 . + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include "defun.h" +#include "error.h" +#include "gripes.h" +#include "pr-output.h" +#include "utils.h" +#include "ovl.h" + +#include "CmplxGSVD.h" +#include "dbleGSVD.h" + +DEFUN (gsvd, args, nargout, + doc: /* -*- texinfo -*- +@deftypefn {Loadable Function} {@var{s} =} gsvd (@var{a}, @var{b}) +@deftypefnx {Loadable Function} {[@var{u}, @var{v}, @var{c}, @var{s}, @var{x} [, @var{r}]] =} gsvd (@var{a}, @var{b}) +@cindex generalised singular value decomposition +Compute the generalised singular value decomposition of (@var{a}, @var{b}): +@iftex +@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 +@end iftex +@ifinfo + +@example +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 example +@end ifinfo + +The function @code{gsvd} normally returns the vector of generalised singular +values +@iftex +@tex +diag(C)./diag(S). +@end tex +@end iftex +@ifinfo +diag(r)./diag(s). +@end ifinfo +If asked for five return values, it computes +@iftex +@tex +$U$, $V$, and $X$. +@end tex +@end iftex +@ifinfo +U, V, and X. +@end ifinfo +With a sixth output argument, it also returns +@iftex +@tex +R, +@end tex +@end iftex +@ifinfo +r, +@end ifinfo +The common upper triangular right term. Other authors, like S. Van Huffel, +define this transformation as the simulatenous diagonalisation of the +input matrices, this can be achieved by multiplying +@iftex +@tex +X +@end tex +@end iftex +@ifinfo +x +@end ifinfo +by the inverse of +@iftex +@tex +[I 0; 0 R]. +@end tex +@end iftex +@ifinfo +[I 0; 0 r]. +@end ifinfo + +For example, + +@example +gsvd (hilb (3), [1 2 3; 3 2 1]) +@end example + +@noindent +returns + +@example +ans = + + 0.1055705 + 0.0031759 +@end example + +@noindent +and + +@example +[u, v, c, s, x, r] = gsvd (hilb (3), [1 2 3; 3 2 1]) +@end example + +@noindent +returns + +@example +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 + +c = + + 0.10499 0.00000 + 0.00000 0.00318 + +s = + 0.99447 0.00000 + 0.00000 0.99999 +x = + + 0.408248 0.902199 0.139179 + -0.816497 0.429063 -0.386314 + 0.408248 -0.044073 -0.911806 + +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 Lapack dggsvd and zggsvd routines. + +@end deftypefn */) +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 2 || nargin > 2 || (nargout > 1 && (nargout < 5 || nargout > 6))) + { + print_usage (); + return retval; + } + + octave_value argA = args(0), argB = args(1); + + octave_idx_type nr = argA.rows (); + octave_idx_type nc = argA.columns (); + +// octave_idx_type nn = argB.rows (); + octave_idx_type np = argB.columns (); + + if (nr == 0 || nc == 0) + { + if (nargout == 5) + retval = ovl (identity_matrix (nc, nc), identity_matrix (nc, nc), + Matrix (nr, nc), identity_matrix (nr, nr), + identity_matrix (nr, nr)); + else if (nargout == 6) + retval = ovl (identity_matrix (nc, nc), identity_matrix (nc, nc), + Matrix (nr, nc), identity_matrix (nr, nr), + identity_matrix (nr, nr), + identity_matrix (nr, nr)); + else + retval = ovl (Matrix (0, 1)); + } + else + { + if ((nc != np)) + { + print_usage (); + return retval; + } + + GSVD::type type = ((nargout == 0 || nargout == 1) + ? GSVD::sigma_only + : (nargout > 5) ? GSVD::std : GSVD::economy ); + + if (argA.is_real_type () && argB.is_real_type ()) + { + Matrix tmpA = argA.matrix_value (); + Matrix tmpB = argB.matrix_value (); + + if (! error_state) + { + if (tmpA.any_element_is_inf_or_nan ()) + { + error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + return retval; + } + + if (tmpB.any_element_is_inf_or_nan ()) + { + error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + return retval; + } + + GSVD result (tmpA, tmpB, type); + + // DiagMatrix sigma = result.singular_values (); + + if (nargout == 0 || nargout == 1) + { + 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 = ovl (sigA.diag()); + } + else + { + if (nargout > 5) + retval = ovl (result.left_singular_matrix_A (), + result.left_singular_matrix_B (), + result.singular_values_A (), + result.singular_values_B (), + result.right_singular_matrix (), + result.R_matrix ()); + else + retval = ovl (result.left_singular_matrix_A (), + result.left_singular_matrix_B (), + result.singular_values_A (), + result.singular_values_B (), + result.right_singular_matrix ()); + } + } + } + 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: cannot take GSVD of matrix containing Inf or NaN values"); + return retval; + } + if (ctmpB.any_element_is_inf_or_nan ()) + { + error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + return retval; + } + + ComplexGSVD result (ctmpA, ctmpB, type); + + // DiagMatrix sigma = result.singular_values (); + + if (nargout == 0 || nargout == 1) + { + 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 = ovl (sigA.diag()); + } + else + { + if (nargout > 5) + retval = ovl (result.left_singular_matrix_A (), + result.left_singular_matrix_B (), + result.singular_values_A (), + result.singular_values_B (), + result.right_singular_matrix (), + result.R_matrix ()); + else + retval = ovl (result.left_singular_matrix_A (), + result.left_singular_matrix_B (), + result.singular_values_A (), + result.singular_values_B (), + result.right_singular_matrix ()); + } + } + } + else + { + gripe_wrong_type_arg ("gsvd", argA); + gripe_wrong_type_arg ("gsvd", argB); + return retval; + } + } + + 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; +%! # disp('Full rank, 5x3 by 3x3 matrices'); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp('A 5x3 full rank, B 3x3 rank deficient'); +%! B(2, 2) = 0; +%! [U, V, C, S, X, 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) + +%! # disp('A 5x3 rank deficient, B 3x3 full rank'); +%! B = B0; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! [U, V, C, S, X, 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) + +%! # disp("A 5x3, B 3x3, [A' B'] rank deficient"); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! [U, V, C, S, X, 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) + +%! # now, A is 3x5 +%! A = A0.'; B0=diag([1 2 4 8 16]); B = B0; +%! # disp('Full rank, 3x5 by 5x5 matrices'); +%! # disp([rank(A) rank(B) rank([A' B'])]); + +%! [U, V, C, S, X, 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) + +%! # disp('A 5x3 full rank, B 5x5 rank deficient'); +%! B(2, 2) = 0; +%! # disp([rank(A) rank(B) rank([A' B'])]); + +%! [U, V, C, S, X, 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) + +%! # disp('A 3x5 rank deficient, B 5x5 full rank'); +%! B = B0; +%! A(3, :) = 2*A(1, :) - A(2, :); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp("A 5x3, B 5x5, [A' B'] rank deficient"); +%! A = A0.'; B = B0.'; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! A0 = A0 +j * randn(5, 3); B0 = B0=diag([1 2 4]) + j*diag([4 -2 -1]); +%! A = A0; B = B0; +%! # disp('Complex: Full rank, 5x3 by 3x3 matrices'); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp('Complex: A 5x3 full rank, B 3x3 rank deficient'); +%! B(2, 2) = 0; +%! # disp([rank(A) rank(B) rank([A' B'])]); + +%! [U, V, C, S, X, 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) + +%! # disp('Complex: A 5x3 rank deficient, B 3x3 full rank'); +%! B = B0; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp("Complex: A 5x3, B 3x3, [A' B'] rank deficient"); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # now, A is 3x5 +%! A = A0.'; B0=diag([1 2 4 8 16])+j*diag([-5 4 -3 2 -1]); +%! B = B0; +%! # disp('Complex: Full rank, 3x5 by 5x5 matrices'); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp('Complex: A 5x3 full rank, B 5x5 rank deficient'); +%! B(2, 2) = 0; +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp('Complex: A 3x5 rank deficient, B 5x5 full rank'); +%! B = B0; +%! A(3, :) = 2*A(1, :) - A(2, :); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + +%! # disp("Complex: A 5x3, B 5x5, [A' B'] rank deficient"); +%! A = A0.'; B = B0.'; +%! A(:, 3) = 2*A(:, 1) - A(:, 2); +%! B(:, 3) = 2*B(:, 1) - B(:, 2); +%! # disp([rank(A) rank(B) rank([A' B'])]); +%! [U, V, C, S, X, 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) + + */ diff -r 66dd260512a4 -r 63b41167ef1e libinterp/corefcn/module.mk --- a/libinterp/corefcn/module.mk Tue Aug 09 14:41:52 2016 -0400 +++ b/libinterp/corefcn/module.mk Thu Aug 04 07:50:31 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 \ diff -r 66dd260512a4 -r 63b41167ef1e liboctave/numeric/CmplxGSVD.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/CmplxGSVD.cc Thu Aug 04 07:50:31 2016 +0200 @@ -0,0 +1,287 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis +// +// 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 . + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include "CmplxGSVD.h" +#include "f77-fcn.h" +#include "lo-error.h" + + +/* + uncomment those lines to monitor k and l + #include "oct-obj.h" + #include "pager.h" +*/ + +extern "C" +{ + F77_RET_T + F77_FUNC (zggsvd, ZGGSVD) + ( + F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 + const octave_idx_type&, // M (input) INTEGER + const octave_idx_type&, // N (input) INTEGER + const octave_idx_type&, // P (input) INTEGER + octave_idx_type &, // K (output) INTEGER + octave_idx_type &, // L (output) INTEGER + Complex*, // A (input/output) COMPLEX*16 array, dimension (LDA,N) + const octave_idx_type&, // LDA (input) INTEGER + Complex*, // B (input/output) COMPLEX*16 array, dimension (LDB,N) + const octave_idx_type&, // LDB (input) INTEGER + double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) + double*, // BETA (output) DOUBLE PRECISION array, dimension (N) + Complex*, // U (output) COMPLEX*16 array, dimension (LDU,M) + const octave_idx_type&, // LDU (input) INTEGER + Complex*, // V (output) COMPLEX*16 array, dimension (LDV,P) + const octave_idx_type&, // LDV (input) INTEGER + Complex*, // Q (output) COMPLEX*16 array, dimension (LDQ,N) + const octave_idx_type&, // LDQ (input) INTEGER + Complex*, // WORK (workspace) COMPLEX*16 array + double*, // RWORK (workspace) DOUBLE PRECISION array + int*, // IWORK (workspace/output) INTEGER array, dimension (N) + octave_idx_type& // INFO (output)INTEGER + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + ); +} + +ComplexMatrix +ComplexGSVD::left_singular_matrix_A (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("CmplxGSVD: U not computed because type == GSVD::sigma_only"); + return ComplexMatrix (); + } + else + return left_smA; +} + +ComplexMatrix +ComplexGSVD::left_singular_matrix_B (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("CmplxGSVD: V not computed because type == GSVD::sigma_only"); + return ComplexMatrix (); + } + else + return left_smB; +} + +ComplexMatrix +ComplexGSVD::right_singular_matrix (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("CmplxGSVD: X not computed because type == GSVD::sigma_only"); + return ComplexMatrix (); + } + else + return right_sm; +} + +ComplexMatrix +ComplexGSVD::R_matrix (void) const +{ + if (type_computed != GSVD::std) + { + (*current_liboctave_error_handler) + ("CmplxGSVD: R not computed because type != GSVD::std"); + return ComplexMatrix (); + } + else + return R; +} + +octave_idx_type +ComplexGSVD::init (const ComplexMatrix& a, const ComplexMatrix& 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 (); + + ComplexMatrix atmp = a; + Complex *tmp_dataA = atmp.fortran_vec (); + + ComplexMatrix btmp = b; + Complex *tmp_dataB = btmp.fortran_vec (); + + // octave_idx_type min_mn = m < n ? m : n; + + 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::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_computed = gsvd_type; + + if (! (jobu == 'N' || jobu == 'O')) { + left_smA.resize (nrow_u, m); + } + + Complex *u = left_smA.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) { + left_smB.resize (nrow_v, p); + } + + Complex *v = left_smB.fortran_vec (); + + if (! (jobq == 'N' || jobq == 'O')) { + right_sm.resize (nrow_q, n); + } + Complex *q = right_sm.fortran_vec (); + + octave_idx_type lwork = 3*n; + lwork = lwork > m ? lwork : m; + lwork = (lwork > p ? lwork : p) + n; + + Array work (dim_vector (lwork, 1)); + Array alpha (dim_vector (n, 1)); + Array beta (dim_vector (n, 1)); + Array rwork(dim_vector (2*n, 1)); + Array iwork (dim_vector (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, tmp_dataA, m, + tmp_dataB, p, alpha.fortran_vec (), + beta.fortran_vec (), u, nrow_u, + v, nrow_v, q, nrow_q, work.fortran_vec (), + rwork.fortran_vec (), iwork.fortran_vec (), + info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zggsvd"); + + if (info < 0) { + (*current_liboctave_error_handler) ("zggsvd.f: argument %d illegal", -info); + } else { + if (info > 0) { + (*current_liboctave_error_handler) ("zggsvd.f: Jacobi-type procedure failed to converge."); + } else { + octave_idx_type i, j; + + if (GSVD::std == gsvd_type) { + R.resize(k+l, k+l); + int astart = n-k-l; + if (m - k - l >= 0) { + int 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); + } + } + } + /* + uncomment this to monitor k and l + octave_value tmp; + octave_stdout << "CmplxGSVD k: "; + tmp = k; + tmp.print(octave_stdout); + octave_stdout << "\n"; + octave_stdout << "CmplxGSVD l: "; + tmp = l; + tmp.print(octave_stdout); + octave_stdout << "\n"; + */ + + 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); + } + } + } + } + return info; +} diff -r 66dd260512a4 -r 63b41167ef1e liboctave/numeric/CmplxGSVD.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/CmplxGSVD.h Thu Aug 04 07:50:31 2016 +0200 @@ -0,0 +1,94 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis +// +// 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 . + +#if !defined (octave_ComplexGSVD_h) +#define octave_ComplexGSVD_h 1 + +#include "octave-config.h" + +#include + +#include "dDiagMatrix.h" +#include "CMatrix.h" +#include "dbleGSVD.h" + +class +ComplexGSVD +{ +public: + + ComplexGSVD (void) { } + + ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, + GSVD::type gsvd_type = GSVD::economy) + { + init (a, b, gsvd_type); + } + + ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, + octave_idx_type& info, GSVD::type gsvd_type = GSVD::economy) + { + info = init (a, b, gsvd_type); + } + + ComplexGSVD (const ComplexGSVD& a) + : type_computed (a.type_computed), + sigmaA (a.sigmaA), sigmaB (a.sigmaB), + left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm), + R(a.R) { } + + ComplexGSVD& operator = (const ComplexGSVD& a) + { + if (this != &a) + { + type_computed = a.type_computed; + 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; + } + + ~ComplexGSVD (void) { } + + DiagMatrix singular_values_A (void) const { return sigmaA; } + DiagMatrix singular_values_B (void) const { return sigmaB; } + + ComplexMatrix left_singular_matrix_A (void) const; + ComplexMatrix left_singular_matrix_B (void) const; + + ComplexMatrix right_singular_matrix (void) const; + ComplexMatrix R_matrix (void) const; + + friend std::ostream& operator << (std::ostream& os, const ComplexGSVD& a); + +private: + + GSVD::type type_computed; + + DiagMatrix sigmaA, sigmaB; + ComplexMatrix left_smA, left_smB; + ComplexMatrix right_sm, R; + + octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, + GSVD::type gsvd_type = GSVD::economy); +}; + +#endif diff -r 66dd260512a4 -r 63b41167ef1e liboctave/numeric/dbleGSVD.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/dbleGSVD.cc Thu Aug 04 07:50:31 2016 +0200 @@ -0,0 +1,293 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis +// +// 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 . + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include "dbleGSVD.h" +#include "f77-fcn.h" +#include "lo-error.h" + +/* + uncomment those lines to monitor k and l + #include "oct-obj.h" + #include "pager.h" +*/ + +extern "C" +{ + F77_RET_T + F77_FUNC (dggsvd, DGGSVD) + ( + F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 + const octave_idx_type&, // M (input) INTEGER + const octave_idx_type&, // N (input) INTEGER + const octave_idx_type&, // P (input) INTEGER + octave_idx_type &, // K (output) INTEGER + octave_idx_type &, // L (output) INTEGER + double*, // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + const octave_idx_type&, // LDA (input) INTEGER + double*, // B (input/output) DOUBLE PRECISION array, dimension (LDB,N) + const octave_idx_type&, // LDB (input) INTEGER + double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) + double*, // BETA (output) DOUBLE PRECISION array, dimension (N) + double*, // U (output) DOUBLE PRECISION array, dimension (LDU,M) + const octave_idx_type&, // LDU (input) INTEGER + double*, // V (output) DOUBLE PRECISION array, dimension (LDV,P) + const octave_idx_type&, // LDV (input) INTEGER + double*, // Q (output) DOUBLE PRECISION array, dimension (LDQ,N) + const octave_idx_type&, // LDQ (input) INTEGER + double*, // WORK (workspace) DOUBLE PRECISION array + int*, // IWORK (workspace/output) INTEGER array, dimension (N) + octave_idx_type& // INFO (output)INTEGER + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + ); +} + +Matrix +GSVD::left_singular_matrix_A (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("dbleGSVD: U not computed because type == GSVD::sigma_only"); + return Matrix (); + } + else + return left_smA; +} + +Matrix +GSVD::left_singular_matrix_B (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("dbleGSVD: V not computed because type == GSVD::sigma_only"); + return Matrix (); + } + else + return left_smB; +} + +Matrix +GSVD::right_singular_matrix (void) const +{ + if (type_computed == GSVD::sigma_only) + { + (*current_liboctave_error_handler) + ("dbleGSVD: X not computed because type == GSVD::sigma_only"); + return Matrix (); + } + else + return right_sm; +} +Matrix +GSVD::R_matrix (void) const +{ + if (type_computed != GSVD::std) + { + (*current_liboctave_error_handler) + ("dbleGSVD: R not computed because type != GSVD::std"); + return Matrix (); + } + else + return R; +} + +octave_idx_type +GSVD::init (const Matrix& a, const Matrix& 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 (); + + Matrix atmp = a; + double *tmp_dataA = atmp.fortran_vec (); + + Matrix btmp = b; + double *tmp_dataB = btmp.fortran_vec (); + + // octave_idx_type min_mn = m < n ? m : n; + + 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::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_computed = gsvd_type; + + if (! (jobu == 'N' || jobu == 'O')) { + left_smA.resize (nrow_u, m); + } + + double *u = left_smA.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) { + left_smB.resize (nrow_v, p); + } + + double *v = left_smB.fortran_vec (); + + if (! (jobq == 'N' || jobq == 'O')) { + right_sm.resize (nrow_q, n); + } + double *q = right_sm.fortran_vec (); + + octave_idx_type lwork = 3*n; + lwork = lwork > m ? lwork : m; + lwork = (lwork > p ? lwork : p) + n; + + Array work (dim_vector (lwork, 1)); + Array alpha (dim_vector (n, 1)); + Array beta (dim_vector (n, 1)); + Array iwork (dim_vector (n, 1)); + + 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, m, + tmp_dataB, p, alpha.fortran_vec (), + beta.fortran_vec (), u, nrow_u, + v, nrow_v, q, nrow_q, work.fortran_vec (), + iwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dggsvd"); + + if (info < 0) { + (*current_liboctave_error_handler) ("dggsvd.f: argument %d illegal", -info); + } else { + if (info > 0) { + (*current_liboctave_error_handler) ("dggsvd.f: Jacobi-type procedure failed to converge."); + } else { + octave_idx_type i, j; + + if (GSVD::std == gsvd_type) { + R.resize(k+l, k+l); + int astart = n-k-l; + if (m - k - l >= 0) { + int 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); + } + } + } + /* + uncomment this to monitor k and l + octave_value tmp; + octave_stdout << "dbleGSVD k: "; + tmp = k; + tmp.print(octave_stdout); + octave_stdout << "\n"; + octave_stdout << "dbleGSVD l: "; + tmp = l; + tmp.print(octave_stdout); + octave_stdout << "\n"; + */ + + 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); + } + } + } + } + return info; +} + +std::ostream& +operator << (std::ostream& os, const GSVD& a) +{ + os << a.left_singular_matrix_A () << "\n"; + os << a.left_singular_matrix_B () << "\n"; + os << a.singular_values_A () << "\n"; + os << a.singular_values_B () << "\n"; + os << a.right_singular_matrix () << "\n"; + + return os; +} diff -r 66dd260512a4 -r 63b41167ef1e liboctave/numeric/dbleGSVD.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/dbleGSVD.h Thu Aug 04 07:50:31 2016 +0200 @@ -0,0 +1,91 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis +// +// 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 . + +#if !defined (octave_GSVD_h) +#define octave_GSVD_h 1 + +#include "octave-config.h" + +#include "dDiagMatrix.h" +#include "dMatrix.h" + +class +GSVD +{ +public: + + enum type + { + std, + economy, + sigma_only + }; + + GSVD (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { } + GSVD (const Matrix& a, const Matrix& b, type gsvd_type = GSVD::economy) { init (a, b, gsvd_type); } + + GSVD (const Matrix& a, const Matrix& b, octave_idx_type& info, type gsvd_type = GSVD::economy) + { + info = init (a, b, gsvd_type); + } + + GSVD (const GSVD& a) + : type_computed (a.type_computed), + 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_computed = a.type_computed; + 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) { } + + DiagMatrix singular_values_A (void) const { return sigmaA; } + DiagMatrix singular_values_B (void) const { return sigmaB; } + + Matrix left_singular_matrix_A (void) const; + Matrix left_singular_matrix_B (void) const; + + Matrix right_singular_matrix (void) const; + Matrix R_matrix (void) const; + + friend std::ostream& operator << (std::ostream& os, const GSVD& a); + +private: + + GSVD::type type_computed; + + DiagMatrix sigmaA, sigmaB; + Matrix left_smA, left_smB; + Matrix right_sm, R; + + octave_idx_type init (const Matrix& a, const Matrix& b, type gsvd_type = economy); +}; + +#endif diff -r 66dd260512a4 -r 63b41167ef1e liboctave/numeric/module.mk --- a/liboctave/numeric/module.mk Tue Aug 09 14:41:52 2016 -0400 +++ b/liboctave/numeric/module.mk Thu Aug 04 07:50:31 2016 +0200 @@ -8,6 +8,7 @@ LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in) NUMERIC_INC = \ + liboctave/numeric/CmplxGSVD.h \ liboctave/numeric/CollocWt.h \ liboctave/numeric/DAE.h \ liboctave/numeric/DAEFunc.h \ @@ -16,6 +17,7 @@ liboctave/numeric/DASPK.h \ liboctave/numeric/DASRT.h \ liboctave/numeric/DASSL.h \ + liboctave/numeric/dbleGSVD.h \ liboctave/numeric/DET.h \ liboctave/numeric/EIG.h \ liboctave/numeric/LSODE.h \ @@ -56,10 +58,12 @@ liboctave/numeric/svd.h NUMERIC_SRC = \ + liboctave/numeric/CmplxGSVD.cc \ liboctave/numeric/CollocWt.cc \ liboctave/numeric/DASPK.cc \ liboctave/numeric/DASRT.cc \ liboctave/numeric/DASSL.cc \ + liboctave/numeric/dbleGSVD.cc \ liboctave/numeric/EIG.cc \ liboctave/numeric/LSODE.cc \ liboctave/numeric/ODES.cc \ diff -r 66dd260512a4 -r 63b41167ef1e liboctave/operators/mx-defs.h --- a/liboctave/operators/mx-defs.h Tue Aug 09 14:41:52 2016 -0400 +++ b/liboctave/operators/mx-defs.h Thu Aug 04 07:50:31 2016 +0200 @@ -66,6 +66,8 @@ class EIG; +class GSVD; + template class hess; template class schur; diff -r 66dd260512a4 -r 63b41167ef1e liboctave/operators/mx-ext.h --- a/liboctave/operators/mx-ext.h Tue Aug 09 14:41:52 2016 -0400 +++ b/liboctave/operators/mx-ext.h Thu Aug 04 07:50:31 2016 +0200 @@ -57,6 +57,12 @@ #include "EIG.h" +// Result of a Generalized Singular Value Decomposition. + +#include "dbleGSVD.h" +#include "CmplxGSVD.h" + + // Result of an LU decomposition. #include "lu.h" diff -r 66dd260512a4 -r 63b41167ef1e scripts/help/__unimplemented__.m --- a/scripts/help/__unimplemented__.m Tue Aug 09 14:41:52 2016 -0400 +++ b/scripts/help/__unimplemented__.m Thu Aug 04 07:50:31 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 ", ... @@ -658,7 +653,6 @@ "graph", "graymon", "griddedInterpolant", - "gsvd", "guide", "h5create", "h5disp",