[SciPy-dev] lu_solve problem -- f2py related

Travis Oliphant oliphant.travis at ieee.org
Wed Apr 3 23:38:49 CST 2002


Pearu,


I discovered an interesting behavior in f2py's wrapping of clapack_dgetrs.  
This is the function that solves a system of equations.

It results in wrong answers when the argument b is a one-dimensional array 
that is not contiguous and the functions clapack_dgetrs (or clapack_gesv) are 
called.  This may have implications for other functions as well, I'm not sure.

The easiest way to show the problem is to look at linalg2.solve for the 
following problem.

A = rand(5,5)
b1 = rand(5,2)

x1 = linalg2.solve(A,b1)
x2 = linalg2.solve(A,b1[:,0])
x3 = linalg2.solve(A,b1[:,0].copy())

print x1[:,0]
print x2
print x3


Notice that x2 is not correct.  The problem I think is the fact that in 
array_from_pyobj there is a check for (rank > 1) arrays before 
lazy-transposes are done.  I think a check should be made for rank-1 arrays 
as well, to handle this case.

For, now I've just made a hack in the solve routines, to check for this case.

I'd love to hear your perspective on this problem, though, Pearu.

-Travis



More information about the Scipy-dev mailing list