scipy.linalg.lapack.dgesjv overwrites original arrays if called on a transposed array

[highlando]

The wrapper for dgejsv in scipy.linalg.lapack (cp. #11354) seems to have troubles with the transpose (even if deep copies are used). If dgejsv is called on an array that had been transposed, it overwrites it. This does not happen for the the untransposed.

The computed values are not affected.

Reproducing code example:

import numpy as np
from scipy.linalg.lapack import dgejsv

n, m = 10, 10
amat = np.sin(np.arange(m*n).reshape(m, n))
amat = amat + amat.T  # make it symmetric

print('norm of A: ', np.linalg.norm(amat))
sva, u, v, wo, iwo, info = dgejsv(amat)
print('norm of A: ', np.linalg.norm(amat))
svad = np.diag(sva)
print('SVD error: ', np.linalg.norm(amat - u.dot(svad.dot(v.T))))

amatt = np.copy(amat.T)
bmatt = np.copy(amatt)

print('norm of A.T: ', np.linalg.norm(amatt))
sva, u, v, wo, iwo, info = dgejsv(amatt)
print('norm of A.T: ', np.linalg.norm(amatt))
svad = np.diag(sva)
print('SVD error: ', np.linalg.norm(bmatt - u.dot(svad.dot(v.T))))

Error message:

Output

norm of A:  10.035024236333806
norm of A:  10.035024236333806
SVD error:  4.5248923239679534e-15
norm of A.T:  10.035024236333806
norm of A.T:  1.2350051981453427e+154
SVD error:  4.5248923239679534e-15

Scipy/Numpy/Python version information:

python3 -c "import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)"
1.5.2 1.19.4 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)

[ilayn]

Transpose is returning a view and copy is still taking the original array. If you replace the line with

amatt = amat.T.copy()

then probably it should work but I'll check it properly later when I have some time in case others didn't already.

[peterbell10]

To be consistent with the note at the top of the docs, shouldn't this function have an overwrite_a parameter that defaults to False?

[ilayn]

ouch, if it is not there, then we forgot to add it in the wrapper.

[peterbell10]

Doesn't seem to be there:

>>> lapack.dgejsv(amat, overwrite_a=False)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-29-c2555ce2437c> in <module>
----> 1 lapack.dgejsv(amat, overwrite_a=False)

TypeError: 'overwrite_a' is an invalid keyword argument for _flapack.dgejsv()

[highlando]

Another hint was given to me by @pmli

[...] I think the issue is in C vs. Fortran ordering. Quoting from scipy.linalg.lapack docs:

Similarly, if a C-contiguous array is passed, f2py will pass a FORTRAN-contiguous array internally. Please make sure that these details are satisfied. More information can be found in the f2py documentation.

Also, this might be related to double vs. single precision (as mentioned ibid. in the second paragraph). For arrays that only had integer values, this issue did not occur.

[highlando]

Transpose is returning a view and copy is still taking the original array. If you replace the line with

amatt = amat.T.copy()

then probably it should work but I'll check it properly later when I have some time in case others didn't already.

This reformulation, in fact, avoids this issue.

[ilayn]

This is an issue in the wrapper. We usually add a copying option and f2py creates an overwrite_a option automatically being False as the default. Here it is True by default and gets overwritten because transpose of a C-contiguous array is Fortran contiguous and gets overwritten.

@tylerjereddy I'm going to fix this asap. Could you please tag accordingly whether we can squeeze it into 1.6 otherwise 1.7 is fine.