[SciPy-dev] problem with linalg.cholesky?

[Andrew Jaffe]

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