[SciPy-dev] problem with linalg.cholesky?
Andrew Jaffe
a.h.jaffe at gmail.com
Fri Dec 2 15:31:09 CST 2005
Travis Oliphant wrote:
> Andrew Jaffe wrote:
>>
>>OK, drilling down a bit further... The problem lies in the use of
>>scipy.transpose() in _castCopyAndTranspose, since transpose() makes a
>>contiguous, non-Fortran array into a non-contiguous fortran array (since
>>it appears that at present it doesn't copy, just makes a new view). One
>>solution is just to actually copy the array. The other function,
>>_fastCopyAndTranspose(), seems to work fine in this situation, but I'm
>>not sure how and why this differs from _castCopyAndTranspose (in
>>particular, it doesn't seem to set the Fortran flag).
>>
>>[However, the cholesky decomposition only makes sense for a symmetric
>>(or Hermitian) array, so, if we're not changing the type, for real
>>matrices really all we need to do is the cast and copy; for complex
>>matrices we need to transpose (not hermitian conjugate since we really
>>want to change the ordering), but this is equivalent to just taking the
>>complex conjugate of every element, which must be faster than copying.
>>If we are changing the type, I assume the speed difference isn't as
>>important.]
>>
>>Any ideas on what the best change(s) would be?
>>
> The problem with making changes based on what's been posted so far is
> that the code works for other platforms. If this is a problem with OS
> X, is this a problem for just your installation or does everybody else
> see it.
>
> I have not seen enough data to know. One strong possibility is that
> the error is in the lapack_lite and most people are using the optimized
> lapack. We should definitely check that direction.
I still do get the problem with the latest version (rev 1550). I thought
I was using the optimized lapack -- it's automatic on OS X, isn't it?
Does the fact that the problem is occuring during a lapack_lite call
imply something else?
If there's no problem there, I guess the questions are: does
_castCopyAndTranspose result in a non-contiguous array on all platforms?
And does lapack_lite.dpotrf() balk at that, and why?
Andrew
More information about the Scipy-dev
mailing list