Sci.linalg.lu permuation error

[Luke]

I'm trying to use Scipy's LU factorization. Here is what I've got:

from numpy import *
import scipy as Sci
import scipy.linalg

A=array([[3., -2., 1., 0., 0.],[-1., 1., 0., 1., 0.],[4., 1., 0., 0.,
1.]])
p,l,u=Sci.linalg.lu(A,permute_l = 0)

now, according to the documentation:
**************************************************
Definition: Sci.linalg.lu(a, permute_l=0, overwrite_a=0)
Docstring:
Return LU decompostion of a matrix.

Inputs:

a -- An M x N matrix.
permute_l -- Perform matrix multiplication p * l [disabled].

Outputs:

p,l,u -- LU decomposition matrices of a [permute_l=0]
pl,u -- LU decomposition matrices of a [permute_l=1]

Definitions:

a = p * l * u [permute_l=0]
a = pl * u [permute_l=1]

p - An M x M permutation matrix
l - An M x K lower triangular or trapezoidal matrix
with unit-diagonal
u - An K x N upper triangular or trapezoidal matrix
K = min(M,N)
*********************************************


So it would seem that from my above commands, I should have:
A == dot(p,dot(l,u))

but instead, this results in:
In [21]: dot(p,dot(l,u))
Out[21]:
array([[-1., 1., 0., 1., 0.],
[ 4., 1., 0., 0., 1.],
[ 3., -2., 1., 0., 0.]])

Which isn't equal to A!!!
In [23]: A
Out[23]:
array([[ 3., -2., 1., 0., 0.],
[-1., 1., 0., 1., 0.],
[ 4., 1., 0., 0., 1.]])

However, if I do use p.T instead of p:
In [22]: dot(p.T,dot(l,u))
Out[22]:
array([[ 3., -2., 1., 0., 0.],
[-1., 1., 0., 1., 0.],
[ 4., 1., 0., 0., 1.]])

I get the correct answer.

Either the documentation is wrong, or somehow Scipy is returning the
wrong permutation matrix... anybody have any experience with this or
tell me how to submit a bug report?

Thanks.
~Luke

 

[bernhard.voigt]

Hi Luke,

you should send this to the scipy user list: (e-mail address removed)

Bernhard