Error in U matrix calculation when one matrix is rectangular and the other square

[jacklovell]

Problem description

Trying to compute the GSVD of a rectangular and a square matrix results in an error when the U matrix is requested:

import pygsvd
import numpy as np
a = np.random.poisson(1, (4, 15))
b = np.random.poisson(1, (15, 15))
C, S, X, U, V = pygsvd.gsvd(a, b, extras='uv')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jlovell/Software/pygsvd/pygsvd.py", line 116, in gsvd
    U[:, :rank] = U[:, ix]
IndexError: index 14 is out of bounds for axis 1 with size 4

There is no error if only the V matrix is requested:

C, S, X, V = pygsvd.gsvd(a, b, extras='v')
C.shape, S.shape, X.shape, V.shape
((15,), (15,), (15, 15), (15, 15))

There is also no error if both a and b are square. But the error is present if b is rectangular and a is square.

If both a and b are rectangular, a different error occurs:

import pygsvd
import numpy as np
a = np.random.poisson(1, (4, 15))
b = np.random.poisson(1, (6, 15))
C, S, X, U, V = pygsvd.gsvd(a, b, extras='uv')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jlovell/Software/pygsvd/pygsvd.py", line 106, in gsvd
    R = _extract_R(Ac, Bc, k, l)
  File "/home/jlovell/Software/pygsvd/pygsvd.py", line 148, in _extract_R
    R[:m, :] = A
ValueError: could not broadcast input array from shape (4,15) into shape (4,10)

Expected output

The GSVD of these two matrices should be successfully computed, since Octave manages with no problem:

pkg load linear-algebra
a = poissrnd(1, [4, 15]);
b = poissrnd(1, [15, 15]);
[u, v, c, s, x, r] = gsvd(a, b);

Additional details

I'm using Intel MKL 2019, with the mkl 2019.0 and mkl-include 2019.0 packages provided by Anaconda, and Python 3.6.

[bnaecker]

@jacklovell Thanks for filing the issue. I'll take a look, and try to fix it soon. This has to do with how LAPACK returns the computed arrays, with singular values unsorted. The routine tries to arrange them in decreasing order, but I'm clearly doing this incorrectly.