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 -- http://mail.python.org/mailman/listinfo/python-list