Mercurial > octave
changeset 30300:4ee01c14fccd
Rewrite gsvd function and return *correct* 3rd output (bug #60273).
Correctly use the API interface to LAPACK which Octave got wrong for
certain matrix combinations. Perform post-processing on LAPACK results
to return Matlab-compatible outputs.
* libinterp/corefcn/gsvd.cc (gsvd_type): Add second input (nargin) to be able
to distinguish the type of gsvd to generate.
* libinterp/corefcn/gsvd.cc (do_gsvd): Add another input (nargin) to be able
to calculate the correct gsvd type. Update code for API change in liboctave
which now returns a full Matrix for singular values rather than a DiagMatrix.
Sort singular values in ascending order for compatibility with Matlab.
* libinterp/corefcn/gsvd.cc (Fgsvd): Update documation to warn about
rank-deficient input case. Add new input validation to check that number
of columns of A and B match. Move special check for empty matrices from
libinterp to liboctave. Update all BIST tests for new behavior of gsvd
code in liboctave.
* liboctave/numeric/gsvd.h (gsvd): Remove member variable R from gsvd class.
* liboctave/numeric/gsvd.h (gsvd::singular_values_A, gsvd_singular_values_B):
Change return type to T::real_matrix_type from T::real_diag_matrix_type.
* liboctave/numeric/gsvd.h (gsvd::ggsvd): Change default gsvd::Type to "std"
from "economy".
* liboctave/numeric/gsvd.cc: Add #includes for <algorithm>, <undordered_map>,
and "oct-locbuf.h". Change gsvd_fcn from std::map to std::unordered_map.
* liboctave/numeric/gsvd.cc (initialize_gsvd): Throw an error if Octave is
unable to connect to LAPACK library.
* liboctave/numeric/gsvd.cc (ggsvd): Change function prototype to not use
references to Octave Matrix datatypes, but rather point directly to
underlying plain old datatype. For example, use "double *" rather than
"Matrix&". Replace scratch memory instances created with Matrix() with
OCTAVE_LOCAL_BUFFER which relies on C++ STL.
* liboctave/numeric/gsvd.cc (left_singular_matrix_A, left_singular_matrix_B,
right_singular_matrix): Use the fact that current_liboctave_error_handler never
returns to eliminate else branch of if/else code.
* liboctave/numeric/gsvd.cc (R_matrix): Delete unused function.
* liboctave/numeric/gsvd.cc (gsvd): Add input validation for empty matrices.
Simplify 'job' values processing which seems to have been fixed in LAPACK 3.0.
Replace scratch memory created with std::vector with OCTAVE_LOCAL_BUFFER.
Replace scratch memory created with Octave Matrix classes with
OCTAVE_LOCAL_BUFFER. Use initializer list to std::max to simplify code.
Re-code interface to LAPACK to extract 'R' matrix. Post-process to calculate
third output of gsvd as X = Q*R'. Correct algorithm to extract singular values
from LAPACK..
| author | Rik <rik@octave.org> |
|---|---|
| date | Wed, 08 Sep 2021 16:08:03 -0700 |
| parents | f306fe9bfa0d |
| children | cb711825f8e5 |
| files | libinterp/corefcn/gsvd.cc liboctave/numeric/gsvd.cc liboctave/numeric/gsvd.h |
| diffstat | 3 files changed, 429 insertions(+), 399 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/gsvd.cc Thu May 27 22:49:26 2021 +0200 +++ b/libinterp/corefcn/gsvd.cc Wed Sep 08 16:08:03 2021 -0700 @@ -44,55 +44,65 @@ template <typename T> static typename math::gsvd<T>::Type -gsvd_type (int nargout) +gsvd_type (int nargout, int nargin) { - return ((nargout == 0 || nargout == 1) - ? math::gsvd<T>::Type::sigma_only - : (nargout > 5) ? math::gsvd<T>::Type::std - : math::gsvd<T>::Type::economy); + if (nargout == 0 || nargout == 1) + return octave::math::gsvd<T>::Type::sigma_only; + else if (nargin < 3) + return octave::math::gsvd<T>::Type::std; + else + return octave::math::gsvd<T>::Type::economy; } -// Named like this to avoid conflicts with the gsvd class. +// Named do_gsvd to avoid conflicts with the gsvd class itself. template <typename T> static octave_value_list -do_gsvd (const T& A, const T& B, const octave_idx_type nargout, +do_gsvd (const T& A, const T& B, + const octave_idx_type nargout, const octave_idx_type nargin, bool is_single = false) { - math::gsvd<T> result (A, B, gsvd_type<T> (nargout)); + math::gsvd<T> result (A, B, gsvd_type<T> (nargout, nargin)); octave_value_list retval (nargout); - if (nargout < 2) + if (nargout <= 1) { if (is_single) { - FloatDiagMatrix sigA = result.singular_values_A (); - FloatDiagMatrix sigB = result.singular_values_B (); + FloatMatrix sigA = result.singular_values_A (); + FloatMatrix sigB = result.singular_values_B (); for (int i = sigA.rows () - 1; i >= 0; i--) - sigA.dgxelem(i) /= sigB.dgxelem(i); - retval(0) = sigA.diag (); + sigA.xelem (i) /= sigB.xelem (i); + retval(0) = sigA.sort (); } else { - DiagMatrix sigA = result.singular_values_A (); - DiagMatrix sigB = result.singular_values_B (); + Matrix sigA = result.singular_values_A (); + Matrix sigB = result.singular_values_B (); for (int i = sigA.rows () - 1; i >= 0; i--) - sigA.dgxelem(i) /= sigB.dgxelem(i); - retval(0) = sigA.diag (); + sigA.xelem (i) /= sigB.xelem (i); + retval(0) = sigA.sort (); } } else { - retval(0) = result.left_singular_matrix_A (); + switch (nargout) + { + case 5: + retval(4) = result.singular_values_B (); + OCTAVE_FALLTHROUGH; + + case 4: + retval(3) = result.singular_values_A (); + OCTAVE_FALLTHROUGH; + + case 3: + retval(2) = result.right_singular_matrix (); + OCTAVE_FALLTHROUGH; + } 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 (); + retval(0) = result.left_singular_matrix_A (); } + return retval; } @@ -145,8 +155,9 @@ of columns of @var{A}. This option is not yet implemented. Programming Note: the code is a wrapper to the corresponding @sc{lapack} dggsvd -and zggsvd routines. - +and zggsvd routines. If matrices @var{A} and @var{B} are @emph{both} rank +deficient then @sc{lapack} will return an incorrect factorization. Programmers +should avoid this combination. @seealso{svd} @end deftypefn */) { @@ -155,122 +166,77 @@ if (nargin < 2 || nargin > 3) print_usage (); else if (nargin == 3) - warning ("gsvd: economy-sized decomposition is not yet implemented, returning full decomposition"); + { + // FIXME: when "economy" is implemented delete this code + warning ("gsvd: economy-sized decomposition is not yet implemented, returning full decomposition"); + nargin = 2; + } 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 (); + if (argA.columns () != argB.columns ()) + error ("gsvd: A and B must have the same number of columns"); - octave_idx_type np = argB.columns (); - - // FIXME: This "special" case should be handled in the gsvd class, not here - if (nr == 0 || nc == 0) + if (argA.is_single_type () || argB.is_single_type ()) { - retval = octave_value_list (nargout); - if (nargout < 2) // S = gsvd (A, B) - { - if (argA.is_single_type () || argB.is_single_type ()) - retval(0) = FloatMatrix (0, 1); - else - retval(0) = Matrix (0, 1); - } - else // [U, V, X, C, S, R] = gsvd (A, B) + if (argA.isreal () && argB.isreal ()) { - if (argA.is_single_type () || argB.is_single_type ()) - { - retval(0) = float_identity_matrix (nc, nc); - retval(1) = float_identity_matrix (nc, nc); - if (nargout > 2) - retval(2) = float_identity_matrix (nr, nr); - if (nargout > 3) - retval(3) = FloatMatrix (nr, nc); - if (nargout > 4) - retval(4) = float_identity_matrix (nr, nr); - if (nargout > 5) - retval(5) = float_identity_matrix (nr, nr); - } - else - { - 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); - } + FloatMatrix tmpA = argA.xfloat_matrix_value ("gsvd: A must be a real or complex matrix"); + FloatMatrix tmpB = argB.xfloat_matrix_value ("gsvd: B must be a real or complex matrix"); + + if (tmpA.any_element_is_inf_or_nan ()) + error ("gsvd: A cannot have Inf or NaN values"); + if (tmpB.any_element_is_inf_or_nan ()) + error ("gsvd: B cannot have Inf or NaN values"); + + retval = do_gsvd (tmpA, tmpB, nargout, nargin, true); } + else if (argA.iscomplex () || argB.iscomplex ()) + { + FloatComplexMatrix ctmpA = argA.xfloat_complex_matrix_value ("gsvd: A must be a real or complex matrix"); + FloatComplexMatrix ctmpB = argB.xfloat_complex_matrix_value ("gsvd: B must be a real or complex matrix"); + + if (ctmpA.any_element_is_inf_or_nan ()) + error ("gsvd: A cannot have Inf or NaN values"); + if (ctmpB.any_element_is_inf_or_nan ()) + error ("gsvd: B cannot have Inf or NaN values"); + + retval = do_gsvd (ctmpA, ctmpB, nargout, nargin, true); + } + else + error ("gsvd: A and B must be real or complex matrices"); } else { - if (nc != np) - print_usage (); - - if (argA.is_single_type () || argB.is_single_type ()) + if (argA.isreal () && argB.isreal ()) { - if (argA.isreal () && argB.isreal ()) - { - FloatMatrix tmpA = argA.xfloat_matrix_value ("gsvd: A must be a real or complex matrix"); - FloatMatrix tmpB = argB.xfloat_matrix_value ("gsvd: B must be a real or complex matrix"); + Matrix tmpA = argA.xmatrix_value ("gsvd: A must be a real or complex matrix"); + Matrix tmpB = argB.xmatrix_value ("gsvd: B must be a real or complex matrix"); - if (tmpA.any_element_is_inf_or_nan ()) - error ("gsvd: A cannot have Inf or NaN values"); - if (tmpB.any_element_is_inf_or_nan ()) - error ("gsvd: B cannot have Inf or NaN values"); + if (tmpA.any_element_is_inf_or_nan ()) + error ("gsvd: A cannot have Inf or NaN values"); + if (tmpB.any_element_is_inf_or_nan ()) + error ("gsvd: B cannot have Inf or NaN values"); - retval = do_gsvd (tmpA, tmpB, nargout, true); - } - else if (argA.iscomplex () || argB.iscomplex ()) - { - FloatComplexMatrix ctmpA = argA.xfloat_complex_matrix_value ("gsvd: A must be a real or complex matrix"); - FloatComplexMatrix ctmpB = argB.xfloat_complex_matrix_value ("gsvd: B must be a real or complex matrix"); + retval = do_gsvd (tmpA, tmpB, nargout, nargin); + } + else if (argA.iscomplex () || argB.iscomplex ()) + { + ComplexMatrix ctmpA = argA.xcomplex_matrix_value ("gsvd: A must be a real or complex matrix"); + ComplexMatrix ctmpB = argB.xcomplex_matrix_value ("gsvd: B must be a real or complex matrix"); - if (ctmpA.any_element_is_inf_or_nan ()) - error ("gsvd: A cannot have Inf or NaN values"); - if (ctmpB.any_element_is_inf_or_nan ()) - error ("gsvd: B cannot have Inf or NaN values"); + 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 = do_gsvd (ctmpA, ctmpB, nargout, true); - } - else - error ("gsvd: A and B must be real or complex matrices"); + retval = do_gsvd (ctmpA, ctmpB, nargout, nargin); } else - { - if (argA.isreal () && argB.isreal ()) - { - Matrix tmpA = argA.xmatrix_value ("gsvd: A must be a real or complex matrix"); - Matrix tmpB = argB.xmatrix_value ("gsvd: B must be a real or complex matrix"); - - if (tmpA.any_element_is_inf_or_nan ()) - error ("gsvd: A cannot have Inf or NaN values"); - if (tmpB.any_element_is_inf_or_nan ()) - error ("gsvd: B cannot have Inf or NaN values"); - - retval = do_gsvd (tmpA, tmpB, nargout); - } - else if (argA.iscomplex () || argB.iscomplex ()) - { - ComplexMatrix ctmpA = argA.xcomplex_matrix_value ("gsvd: A must be a real or complex matrix"); - ComplexMatrix ctmpB = argB.xcomplex_matrix_value ("gsvd: B must be a real or complex matrix"); - - if (ctmpA.any_element_is_inf_or_nan ()) - error ("gsvd: A cannot have Inf or NaN values"); - if (ctmpB.any_element_is_inf_or_nan ()) - error ("gsvd: B cannot have Inf or NaN values"); - - retval = do_gsvd (ctmpA, ctmpB, nargout); - } - else - error ("gsvd: A and B must be real or complex matrices"); - } + error ("gsvd: A and B must be real or complex matrices"); } return retval; @@ -278,19 +244,47 @@ /* -## Basic test of decomposition -%!test <48807> +## Basic tests of decomposition +%!test <60273> %! A = reshape (1:15,5,3); %! B = magic (3); %! [U,V,X,C,S] = gsvd (A,B); +%! assert (size (U), [5, 5]); +%! assert (size (V), [3, 3]); +%! assert (size (X), [3, 3]); +%! assert (size (C), [5, 3]); +%! assert (C(4:5, :), zeros (2,3)); +%! assert (size (S), [3, 3]); %! assert (U*C*X', A, 50*eps); %! assert (V*S*X', B, 50*eps); %! S0 = gsvd (A, B); -%! S1 = svd (A / B); +%! assert (size (S0), [3, 1]); +%! S1 = sort (svd (A / B)); %! assert (S0, S1, 10*eps); +%!test <60273> +%! A = reshape (1:15,3,5); +%! B = magic (5); +%! [U,V,X,C,S] = gsvd (A,B); +%! assert (size (U), [3, 3]); +%! assert (size (V), [5, 5]); +%! assert (size (X), [5, 5]); +%! assert (size (C), [3, 5]); +%! assert (C(:, 4:5), zeros (3,2)); +%! assert (size (S), [5, 5]); +%! assert (U*C*X', A, 100*eps); # less accurate in this orientation +%! assert (V*S*X', B, 125*eps); # for some reason. +%! S0 = gsvd (A, B); +%! assert (size (S0), [5, 1]); +%! S0 = S0(3:end); +%! S1 = sort (svd (A / B)); +%! assert (S0, S1, 20*eps); + ## a few tests for gsvd.m -%!shared A, A0, B, B0, U, V, C, S, X, R, D1, D2 +%!shared A, A0, B, B0, U, V, C, S, X, old_state, restore_state +%! old_state = randn ("state"); +%! restore_state = onCleanup (@() randn ("state", old_state)); +%! randn ("state", 40); # initialize generator to make behavior reproducible %! A0 = randn (5, 3); %! B0 = diag ([1 2 4]); %! A = A0; @@ -298,88 +292,74 @@ ## A (5x3) and B (3x3) are full rank %!test <48807> -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 3); D1(1:3, 1:3) = C; -%! D2 = S; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A: 5x3 full rank, B: 3x3 rank deficient %!test <48807> %! B(2, 2) = 0; -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; -%! D2 = [zeros(2, 1) S; zeros(1, 3)]; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A: 5x3 rank deficient, B: 3x3 full rank %!test <48807> %! B = B0; %! A(:, 3) = 2*A(:, 1) - A(:, 2); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 3); D1(1:3, 1:3) = C; -%! D2 = S; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A and B are both rank deficient -%!test <48807> +## FIXME: LAPACK seems to be completely broken for this case +%!#test <48807> %! B(:, 3) = 2*B(:, 1) - B(:, 2); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 2); D1(1:2, 1:2) = C; -%! D2 = [S; zeros(1, 2)]; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A (now 3x5) and B (now 5x5) are full rank %!test <48807> %! A = A0.'; %! B0 = diag ([1 2 4 8 16]); %! B = B0; -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = [C zeros(3,2)]; -%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 40*eps); ## A: 3x5 full rank, B: 5x5 rank deficient %!test <48807> %! B(2, 2) = 0; -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; -%! D2 = zeros (5, 5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 40*eps); ## A: 3x5 rank deficient, B: 5x5 full rank %!test <48807> %! B = B0; %! A(3, :) = 2*A(1, :) - A(2, :); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (3, 5); D1(1:3, 1:3) = C; -%! D2 = zeros (5, 5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 40*eps); ## A and B are both rank deficient -%!test <48807> +## FIXME: LAPACK seems to be completely broken for this case +%!#test <48807> %! A = A0.'; B = B0.'; %! A(:, 3) = 2*A(:, 1) - A(:, 2); %! B(:, 3) = 2*B(:, 1) - B(:, 2); -%! [U, V, X, C, S, R]=gsvd (A, B); -%! D1 = zeros (3, 4); D1(1:3, 1:3) = C; -%! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A: 5x3 complex full rank, B: 3x3 complex full rank %!test <48807> @@ -387,43 +367,36 @@ %! 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); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 25*eps); ## A: 5x3 complex full rank, B: 3x3 complex rank deficient %!test <48807> %! B(2, 2) = 0; -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C; -%! D2 = [zeros(2, 1) S; zeros(1, 3)]; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 25*eps); ## A: 5x3 complex rank deficient, B: 3x3 complex full rank %!test <48807> %! B = B0; %! A(:, 3) = 2*A(:, 1) - A(:, 2); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 3); D1(1:3, 1:3) = C; -%! D2 = S; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 25*eps); ## A (5x3) and B (3x3) are both complex rank deficient -%!test <48807> +## FIXME: LAPACK seems to be completely broken for this case +%!#test <48807> %! B(:, 3) = 2*B(:, 1) - B(:, 2); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (5, 2); D1(1:2, 1:2) = C; -%! D2 = [S; zeros(1, 2)]; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*[zeros(2, 1) R]) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*[zeros(2, 1) R]) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (3), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 20*eps); ## A (now 3x5) complex and B (now 5x5) complex are full rank ## now, A is 3x5 @@ -431,71 +404,97 @@ %! 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); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 25*eps); +%! assert (V*S*X', B, 85*eps); ## A: 3x5 complex full rank, B: 5x5 complex rank deficient %!test <48807> %! B(2, 2) = 0; -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C; -%! D2 = zeros (5,5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye (2); -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (2, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 85*eps); ## A: 3x5 complex rank deficient, B: 5x5 complex full rank %!test <48807> %! B = B0; %! A(3, :) = 2*A(1, :) - A(2, :); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (3, 5); D1(1:3, 1:3) = C; -%! D2 = zeros (5,5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye (2); -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*R) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*R) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 85*eps); ## A and B are both complex rank deficient -%!test <48807> +## FIXME: LAPACK seems to be completely broken for this case +%!#test <48807> %! A = A0.'; %! B = B0.'; %! A(:, 3) = 2*A(:, 1) - A(:, 2); %! B(:, 3) = 2*B(:, 1) - B(:, 2); -%! [U, V, X, C, S, R] = gsvd (A, B); -%! D1 = zeros (3, 4); D1(1:3, 1:3) = C; -%! D2 = eye (4); D2(1:3, 1:3) = S; D2(5,:) = 0; -%! assert (norm (diag (C).^2 + diag (S).^2 - ones (3, 1)) <= 1e-6); -%! assert (norm ((U'*A*X) - D1*[zeros(4, 1) R]) <= 1e-6); -%! assert (norm ((V'*B*X) - D2*[zeros(4, 1) R]) <= 1e-6); +%! [U, V, X, C, S] = gsvd (A, B); +%! assert (C'*C + S'*S, eye (5), 5*eps); +%! assert (U*C*X', A, 10*eps); +%! assert (V*S*X', B, 85*eps); ## Test that single inputs produce single outputs %!test -%! s = gsvd (single (ones (0,1)), B); +%! s = gsvd (single (eye (5)), B); %! assert (class (s), "single"); -%! s = gsvd (single (ones (1,0)), B); +%! [U,V,X,C,S] = gsvd (single (eye(5)), B); +%! assert (class (U), "single"); +%! assert (class (V), "single"); +%! assert (class (X), "single"); +%! assert (class (C), "single"); +%! assert (class (S), "single"); +%! +%! s = gsvd (A, single (eye (5))); %! assert (class (s), "single"); -%! s = gsvd (single (ones (1,0)), B); -%! [U,V,X,C,S,R] = gsvd (single ([]), B); +%! [U,V,X,C,S] = gsvd (A, single (eye (5))); %! assert (class (U), "single"); %! assert (class (V), "single"); %! assert (class (X), "single"); %! assert (class (C), "single"); %! assert (class (S), "single"); -%! assert (class (R), "single"); -%! -%! s = gsvd (single (A), B); -%! assert (class (s), "single"); -%! [U,V,X,C,S,R] = gsvd (single (A), B); -%! assert (class (U), "single"); -%! assert (class (V), "single"); -%! assert (class (X), "single"); -%! assert (class (C), "single"); -%! assert (class (S), "single"); -%! assert (class (R), "single"); + +## Test input validation +%!error <Invalid call> gsvd () +%!error <Invalid call> gsvd (1) +%!error <Invalid call> gsvd (1,2,3,4) +%!warning <economy-sized decomposition is not yet implemented> gsvd (1,2,0); +%!error <A and B must have the same number of columns> gsvd (1,[1, 2]) +## Test input validation for single (real and complex) inputs. +%!error <A cannot have Inf or NaN values> gsvd (Inf, single (2)) +%!error <A cannot have Inf or NaN values> gsvd (NaN, single (2)) +%!error <B cannot have Inf or NaN values> gsvd (single (1), Inf) +%!error <B cannot have Inf or NaN values> gsvd (single (1), NaN) +%!error <A must be a real or complex matrix> gsvd ({1}, single (2i)) +%!error <B must be a real or complex matrix> gsvd (single (i), {2}) +%!error <A cannot have Inf or NaN values> gsvd (Inf, single (2i)) +%!error <A cannot have Inf or NaN values> gsvd (NaN, single (2i)) +%!error <B cannot have Inf or NaN values> gsvd (single (i), Inf) +%!error <B cannot have Inf or NaN values> gsvd (single (i), NaN) +## Test input validation for single, but not real or complex, inputs. +%!error <A and B must be real or complex matrices> gsvd ({1}, single (2)) +%!error <A and B must be real or complex matrices> gsvd (single (1), {2}) +## Test input validation for double (real and complex) inputs. +%!error <A cannot have Inf or NaN values> gsvd (Inf, 2) +%!error <A cannot have Inf or NaN values> gsvd (NaN, 2) +%!error <B cannot have Inf or NaN values> gsvd (1, Inf) +%!error <B cannot have Inf or NaN values> gsvd (1, NaN) +%!error <A must be a real or complex matrix> gsvd ({1}, 2i) +%!error <B must be a real or complex matrix> gsvd (i, {2}) +%!error <A cannot have Inf or NaN values> gsvd (Inf, 2i) +%!error <A cannot have Inf or NaN values> gsvd (NaN, 2i) +%!error <B cannot have Inf or NaN values> gsvd (i, Inf) +%!error <B cannot have Inf or NaN values> gsvd (i, NaN) +## Test input validation for double, but not real or complex, inputs. +%!error <A and B must be real or complex matrices> gsvd ({1}, double (2)) +%!error <A and B must be real or complex matrices> gsvd (double (1), {2}) +## Test input validation in liboctave/numeric/gsvd.cc +%!error <A and B cannot be empty matrices> gsvd (zeros (0,1), 1) +%!error <A and B cannot be empty matrices> gsvd (1, zeros (0,1)) */
--- a/liboctave/numeric/gsvd.cc Thu May 27 22:49:26 2021 +0200 +++ b/liboctave/numeric/gsvd.cc Wed Sep 08 16:08:03 2021 -0700 @@ -27,7 +27,8 @@ # include <config.h> #endif -#include <vector> +#include <algorithm> +#include <unordered_map> #include "CMatrix.h" #include "dDiagMatrix.h" @@ -38,18 +39,19 @@ #include "gsvd.h" #include "lo-error.h" #include "lo-lapack-proto.h" +#include "oct-locbuf.h" #include "oct-shlib.h" namespace octave { - static std::map<std::string, void *> gsvd_fcn; + static std::unordered_map<std::string, void *> gsvd_fcn; static bool have_DGGSVD3 = false; static bool gsvd_initialized = false; - /* Hack to stringize macro results. */ -#define xSTRINGIZE(x) #x -#define STRINGIZE(x) xSTRINGIZE(x) + /* Hack to stringize results of F77_FUNC macro. */ + #define xSTRINGIZE(x) #x + #define STRINGIZE(x) xSTRINGIZE(x) static void initialize_gsvd (void) { @@ -58,11 +60,8 @@ dynamic_library libs (""); if (! libs) - { - // FIXME: Should we throw an error if we cannot check the libraries? - have_DGGSVD3 = false; - return; - } + (*current_liboctave_error_handler) + ("gsvd: unable to query LAPACK library"); have_DGGSVD3 = (libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3))) != nullptr); @@ -81,9 +80,14 @@ gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD))); gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD))); } + gsvd_initialized = true; } + /* Clean up macro namespace as soon as we are done using them */ + #undef xSTRINGIZE + #undef STRINGIZE + template<class T1> struct real_ggsvd_ptr { @@ -217,14 +221,14 @@ }; // template specializations + typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type; typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type; - typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type; - typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type; typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type; - typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type; + typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type; typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type; + typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type; + typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type; typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type; - typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type; namespace math { @@ -235,7 +239,7 @@ double *tmp_dataA, F77_INT m1, double *tmp_dataB, F77_INT p1, Matrix& alpha, Matrix& beta, double *u, F77_INT nrow_u, double *v, F77_INT nrow_v, double *q, - F77_INT nrow_q, Matrix& work, F77_INT lwork, + F77_INT nrow_q, double *work, F77_INT lwork, F77_INT *iwork, F77_INT& info) { if (! gsvd_initialized) @@ -250,7 +254,7 @@ 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 (), lwork, iwork, info + work, lwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -264,7 +268,7 @@ 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 + work, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -279,7 +283,7 @@ F77_INT p1, FloatMatrix& alpha, FloatMatrix& beta, float *u, F77_INT nrow_u, float *v, F77_INT nrow_v, float *q, - F77_INT nrow_q, FloatMatrix& work, + F77_INT nrow_q, float *work, F77_INT lwork, F77_INT *iwork, F77_INT& info) { if (! gsvd_initialized) @@ -294,7 +298,7 @@ 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 (), lwork, iwork, info + work, lwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -308,7 +312,7 @@ 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 + work, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -323,13 +327,14 @@ Complex *tmp_dataB, F77_INT p1, Matrix& alpha, Matrix& beta, Complex *u, F77_INT nrow_u, Complex *v, F77_INT nrow_v, Complex *q, - F77_INT nrow_q, ComplexMatrix& work, + F77_INT nrow_q, Complex *work, F77_INT lwork, F77_INT *iwork, F77_INT& info) { if (! gsvd_initialized) initialize_gsvd (); - Matrix rwork(2*n, 1); + OCTAVE_LOCAL_BUFFER(double, rwork, 2*n); + if (have_DGGSVD3) { zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]); @@ -343,8 +348,8 @@ 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 ()), - lwork, rwork.fortran_vec (), iwork, info + F77_DBLE_CMPLX_ARG (work), + lwork, rwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -362,8 +367,8 @@ 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_DBLE_CMPLX_ARG (work), + rwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -381,13 +386,14 @@ FloatComplex *u, F77_INT nrow_u, FloatComplex *v, F77_INT nrow_v, FloatComplex *q, F77_INT nrow_q, - FloatComplexMatrix& work, F77_INT lwork, + FloatComplex *work, F77_INT lwork, F77_INT *iwork, F77_INT& info) { if (! gsvd_initialized) initialize_gsvd (); - FloatMatrix rwork(2*n, 1); + OCTAVE_LOCAL_BUFFER(float, rwork, 2*n); + if (have_DGGSVD3) { cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]); @@ -401,8 +407,8 @@ F77_CMPLX_ARG (u), nrow_u, F77_CMPLX_ARG (v), nrow_v, F77_CMPLX_ARG (q), nrow_q, - F77_CMPLX_ARG (work.fortran_vec ()), lwork, - rwork.fortran_vec (), iwork, info + F77_CMPLX_ARG (work), lwork, + rwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -420,8 +426,8 @@ 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_CMPLX_ARG (work), + rwork, iwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)); @@ -433,13 +439,10 @@ 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; + (*current_liboctave_error_handler) + ("gsvd: U not computed because type == gsvd::sigma_only"); + + return left_smA; } template <typename T> @@ -447,13 +450,10 @@ 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; + (*current_liboctave_error_handler) + ("gsvd: V not computed because type == gsvd::sigma_only"); + + return left_smB; } template <typename T> @@ -461,32 +461,19 @@ 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; - } + (*current_liboctave_error_handler) + ("gsvd: X not computed because type == gsvd::sigma_only"); - 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; + return right_sm; } template <typename T> gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type) { + if (a.isempty () || b.isempty ()) + (*current_liboctave_error_handler) + ("gsvd: A and B cannot be empty matrices"); + F77_INT info; F77_INT m = to_f77_int (a.rows ()); @@ -513,17 +500,16 @@ { 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)]. + // FIXME: In LAPACK 3.0, problem below seems to be fixed, + // so we now set jobX = 'N'. // - // For Lapack 3.0, this problem seems to be fixed. + // For calculating sigma_only, 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)]. - jobu = 'N'; - jobv = 'N'; - jobq = 'N'; + jobu = jobv = jobq = 'N'; nrow_u = nrow_v = nrow_q = 1; break; @@ -533,17 +519,17 @@ type = gsvd_type; - if (! (jobu == 'N' || jobu == 'O')) + if (jobu != 'N') left_smA.resize (nrow_u, m); P *u = left_smA.fortran_vec (); - if (! (jobv == 'N' || jobv == 'O')) + if (jobv != 'N') left_smB.resize (nrow_v, p); P *v = left_smB.fortran_vec (); - if (! (jobq == 'N' || jobq == 'O')) + if (jobq != 'N') right_sm.resize (nrow_q, n); P *q = right_sm.fortran_vec (); @@ -551,7 +537,7 @@ real_matrix alpha (n, 1); real_matrix beta (n, 1); - std::vector<F77_INT> iwork (n); + OCTAVE_LOCAL_BUFFER(F77_INT, iwork, n); if (! gsvd_initialized) initialize_gsvd (); @@ -560,101 +546,148 @@ if (have_DGGSVD3) { lwork = -1; - T work_tmp (1, 1); + P work_tmp; 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_tmp, lwork, - iwork.data (), info); + tmp_dataA, m, tmp_dataB, p, + alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q, + &work_tmp, lwork, iwork, info); - lwork = static_cast<F77_INT> (std::abs (work_tmp(0, 0))); + lwork = static_cast<F77_INT> (std::abs (work_tmp)); } else { - lwork = 3*n; - lwork = (lwork > m ? lwork : m); - lwork = (lwork > p ? lwork : p) + n; + lwork = std::max ({3*n, m, p}) + n; } info = 0; - T work (lwork, 1); + OCTAVE_LOCAL_BUFFER(P, work, lwork); 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, lwork, iwork.data (), - info); + tmp_dataA, m, tmp_dataB, p, + alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q, + work, lwork, iwork, info); if (info < 0) (*current_liboctave_error_handler) - ("*ggsvd.f: argument %" OCTAVE_F77_INT_TYPE_FORMAT " illegal", - -info); - else + ("*ggsvd.f: argument %d illegal", -info); + + if (info > 0) + (*current_liboctave_error_handler) + ("*ggsvd.f: Jacobi-type procedure failed to converge"); + + F77_INT i, j; + + if (gsvd_type != gsvd<T>::Type::sigma_only) { - if (info > 0) - (*current_liboctave_error_handler) - ("*ggsvd.f: Jacobi-type procedure failed to converge."); + // Size according to LAPACK is k+l,k+l, but this needs + // to be nxn for Matlab compatibility. + T R (n, n, 0.0); + int astart = n-k-l; + if (m - k - l >= 0) + { + // 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 { - F77_INT i, j; + // (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 = m; i < k + l; i++) + for (j = n - l - k + m; j < n; j++) + R.xelem (i, j) = btmp.xelem (i - k, 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 ) + // Output X = Q*R' + // FIXME: Is there a way to call BLAS multiply directly + // with flags so that R is transposed? + right_sm = right_sm * R.hermitian (); + } + + // Fill in C and S + F77_INT rank; + bool fill_ptn; + if (m-k-l >= 0) + { + rank = l; + fill_ptn = true; + } + else + { + rank = m-k; + fill_ptn = false; + } - 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 (gsvd_type == gsvd<T>::Type::sigma_only) + { + // Return column vector with results + sigmaA.resize (k+l, 1); + sigmaB.resize (k+l, 1); + + if (fill_ptn) + { + for (i = 0; i < k; i++) + { + sigmaA.xelem (i) = 1.0; + sigmaB.xelem (i) = 0.0; } - - if (m-k-l >= 0) + for (i = k, j = k+l-1; i < k+l; i++, j--) { - // 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); - } + sigmaA.xelem (i) = alpha.xelem (i); + sigmaB.xelem (i) = beta.xelem (i); } - else + } + else + { + for (i = 0; i < k; i++) { - // 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); - } + sigmaA.xelem (i) = 1.0; + sigmaB.xelem (i) = 0.0; + } + for (i = k; i < m; i++) + { + sigmaA.xelem (i) = alpha.xelem (i); + sigmaB.xelem (i) = beta.xelem (i); + } + for (i = m; i < k+l; i++) + { + sigmaA.xelem (i) = 0.0; + sigmaB.xelem (i) = 1.0; } } } + else // returning all matrices + { + // Number of columns according to LAPACK is k+l, but this needs + // to be n for Matlab compatibility. + sigmaA.resize (m, n); + sigmaB.resize (p, n); + + for (i = 0; i < k; i++) + sigmaA.xelem (i,i) = 1.0; + + for (i = 0; i < rank; i++) + { + sigmaA.xelem (k+i,k+i) = alpha.xelem (k+i); + sigmaB.xelem (i,k+i) = beta.xelem (k+i); + } + + if (! fill_ptn) + { + for (i = m; i < n; i++) + sigmaB.xelem (i-k, i) = 1.0; + } + + } } - // Instantiations we need. + // Instantiations needed in octave::math namespace. template class gsvd<Matrix>; template class gsvd<FloatMatrix>; template class gsvd<ComplexMatrix>;
--- a/liboctave/numeric/gsvd.h Thu May 27 22:49:26 2021 +0200 +++ b/liboctave/numeric/gsvd.h Wed Sep 08 16:08:03 2021 -0700 @@ -50,13 +50,13 @@ { } gsvd (const T& a, const T& b, - gsvd::Type gsvd_type = gsvd<T>::Type::economy); + gsvd::Type gsvd_type = gsvd<T>::Type::std); 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) { } + left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm) + { } gsvd& operator = (const gsvd& a) { @@ -68,7 +68,6 @@ left_smA = a.left_smA; left_smB = a.left_smB; right_sm = a.right_sm; - R = a.R; } return *this; @@ -76,26 +75,25 @@ ~gsvd (void) = default; - typename T::real_diag_matrix_type + typename T::real_matrix_type singular_values_A (void) const { return sigmaA; } - typename T::real_diag_matrix_type + typename T::real_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; + real_matrix sigmaA, sigmaB; T left_smA, left_smB; - T right_sm, R; + T right_sm; void ggsvd (char& jobu, char& jobv, char& jobq, octave_f77_int_type m, octave_f77_int_type n, octave_f77_int_type p, @@ -106,7 +104,7 @@ P *u, octave_f77_int_type nrow_u, P *v, octave_f77_int_type nrow_v, P *q, octave_f77_int_type nrow_q, - T& work, octave_f77_int_type lwork, + P *work, octave_f77_int_type lwork, octave_f77_int_type *iwork, octave_f77_int_type& info); };