sparse matrices

[dschouten]

Hi,

Many thanks for wrapping the GSVD algorithm in Python.

I am using LAPACK 3.4.2 (CentOS 7) so I had to modify the .c file slightly to use DGGSVD instead of DGGSVD3 (as in LAPACK 3.8.0)

But when I run it for the matrices I am interested int, I get an error I do not understand.

My code is below.

import numpy as np
import scipy
import scipy.sparse
import scipy.sparse.linalg
import pygsvd

m = 50
n = 100
p = n-1

# A[m x n] is a random sparse matrix 
A = scipy.sparse.random(m, n, density=5e-2, dtype=float, format='csr')

# L[p x n], p = n-1 looks like
#
#   1 -1  0  0  0  0  0 ... 0 
#   0  1 -1  0  0  0  0 ... 0 
#   0  0  1 -1  0  0  0 ... 0 
#            .                  
#                 .              
#                      .         
#   0  0  0  ...   0  0  1 -1 
# 

L = scipy.sparse.csr_matrix(
        (np.concatenate((np.ones(shape=(p,)), -np.ones(shape=(p,)))),
        (np.concatenate((np.arange(0, p), np.arange(0, p  ))),
         np.concatenate((np.arange(0, p), np.arange(1, p+1))))), shape=(p, n))

C,S,X,U,V = pygsvd.gsvd(L.todense(), A.todense(), full_matrices=True) 

The error I get is

L116 IndexError: index 99 is out of bounds for axis 1 with size 99

Any ideas what is going on?

I should add that when I do

C,S,X,U,V = pygsvd.gsvd(L.todense(), L.todense(), full_matrices=True) 

or

C,S,X,U,V = pygsvd.gsvd(A.todense(), A.todense(), full_matrices=True) 

then I get the desired result.