# [Scipy-tickets] [SciPy] #1458: lu_factor gives different results, depending on whether the input is C or Fortran contiguous (also affects lu_solve)

SciPy Trac scipy-tickets@scipy....
Thu Jun 9 08:52:07 CDT 2011

```#1458: lu_factor gives different results, depending on whether the input is C or
Fortran contiguous (also affects lu_solve)
--------------------------+-------------------------------------------------
Reporter:  mike.wimmer   |       Owner:  somebody
Type:  defect        |      Status:  new
Priority:  normal        |   Milestone:
Component:  scipy.linalg  |     Version:  0.9.0
Keywords:                |
--------------------------+-------------------------------------------------
scipy.inalg.lu_factor returns different results when called with a C or
Fortran contiguous array. The behavior can easily be seen for the
following small example:

import numpy as np[[BR]]
import scipy.linalg as la[[BR]]

A=np.array([[2,1],[0,1]])[[BR]]
la.lu_factor(A)[[BR]]
la.lu_factor(np.asfortranarray(A))

which for me yields the output

(array([[ 2. ,  0.5],[[BR]]
[ 0. ,  1. ]]), array([0, 1], dtype=int32))
(array([[ 2.,  1.],[[BR]]
[ 0.,  1.]]), array([0, 1], dtype=int32))

The matrices lu are obviously different, and only the second output (the
one from la.lu_factor(np.asfortranarray(A))) are of the form as described
by the documentation of lu_factor, with unit triangular L in the strictly
lower and U in the upper triangle of lu.

I suspect that lu_factor acutally factors the transposed system when a C
contiguous array is used.

In fact, the error cancels out when lu_solve is called:

lu = la.lu_factor(A)[[BR]]
la.lu_solve(lu, [1,2])[[BR]]
lu = la.lu_factor(np.asfortranarray(A))[[BR]]
la.lu_solve(lu, [1,2])[[BR]]

give identical results:

array([-0.5,  2. ])[[BR]]
array([-0.5,  2. ])

I suspect again this is because lu_solve actually solves the transposed
system, if it gets as an input a C contiguous LU factorization. It is
however easy to mess up lu_solve, too:

lu = la.lu_factor(A)[[BR]]
la.lu_solve((np.asfortranarray(lu[0]), lu[1]), [1,2])

gives the wrong result

array([ 0.,  2.])

It seems to me that lu_factor and lu_solve use the way the array is stored
as a flag in whether to internally transpoe the system or not.

The function lu() is not affected by this problem.

In my opinion this behavior is not appropriate: Results shouldn't depend
on what kind of internal memory layout one has.

I see two ways to fix this: Either always return the same result - this
will require internally converting C to Fortran contiguous. Or change the
documentation of lu_factor and lu_solve to indicate that the return value
of lu_factor may have different meanings depending on the memory layout of
the input, and that the return value of lu_factor should not be touched
before passing it to lu_solve.

--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1458>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
```

More information about the Scipy-tickets mailing list