[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