[SciPy-user] problems with scipy.test() on Opteron 64-bit

[Pearu Peterson]

On Thu, 25 Nov 2004, Giovanni Samaey wrote:

> Pearu Peterson wrote:
>
>> 
>> 
>> On Thu, 25 Nov 2004, Giovanni Samaey wrote:
>> 
>>> What is remarkable is that, the error is the same in the 3 cases. 
>>> Moreover, if you add up the first two components of array1 you get the 
>>> first component of array2.
>> 
>> 
>> It would look like atlas bug but since the same results are obtained with 
>> Fortran blas/lapack, then I don't know what is happening. I presume that
>> you did
>>   rm -rf build
>
> Of course. If I would forget, I would notice that nothing would be built, 
> anyway...
>
>> What happens when you call geev routine directly? Try for example:
>> 
>> In [1]: import scipy
>> 
>> In [2]: scipy.linalg.flapack.dgeev([[1,2,3],[1,2,3],[2,5,6]],0,0)
>> Out[2]:
>> (array([  9.32182538e+00,  -6.03426271e-16,  -3.21825380e-01]),
>>  array([ 0.,  0.,  0.]),
>>  array([       [ 0.,  0.,  0.]]),
>>  array([       [ 0.,  0.,  0.]]),
>>  0)
>> 
>> In [3]: scipy.linalg.flapack.sgeev([[1,2,3],[1,2,3],[2,5,6]],0,0)
>> Out[3]:
>> (array([  9.32182503e+00,  -2.82419336e-07,  -3.21825117e-01],'f'),
>>  array([ 0.,  0.,  0.],'f'),
>>  array([       [ 0.,  0.,  0.]],'f'),
>>  array([       [ 0.,  0.,  0.]],'f'),
>>  0)
>> 
>> In [4]: scipy.linalg.calc_lwork.geev('d',3)
>> Out[4]: (12, 102)

This was obtained on a 32-bit system.

>> Is there any difference in your case?
>
> Yes there is; in fact there is a difference in all three.
>
> In[1]: import scipy
> In[2]: scipy.linalg.flapack.dgeev([[1,2,3],[1,2,3],[2,5,6]],0,0)
> Out[2]: (array([ 9.43719064, -0.11536526, -0.32182538]),
> array([ 0.,  0.,  0.]),
> array([       [ 0.,  0.,  0.]]),
> array([       [ 0.,  0.,  0.]]),
> 0)

Hmm, try also

scipy.linalg.flapack.zgeev([[1,2,3],[1,2,3],[2,5,6]],0,0)
scipy.linalg.flapack.dgeev(scipy.array([[1,2,3],[1,2,3],[2,5,6]],'d'),0,0)
scipy.linalg.flapack.dgeev(Numeric.array([[1,2,3],[1,2,3],[2,5,6]],'d'),0,0)

If those fail as well then try a C program that uses dgeev to solve the 
same eigenvalue problem:

/* main.c */
#include <stdio.h>
extern 
dgeev_(char*,char*,int*,double*,int*,double*,double*,double*,int*,double*,int*,double*,int*,int*);
main() {
   int n=3,lwork=12,info=0;
   double a[] = {1,1,2,2,2,5,3,3,6};
   double wr[]={0,0,0},wi[]={0,0,0},*vl,*vr,work[lwork];
   dgeev_("N","N",&n,a,&n,wr,wi,vl,&n,vr,&n,work,&lwork,&info);
   printf("wr=%g,%g,%g,wi=%g,%g,%g\n",wr[0],wr[1],wr[2],wi[0],wi[1],wi[2]);
}
/*eof*/

$ gcc main.c -llapack
$ ./a.out
wr=9.32183,-6.20979e-16,-0.321825,wi=0,0,0


> In[3]: scipy.linalg.flapack.sgeev([[1,2,3],[1,2,3],[2,5,6]],0,0)
> Out[3]: (array([ 0.,  0.,  0.],'f'),
> array([ 0.,  0.,  0.],'f'),
> array([       [ 0.,  0.,  0.]],'f'),
> array([       [ 0.,  0.,  0.]],'f'), 3)

Here info==3 which means that the QR algorithm  failed  to  compute  all 
the eigenvalues and 0:info of eigenvalues are such that converged. But
that does not makes sense because info==n...

> In[4]: scipy.linalg.calc_lwork.geev('d',3)
> Out[4]: (12, 174)

I get the same result on an Opteron box.

Pearu