[SciPy-dev] lapack and blas

[Pearu Peterson]

Hi Eric,

On Sat, 19 Jan 2002, eric wrote:

> > However, a bad news is that lapack and atlas lapack functions are not
> > compatible. 
> 
> > For example,
> > >>> LU,ipiv,X,info=flapack.c_sgesv([[1,2],[3,4]],[1,1])
> > >>> print X
> > [-0.5  0.5]
> 
> For some reason, this is solving using the transpose of the matrix. 
> 
>    1   3  x1 = 1  -->   x1 = -0.5
>    2   4  x2   1        x2 =  0.5
> 
> Is this a problem with c_sgesv, or is f2py doing something incorrectly
> with the matrix transpose stuff?  I need to look at the magic your doing
> under the covers to handle all this, since I haven't figured out how
> you'd do it completely transparently in a way that never leads to
> ambiguity and also maintains efficiency.     

I'd appreciate if you could look at the magic that f2py generates. I
welcome any critisism and I am happy to explain the reasons why I did this
or that. You may find there something that I have missed...

However, this issue should be solved in the latest clapack and
flapack interfaces (the problem was related to sgesv handeling the rhs
matrix in the column-major storage order).

Here is an example:

>>> clapack.dgesv([[1,2],[3,4]],[[1,2,1],[2,0,2]])[2]
array([[ 0. , -4. ,  0. ],
       [ 0.5,  3. ,  0.5]])
>>> flapack.dgesv([[1,2],[3,4]],[[1,2,1],[2,0,2]])[2]
array([[ 0. , -4. ,  0. ],
       [ 0.5,  3. ,  0.5]])

that show the solution matrix X for

 [ 1  2 ]         [ 1  2  1 ]
 |      |  *  X = |         |
 [ 3  4 ]         [ 2  0  2 ]

calculated using atlas C lapack and Fortran lapack, respectively.

Regards,
	Pearu