pygsl.linalg.SV_decomp appears to return it's input instead of U when it's input has been transposed. The below script shows the problem.
#!/usr/bin/env python
import numpy as np
import pygsl
import random
pygsl.import_all()
def randomSymmetricMatrix(n):
foo = np.matrix(np.zeros((n,n),dtype='float64'))
for i in range(0,n):
foo[(i,i)] = 1
for j in range(i+1,n):
r = random.random()
foo[(i,j)] = r
foo[(j,i)] = r
return foo
foo = randomSymmetricMatrix(2)
oldFoo = foo.copy()
print "Is foo equal to transpose", foo == foo.transpose()
print foo.trace()
print
def test(tag, bar):
U,V,s = pygsl.linalg.SV_decomp(bar)
print tag,"strides ",bar.strides
print tag,"original",foo
print tag,"U ",U==foo,U
print tag,"V ",V
print tag,"s ",s
print
bar = np.array(foo,dtype='float64')
test("Good",bar)
bar = (np.array(foo,dtype='float64').transpose())
test("Bad ",bar)
bar = (np.array(foo.transpose(),dtype='float64'))
test("Bad2",bar)
bar = (np.matrix(foo.transpose(),dtype='float64'))
test("Bad3",bar)
bar = (np.matrix(foo.transpose().transpose(),dtype='float64'))
test("Good",bar)
bar = (np.matrix(foo.transpose(),dtype='float64').transpose())
test("Good",bar)
bar = (np.matrix(foo.transpose().copy(),dtype='float64'))
test("Good",bar)
bar = (foo.transpose())
test("Bad4",bar)
print foo.trace()
print foo == oldFoo
print